-
Notifications
You must be signed in to change notification settings - Fork 0
/
write_norm_model.f90
176 lines (138 loc) · 4.58 KB
/
write_norm_model.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
program write_norm
implicit none
integer, parameter :: CUSTOM_REAL = 4
integer, parameter :: IIN = 4
integer :: NITER, nline, i, j
real(kind=CUSTOM_REAL), allocatable, dimension(:) :: target_rho, target_vp, target_vs, &
model_rho, model_vp, model_vs, error_rho, error_vp, error_vs
double precision, allocatable, dimension(:) :: error_rho_L1, error_vp_L1, &
error_vs_L1, error_rho_L2, error_vp_L2, error_vs_L2
character(len=150) :: dir, filename, arg
integer :: ios, filesize
character(len=150) :: suffix
! read input arguments
if( iargc() .ne. 2 ) stop 'USAGE: write_misfit NITER DIR'
j=1; call getarg(j, arg); read(arg,*,iostat=ios) NITER;
if (ios /= 0) stop 'Error reading NITER'
j=2; call getarg(j, dir);
if (ios /= 0) stop 'Error reading DIR'
dir = trim(dir)
! dynamically allocate arrays
filename = trim(dir) // '/output/model_true/proc000000_rho.bin'
open(unit=3,file=filename,status='old',action='read')
inquire(3, size=filesize)
close(3)
nline=(filesize-8)/4
allocate( target_rho(nline) )
allocate( target_vp(nline) )
allocate( target_vs(nline) )
allocate( model_rho(nline) )
allocate( model_vp(nline) )
allocate( model_vs(nline) )
allocate( error_rho(nline) )
allocate( error_vp(nline) )
allocate( error_vs(nline) )
allocate( error_rho_L1(0:NITER) )
allocate( error_vp_L1(0:NITER) )
allocate( error_vs_L1(0:NITER) )
allocate( error_rho_L2(0:NITER) )
allocate( error_vp_L2(0:NITER) )
allocate( error_vs_L2(0:NITER) )
target_rho(:) = 0.0d0
target_vp(:) = 0.0d0
target_vs(:) = 0.0d0
model_rho(:) = 0.0d0
model_vp(:) = 0.0d0
model_vs(:) = 0.0d0
error_rho(:) = 0.0d0
error_vp(:) = 0.0d0
error_vs(:) = 0.0d0
error_rho_L1(:) = 0.0d0
error_vp_L1(:) = 0.0d0
error_vs_L1(:) = 0.0d0
error_rho_L2(:) = 0.0d0
error_vp_L2(:) = 0.0d0
error_vs_L2(:) = 0.0d0
! read target model
filename = trim(dir) // '/output/model_true/proc000000_rho.bin'
open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
read(IIN) target_rho
close(IIN)
filename = trim(dir) // '/output/model_true/proc000000_vp.bin'
open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
read(IIN) target_vp
close(IIN)
filename = trim(dir) // '/output/model_true/proc000000_vs.bin'
open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
read(IIN) target_vs
close(IIN)
! read subsequent models
do i = 0, NITER
if (i==0) then
suffix = 'init'
else
write(suffix,'(i4.4)') i
endif
!filename =trim(dir) // '/output/model_'//trim(suffix) // '/proc000000_rho.bin'
!open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
!read(IIN) model_rho
!close(IIN)
filename =trim(dir) // '/output/model_'//trim(suffix) //'/proc000000_vp.bin'
open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
read(IIN) model_vp
close(IIN)
filename =trim(dir) // '/output/model_'//trim(suffix) //'/proc000000_vs.bin'
open(unit=IIN,file=filename,status='old',action='read',form='unformatted')
read(IIN) model_vs
close(IIN)
!error_rho(:)=target_rho(:)-model_rho(:)
error_vp(:)=target_vp(:)-model_vp(:)
error_vs(:)=target_vs(:)-model_vs(:)
!error_rho_L1(i) = SUM(abs(error_rho(:)))/nline
error_vp_L1(i) = SUM(abs(error_vp(:)))/nline
error_vs_L1(i) = SUM(abs(error_vs(:)))/nline
!error_rho_L2(i) = (DOT_PRODUCT(error_rho(:),error_rho(:))/nline)**0.5
error_vp_L2(i) = (DOT_PRODUCT(error_vp(:),error_vp(:))/nline)**0.5
error_vs_L2(i) = (DOT_PRODUCT(error_vs(:),error_vs(:))/nline)**0.5
print *, i, &
!error_rho_L2(i),&
error_vp_L2(i),&
error_vs_L2(i)
enddo
!filename = trim(dir) // '/output.stats/error_rho_L1'
!open(unit=6,file=filename,status='unknown',action='write')
!do i = 0, NITER
! write(6,*) error_rho_L1(i)
!enddo
!close(6)
!filename = trim(dir) // '/output.stats/error_vp_L1'
!open(unit=6,file=filename,status='unknown',action='write')
!do i = 0, NITER
write(6,*) error_vp_L1(i)
!enddo
!close(6)
!filename = trim(dir) // '/output.stats/error_vs_L1'
!open(unit=6,file=filename,status='unknown',action='write')
!do i = 0, NITER
! write(6,*) error_vp_L1(i)
!enddo
!close(6)
filename = trim(dir) // '/output.stats/error_rho_L2'
open(unit=6,file=filename,status='unknown',action='write')
do i = 0, NITER
write(6,*) error_rho_L2(i)
enddo
close(6)
filename = trim(dir) // '/output.stats/error_vp_L2'
open(unit=6,file=filename,status='unknown',action='write')
do i = 0, NITER
write(6,*) error_vp_L2(i)
enddo
close(6)
filename = trim(dir) // '/output.stats/error_vs_L2'
open(unit=6,file=filename,status='unknown',action='write')
do i = 0, NITER
write(6,*) error_vs_L2(i)
enddo
close(6)
end program