-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfixedSlovecsEventSolve.f
300 lines (300 loc) · 9.09 KB
/
fixedSlovecsEventSolve.f
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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
C
C Steven J Gibbons
C NGI
C
C fixed slowness vectors event solve
C
C Takes 3 input arguments NEV NSP chrefev
C
PROGRAM fixedSlovecsEventSolve
IMPLICIT NONE
C
INTEGER NIARGS
INTEGER IARGC
INTEGER IARG
INTEGER IFIXLE
C
INTEGER NRELVL
INTEGER NRELVM
PARAMETER ( NRELVM = 50000 )
INTEGER NSP
INTEGER NSPMAX
PARAMETER ( NSPMAX = 120 )
INTEGER NEV
INTEGER NEVMAX
PARAMETER ( NEVMAX = 60 )
C
CHARACTER*(200) CHARG
CHARACTER*(80) CMESS
C
CHARACTER*(8) CSTARR( NSPMAX )
INTEGER LSTARR( NSPMAX )
CHARACTER*(8) CPHARR( NSPMAX )
INTEGER LPHARR( NSPMAX )
REAL*8 DSXARR( NSPMAX )
REAL*8 DSYARR( NSPMAX )
REAL*8 DSTLAT( NSPMAX )
REAL*8 DSTLON( NSPMAX )
REAL*8 DRFLAT( NSPMAX )
REAL*8 DRFLON( NSPMAX )
C
CHARACTER*(24) CEVARR( NEVMAX )
INTEGER LEVARR( NEVMAX )
REAL*8 DRXARR( NEVMAX )
REAL*8 DRYARR( NEVMAX )
C
INTEGER IE1ARR( NRELVM )
INTEGER IE2ARR( NRELVM )
REAL*8 DE1ARR( NRELVM )
REAL*8 DE2ARR( NRELVM )
INTEGER ISPARR( NRELVM )
REAL*8 DCCARR( NRELVM )
INTEGER NUMOU( NRELVM )
REAL*8 DCRES( NRELVM )
C
INTEGER NUMEVO( NEVMAX )
INTEGER NUMSPO( NSPMAX )
INTEGER NLEV
INTEGER NLSP
INTEGER INDLEV( NEVMAX )
INTEGER INDLSP( NSPMAX )
C
C LDA is the maximum number of observation pairs
C
INTEGER LDA
PARAMETER ( LDA = 20000 )
INTEGER ND
C
INTEGER LWORK
PARAMETER ( LWORK = 2*NRELVM )
INTEGER LWOPT
INTEGER NM
PARAMETER ( NM = 3 )
REAL*8 DWORK1( LWORK )
REAL*8 DAMAT( LDA, NM )
REAL*8 DDVEC( LDA )
REAL*8 DWEIG( LDA )
REAL*8 DRESV( LDA )
REAL*8 DOLDRV( LDA )
REAL*8 DTEMP1( LDA )
INTEGER IOBS( LDA )
INTEGER JOBS( LDA )
C
INTEGER NPAIRM
PARAMETER ( NPAIRM = NEVMAX*(NEVMAX-1) )
INTEGER LDATCM
PARAMETER ( LDATCM = NEVMAX )
REAL*8 DATCVM( LDATCM, LDATCM )
INTEGER IABSVL( NPAIRM )
INTEGER JABSVL( NPAIRM )
REAL*8 DRELVX( NPAIRM )
REAL*8 DRELVY( NPAIRM )
REAL*8 DABSVX( NEVMAX )
REAL*8 DABSVY( NEVMAX )
REAL*8 DWGHTA( NPAIRM )
REAL*8 DRELRX( NPAIRM )
REAL*8 DRELRY( NPAIRM )
REAL*8 DWORK2( NPAIRM )
REAL*8 DOLDR2( NPAIRM )
REAL*8 DCVMEX( NEVMAX )
REAL*8 DCVMEY( NEVMAX )
C
REAL*8 DRESPE( NEVMAX )
REAL*8 DRESPS( NSPMAX )
INTEGER NRESPE( NEVMAX )
INTEGER NRESPS( NSPMAX )
C
REAL*8 DCOVEL
C
INTEGER I
INTEGER IEV
INTEGER ITER
C
C CTARER is the name of the reference event
C LTARER is the length of the character string CTARER
CHARACTER*(24) CTARER
INTEGER LTARER
C
INTEGER IERR
INTEGER LUIN
C
CMESS = ' '
NIARGS = IARGC()
IF ( NIARGS.NE.3 ) THEN
WRITE (6,*) 'Usage: NSP NEV EVENTA '
WRITE (6,*) ' 111 6 DPRK2 '
CALL EXIT(1)
ENDIF
C
CHARG = ' '
IARG = 1
CALL GETARG( IARG, CHARG )
CMESS = 'Error reading integer NSP'
READ ( CHARG, *, ERR=99, END=99 ) NSP
IF ( NSP.LT.1 .OR. NSP.GT.NSPMAX ) THEN
WRITE (6,*) 'NSP = ', NSP
WRITE (6,*) 'NSPMAX = ', NSPMAX
CMESS = 'NSPMAX exceeded'
GOTO 99
ENDIF
C
CHARG = ' '
IARG = 2
CALL GETARG( IARG, CHARG )
CMESS = 'Error reading integer NEV'
READ ( CHARG, *, ERR=99, END=99 ) NEV
IF ( NEV.LT.1 .OR. NEV.GT.NEVMAX ) THEN
WRITE (6,*) 'NEV = ', NEV
WRITE (6,*) 'NEVMAX = ', NEVMAX
CMESS = 'NEVMAX exceeded'
GOTO 99
ENDIF
C
CHARG = ' '
IARG = 3
CALL GETARG( IARG, CHARG )
LTARER = 24
DO I = 2, 24
IF ( CHARG(I:I).EQ.' ' .AND. LTARER.EQ.24 ) LTARER = I-1
ENDDO
CTARER = ' '
CTARER(1:LTARER) = CHARG(1:LTARER)
C
print *,' nsp = ', nsp
print *,' nev = ', nev
print *,' ltarer = ', ltarer
print *,' ctarer = ', ctarer
C
C Now read those lines of the input file for the NSP slowness vectors
C
LUIN = 5
CALL NSPSLR( IERR, LUIN, NSPMAX, NSP,
1 CSTARR, LSTARR, CPHARR, LPHARR,
2 DSTLAT, DSTLON, DRFLAT, DRFLON,
3 DSXARR, DSYARR )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'NSPSLR returned IERR = ', IERR
CMESS = 'Error from NSPSLR '
GOTO 99
ENDIF
print *,' NSP slowness vectors read.'
C
C Now read those lines of the input file for the NEV events
C
CALL NEVELR( IERR, LUIN, NEVMAX, NEV,
1 CEVARR, LEVARR, DRXARR, DRYARR )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'NEVELR returned IERR = ', IERR
CMESS = 'Error from NEVELR '
GOTO 99
ENDIF
print *,' NEV events read.'
C
C Now read the remaining lines of the input file for the NRELVL
C observations.
C
CALL DTLNRD( IERR, LUIN, NEV, NSP, NRELVL, NRELVM,
1 CEVARR, LEVARR, CSTARR, LSTARR, CPHARR, LPHARR,
2 IE1ARR, IE2ARR, DE1ARR, DE2ARR, ISPARR, DCCARR)
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'DTLNRD returned IERR = ', IERR
CMESS = 'Error from DTLNRD '
GOTO 99
ENDIF
print *,' NRELVL observations read. NRELVL = ', NRELVL
C
C Initialize the arrays NUMOU and DCRES
C
DO I = 1, NRELVL
NUMOU( I ) = 0
DCRES( I ) = 0.0d0
ENDDO
C
C Now we have read in all of our observations, we need to
C calculate which of our events and slowness vectors are "live"
C i.e. we may have read in some of the NEV and NSP that do not
C have (sufficient) observations.
C So we calculate NLEV and NLSP - the numbers of those that do.
C
CALL NLSPEF( IERR, NEV, NSP, NRELVL,
1 IE1ARR, IE2ARR, ISPARR, NUMEVO, NUMSPO,
2 NLEV, NLSP, INDLEV, INDLSP )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'NLSPEF returned IERR = ', IERR
CMESS = 'Error from NLSPEF '
GOTO 99
ENDIF
print *,' NLEV = ', NLEV
print *,' NLSP = ', NLSP
C
C Now we have to calculate the numbers of the "live" events
C that CTARER and CTAREB correspond to.
C If they are not found and not distinct then we have to abort.
C
CALL INDLEF( IERR, IFIXLE, LTARER, CTARER, NLEV, NEV, NEVMAX,
1 INDLEV, CEVARR, LEVARR )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'INDLEF returned IERR = ', IERR
CMESS = 'Error from INDLEF '
GOTO 99
ENDIF
print *,' IFIXLE = ', IFIXLE
C
ITER = 1
CALL EABSLF( IERR, NLEV, IFIXLE, NEV, INDLEV, NSP, NPAIRM,
1 ITER, NRELVL, IE1ARR, IE2ARR, ISPARR, DE1ARR, DE2ARR,
2 DCCARR, DSXARR, DSYARR, LDA, ND, DAMAT, DDVEC,
3 DWEIG, IOBS, JOBS, LWORK, LWOPT, DWORK1, DRESV,
4 DOLDRV, DTEMP1, IABSVL, JABSVL, DRELVX, DRELVY,
5 DABSVX, DABSVY, DWGHTA, DRELRX, DRELRY,
6 LDATCM, DATCVM, DWORK2, DOLDR2, NUMOU, DCRES,
7 DCVMEX, DCVMEY )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'EABSLF returned IERR = ', IERR
CMESS = 'Error from EABSLF '
GOTO 99
ENDIF
DO I = 1, NRELVL
IF ( NUMOU(I).GT.0 ) THEN
WRITE (6,82) I, NUMOU(I), DCRES(I),
1 DCRES(I)/DBLE( NUMOU(I) )
ENDIF
ENDDO
82 FORMAT('Observation ',I6,' uses ',I8,' cres ',f20.4,
1 ' avg ',f20.4)
DO I = 1, NLEV
IEV = INDLEV( I )
DRXARR( IEV ) = DABSVX( I )
DRYARR( IEV ) = DABSVY( I )
ENDDO
CALL CALCRA( IERR, NLEV, NEV,
1 INDLEV, NSP, NRELVL, IE1ARR, IE2ARR, ISPARR,
2 DE1ARR, DE2ARR, DSXARR, DSYARR,
3 DRXARR, DRYARR, DRESPE, DRESPS,
4 NRESPE, NRESPS )
IF ( IERR.NE.0 ) THEN
WRITE (6,*) 'CALCRA returned IERR = ', IERR
CMESS = 'Error from CALCRA '
GOTO 99
ENDIF
C
DO I = 1, NLEV
IEV = INDLEV( I )
DCOVEL = DRESPE( IEV )/NRESPE( IEV )
c DCOVEL = DSQRT( DCVMEX( I )*DCVMEX( I ) +
c 1 DCVMEY( I )*DCVMEY( I ) )
WRITE (6,81) CEVARR( IEV )(1:LEVARR( IEV ) ),
1 DABSVX( I ), DABSVY( I ),
2 DCOVEL
ENDDO
81 FORMAT(A,' ',f12.6,1X,f12.6,1X,f12.6 )
c WRITE (6,81) CTAREB(1:LTAREB), CTARER(1:LTARER),
c 1 DXBMXA, DYBMYA, ND, DRESVN
c81 FORMAT(A,' minus ',A,1X,f10.4,1X,f10.4,1X,I4,1X,f20.4)
C
CALL EXIT(0)
99 CONTINUE
WRITE (6,'(A)') CMESS
CALL EXIT(1)
END
C