-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtot-pasq.f
14810 lines (14608 loc) · 487 KB
/
tot-pasq.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
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c hymain.for []
program hypoel
c hypoellipse -- john c. lahr
c non-unix version of main routine
c
c************************* notes for programmers **************************
c
c a binary search of the station list is used in the function phaind.
c if the search does not work on your computer, use the version
c of phaind that is commented out, which does a complete search.
c
c subroutines init, opfls, trvcon, trvdrv, and uamag must have
c $notruncate statements added for pc systems so that variable names
c longer than 6 characters may be used.
c
c all non-unix systems use hymain.for and init.for. the equivalent
c routines for unix systems are hypoe.c, setup_server.c,
c listen_serv.c, fdgetstr.c, cleanup.f, initial.f, and
c getbin.f.
c
c subroutines openfl and opfls, which open files, have differences between
c various systems.
c
c subroutine openfl also has differences between the unix-masscomp and
c the unix-sun systems.
c
c subroutines phagt and npunch use a back slash character, which
c must be doubled on unix systems.
c
c subroutines dubl, erset, jdate, and timit use non-standard fortran
c and must be modified for use with vax, pc, or unix computers.
c
c dubl - sets a double precision number equal to a
c single precision number
c erset - resets error limits on vax/vms computers
c jdate - gets current date and time from the operating system
c timit - times the execution on vax/vms computers
c
c alternat versions of the code is enclosed by 'c* unix', 'c* pc',
c or 'c* vax' comment statements in each of the above subroutines.
c
c get filenames, open files and write greeting
intype = 1
call init
c read station list, crustal model, and control records
call input1
c initialize summary of residuals
call lissum(1)
c read arrival times and locate earthquakes
call locate(-1)
stop
end
c end hymain
c adddly.for []
subroutine adddly(nmodel)
c read in delays for model number greater than 5
include 'params.inc'
parameter (ndly = 11)
character*4 iahead*60, msta*5, nsta*5, icard*110
integer punt
common /punt/ punt
common /char/ iahead, msta(npa), nsta(nsn), icard
common /dhip/ inpt,isa,ilis,inmain,injump
common /ilmpu/ ns
common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
common /int/ thk(lmax+1),lbeg(mmax),lend(mmax),vi(lmax),vsq(lmax),
* vsqd(lmmax,lmax),f(lmax,lmmax),kl,ivway,sthk(mmax),sthk1(mmax),
* sdly(ndly,nsn)
common /logfil/ logfil
character*4 stard, dnstrg, rshft
if(ns .eq. 0) then
write(punt, 18)
write(logfil, 62)
18 format('xxxerrorxxx delay models may not preceed the ',
* 'station list')
stop 'abort from adddly'
endif
20 read(inpt, '(a)', end=60) icard
stard = dnstrg(rshft(icard(1:4)))
if(stard .eq. ' end') return
c get station number (i) for station stard
c skip station if not on station list
c ho modificato qui: ho tolto phaind e messo ceck_staz
call ceck_staz(stard, nsta, ns, i, ierr_ceck)
if(ierr_ceck .ne. 0) then
write(punt, 34) stard, icard
write(logfil, 34) stard, icard
34 format(' ***>', a4, ' is not on station list, so these delays ',
* ' will not be used:', /, a)
goto 20
endif
c because standard fortran 77 to can not read an internal file with
c free format, file 14 must be used in the following code!
rewind 14
write(14, '(a)') icard(5:110)
rewind 14
c read(icard(5:110), *) dly(nmodel, i), sdly(nmodel, i)
read(14, *) dly(nmodel, i), sdly(nmodel, i)
goto 20
60 write(punt, 62)
write(logfil, 62)
62 format('xxxerrorxxx end of file found while reading additional',
* ' delays, so stop')
stop 'abort from adddly'
end
c end adddly
c adderr.for [unix]
subroutine adderr(zup, zdn)
c adds zup and zdn to the summary record and adds the event
c to file 4, and file 11
character*117 sumrec
character*117 arcrec
character*6 fmit
common /an/ n14, n15
c* (unix
call flush (14)
call flush (15)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
rewind 14
rewind 15
if(n14 .gt. 0) then
if(n14 .gt. 1) then
print *, ' n14 = ', n14
print *, ' logic error: n14 should be 0 or 1'
endif
do 20 i = 1, n14
read(14, '(a)') sumrec
call formal(zup, izup, 2, 0, fmit, azup)
if (fmit .eq. ' ') then
c write(sumrec(101:102), '(i2)') izup
write(sumrec(103:104), '(i2)') izup
else
c write(sumrec(101:102), fmit) azup
write(sumrec(103:104), fmit) azup
endif
call formal(zdn, izdn, 2, 0, fmit, azdn)
if (fmit .eq. ' ') then
c write(sumrec(103:104), '(i2)') izdn
write(sumrec(105:106), '(i2)') izdn
else
c write(sumrec(103:104), fmit) azdn
write(sumrec(105:106), fmit) azdn
endif
write(4, '(a)') sumrec(1:lentru(sumrec))
20 continue
endif
if(n15 .gt. 1) then
do 30 i = 1, n15
read(15, '(a)') arcrec
c if(arcrec(81:81) .ne. '/') then
if(arcrec(83:83) .ne. '/') then
write(11, '(a)') arcrec(1:lentru(arcrec))
else
c summary record
sumrec = arcrec
call formal(zup, izup, 2, 0, fmit, azup)
if (fmit .eq. ' ') then
c write(sumrec(101:102), '(i2)') izup
write(sumrec(103:105), '(i2)') izup
else
c write(sumrec(101:102), fmit) azup
write(sumrec(103:105), fmit) azup
endif
call formal(zdn, izdn, 2, 0, fmit, azdn)
if (fmit .eq. ' ') then
c write(sumrec(103:104), '(i2)') izdn
write(sumrec(105:106), '(i2)') izdn
else
c write(sumrec(103:104), fmit) azdn
write(sumrec(105:106), fmit) azdn
endif
c print *, 'in adderr, about to write summary rec to unit 11'
c print *, 'sumrec = ', sumrec
c print *, 'lentru = ', lentru(sumrec)
write(11, '(a)') sumrec(1:lentru(sumrec))
endif
30 continue
endif
c* (unix
call flush(4)
call flush(11)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
return
end
c end adderr
c azwtos.for []
subroutine azwtos
c azimuthal weighting of stations by quadrants
include 'params.inc'
character*4 iahead*60, msta*5, nsta*5, icard*110
common /char/ iahead, msta(npa), nsta(nsn), icard
common /amo/ tempg(npa), ntgap
common /gmost/ az(npa)
common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
common /qmost/ wt(npa),z
dimension tx(4),txn(4),ktx(4),kemp(npa),key(npa)
c count and sort by azimuth stations with weight .ne. 0
j=0
do 10 i=1,nr
if (wt(i) .eq. 0.) goto 10
j=j+1
tempg(j)=az(i)
10 continue
c divide into 4 zones with one axis bisecting largest gap between stations
call sort(tempg,key,j)
gap=tempg(1)+360.-tempg(j)
ig=1
do 20 i=2,j
dtemp=tempg(i)-tempg(i-1)
if (dtemp .le. gap) goto 20
gap=dtemp
ig=i
20 continue
tx(1)=tempg(ig)-0.5*gap
tx(2)=tx(1)+90.
tx(3)=tx(1)+180.
tx(4)=tx(1)+270.
do 30 i=1,4
txn(i)=0.
if (tx(i) .lt. 0.) tx(i)=tx(i)+360.
if (tx(i).gt.360.) tx(i)=tx(i)-360.
30 continue
call sort(tx,ktx,4)
do 80 i=1,nr
if (wt(i) .eq. 0.) goto 80
if (az(i) .gt. tx(1)) goto 50
40 txn(1)=txn(1)+1.
kemp(i)=1
goto 80
50 if (az(i) .gt. tx(2)) goto 60
txn(2)=txn(2)+1.
kemp(i)=2
goto 80
60 if (az(i) .gt. tx(3)) goto 70
txn(3)=txn(3)+1.
kemp(i)=3
goto 80
70 if (az(i) .gt. tx(4)) goto 40
txn(4)=txn(4)+1.
kemp(i) = 4
80 continue
xn=4
if (txn(1).eq.0.) xn=xn-1
if (txn(2).eq.0.) xn=xn-1
if (txn(3).eq.0.) xn=xn-1
if (txn(4).eq.0.) xn=xn-1
fj=j/xn
do 90 i=1,nr
if (wt(i) .eq. 0.) goto 90
ki=kemp(i)
c new weight equals old weight times total number of stations
c divided by number of zones contatning stations and divided
c------- again by number of zones containing this station
wt(i)=wt(i)*fj/txn(ki)
90 continue
return
end
c end azwtos
c back.for []
subroutine back (delat, delon, newlat, newlon, slat, slon)
c
c-------- back - calculate geocentric coordinates of secondary point from
c step in latitude (km) and longitude(km)
c
c input: delat change in earthquake latitude in km (northward positive)
c delon change in earthquake longitude in km (westward positive)
c output: newlat new earthquake geocentric latitude in radians
c newlon new earthquake longitude in radians
c
real newlat,newlon
parameter (pi = 3.14159265)
parameter (twopi = 2.0*pi)
parameter (halfpi = 0.5*pi)
parameter (rad = pi/180.)
parameter (deg = 1.0/rad)
parameter (equrad = 6378.2064)
parameter (polrad = 6356.5838)
parameter (flat = (equrad - polrad)/equrad)
c
st0 = cos(slat)
ct0 = sin(slat)
call cvrtop(delat,delon,delta,az1)
if (az1 .lt. 0.0) az1 = az1 + twopi
c use approximation of local radius for derivative of surface
c distance with geocentric latitude.
c more accurate formulation would be:
c drdth = -r**3 * cos(alat)*sin(alat)*( (1.-flat)**(-2) - 1. )/a**2
c dsd = sqrt(drdth**2 + r**2)
radius = (cos(slat)**2/equrad**2 + sin(slat)**2/polrad**2)**(-.5)
sdelt = sin(delta/radius)
cdelt = cos(delta/radius)
cz0 = cos(az1)
ct1 = st0*sdelt*cz0+ct0*cdelt
call cvrtop(st0*cdelt-ct0*sdelt*cz0, sdelt*sin(az1), st1, dlon)
newlat = atan2(ct1, st1)
newlon = slon + dlon
if (abs(newlon) .gt. pi) newlon = newlon - sign(twopi, newlon)
return
end
c end back
c block.for []
block data
c initialize constants in common statements
include 'params.inc'
parameter (ndly = 11)
logical good, eoff, supout
character*4 iahead*60, msta*5, nsta*5, icard*110, uacal*50
common /char/ iahead, msta(npa), nsta(nsn), icard
character*1 iclass
common /dbio/ iclass(0:4)
common /dbiln/ ioldq
common /dhio/ nedit,good
common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
common /dph/ noswt, eoff
common /dhin/ iglob, zup, zdn
common /dhip/ inpt,isa,ilis,inmain,injump
common /dhil/ iq,ilat,kms,icat
character*1 bksrc
common /dipu/ ipkdly, igsum, bksrc
common /dit/ tid(lmax,lmmax),did(lmax,lmmax),lowv,modin(mmax+3)
logical medmag
common /dinx/ imag,imagsv,medmag
common /dix/ iexcal, uacal
common /dmost/ ipun,ivlr,blank
common /imost/ test(100)
character*1 iqcls
common /il1/ iqcls
common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
common /logfil/ logfil
common /idno/ ksel,ksort
common /idt/ v(lmax2)
common /iox/ prr(nsn),iuses
common /it/ vpvs(lmax2),vpvsm(mmax+3),nttab
common /ix/ ir,qspa(9,40)
common /ohq/ gap, supout
common /pot/ ts(npa),ihr,model(npa), keyphi(npa)
common /reloc/ irelo, nreloc
common /rioq/ been,damp,dmpinc,igo
character*1 mgndx
common /xfmno1/ mgndx
common /ip1/ rsew(4)
common /ilmpu/ ns
data uacal/'pub1:[alaska.data]uofacal.dat'/, iexcal/0/
c data rsew/1., 3., 7.5, 15./, inpt/8/, inmain/8/, injump/12/
c revised default relative weights 4/22/89
data rsew/1., 5., 10., 20./, inpt/8/, inmain/8/, injump/12/
data logfil/6/,isa/0/,iqcls/'b'/,damp/0.0010/,bksrc/' '/
data iprun/1/, supout/.false./, ns/0/
data iclass/'0','a','b','c','d'/, nttab/0/
data iahead/' '/,mgndx/' '/,good/.true./
data blank/1.e20/
data test/100*1.23456/, v/lmax2*0.001/, model/npa*1/
data ivlr,kms,iprn,ipun,imag,iq,ksort,ksel,icat/0,1,1,0,0,2,0,1,1/
data ioldq/0/,ir/0/,lowv/1/,imagsv/0/,iuses/0/,ipkdly/1/
cds data tid,did/narray*0.0/, modin/13*0/, iglob/1/
data tid,did/narray*0.0/, modin/mmaxp3*0/, iglob/1/
data nedit/0/, eoff/.false./, igsum/1/, irelo/0/, nreloc/0/
data ilis/1/,medmag/.false./
end
c end block
c boxau.for []
subroutine boxau
c compute rms at auxiliary points and estimate quality
include 'params.inc'
parameter (ndly = 11)
real latep,lonep,latsv,lonsv
save latsv, lonsv, rmssv, zsv
save dez, t6, drsv, kdrsv, diag, kdiag, sum, aalz, aala
save aalo, idirec, idrc
character*1 iqa, ins, iew
character*4 iahead*60, msta*5, nsta*5, icard*110
integer punt
common /punt/ punt
common /char/ iahead, msta(npa), nsta(nsn), icard
common /bz/ na
common /bqz/ avrps,avuse
common /dmost/ ipun,ivlr,blank
character*1 iclass
common /dbio/ iclass(0:4)
common /gmost/ az(npa)
logical freor
common /hopq/ savla,savlo,savez,savor,freor
common /imost/ test(100)
common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
common /obcfn/ ain(npa)
character*1 kwr, iqdo
common /pbnoq/ kwr(npa),iqdo(npa)
common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
common /phoqn/ inst,knst
common /qmost/ wt(npa),z
common /qmost1/ lonep,ni,latep
common /qgnotx/ delta(npa)
common /qgo/ org
common /rob/ v(4,4),noaxp
common /rbno/ idip(3),iaaz(3),a(3,3)
common /rfgnoq/ se(4)
common /tmost/ x(4,npa)
common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
common /zqr/ fno
dimension drsv(14),kdrsv(14),diag(7),kdiag(7)
dimension sum(5)
dimension aalz(14),aala(14),aalo(14)
dimension vec(3)
character*2 idirec(7), idrc(7)
data idirec/' n','se','sw','nw','ne',' z',' e'/
data aala/1.732,1.,1.,1.,1.,0.,0.,0.,0.,-1.,-1.,-1.,-1.,-1.732/,
* aalo/0.,1.,-1.,1.,-1.,0.,1.732,-1.732,0.,1.,-1.,1.,-1.,0./,
* aalz/0.,-1.,-1.,1.,1.,-1.732,0.,0.,1.732,-1.,-1.,1.,1.,0./
kgoon = 0
c
c write heading
17 if( (iprn .gt. 2) .or. (noaxp .eq. 1) ) write(punt,18)
18 format(/50h lat lon z avrps no ,
* 7x,3hrms,25x,4hdrms/)
c
c save current values
rmssv = rms
latsv = latep
lonsv = lonep
zsv = z
orgsv = org
instsv = inst
c set inst = -9 so quakes will compute fixed location rms but will
c not call regres.
inst = -9
freor = .true.
t6 = abs(test(6))/1.732
dez = t6
do 70 na = 1, 14
call back(t6*aala(na),t6*aalo(na),savla,savlo,latsv,lonsv)
savez = zsv +aalz(na)*dez
if(z .lt. 0.) z = 0.0
call quakes
c output rms error af auxiliary points
rmsx = rms
drms = rmsx - rmssv
if(iprn .eq. 4) then
c calcute predicted rms at aux points
tot = 0.0
vec(1) = aalo(na)/1.732
vec(2) = aala(na)/1.732
vec(3) = aalz(na)/1.732
do 48 i = 1,3
do 46 j = 1,3
tot = tot + vec(i)*vec(j)*v(i,j)
46 continue
48 continue
pdrms = sqrt(rmssv**2 + tot*test(6)*test(6)/fno)
drms = rmsx - pdrms
endif
call unfold2(savla,savlo,la,ins,ala,lo,iew,alo)
drsv(na) = drms
if( (iprn .ge. 2) .and. (noaxp .ne. 1) ) then
goto (1,2,3,4,5,6,7,8,9,10,11,4,13,1), na
2 write(punt,801) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
801 format(4f10.2,i10,f10.2,10x,f5.2)
goto 50
3 write(punt,802) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
802 format(4f10.2,i10,f10.2,24x,1h-,8x,f5.2/96x,1h.)
goto 50
1 write(punt,803) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
803 format(/4f10.2,i10,f10.2,23x,f5.2/)
goto 50
4 write(punt,804) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
804 format(4f10.2,i10,f10.2,13x,f5.2,19x,1h.)
goto 50
5 write(punt,805) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
805 format(4f10.2,i10,f10.2,12x,'|',23x,f5.2)
goto 50
6 write(punt,806) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
806 format(4f10.2,i10,f10.2,20x,f5.2,10x,'|')
goto 50
7 write(punt,807) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
807 format(4f10.2,i10,f10.2,3x,f5.2/)
write(punt,815) rmssv
815 format(50x,f10.2,12x,'|',10x,5h 0.00,10x,'|'/95x,'|')
goto 50
8 write(punt,808) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
808 format(4f10.2,i10,f10.2,43x,f5.2)
goto 50
9 write(punt,809) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
809 format(4f10.2,i10,f10.2,26x,f5.2)
goto 50
10 write(punt,810) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
810 format(4f10.2,i10,f10.2,10x,f5.2,23x,'|')
goto 50
11 write(punt,811) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
811 format(4f10.2,i10,f10.2,13x,1h.,10x,1h-,8x,f5.2/
* 74x,1h.,21x,1h.)
goto 50
13 write(punt,813) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
813 format(4f10.2,i10,f10.2,27x,1h-,8x,f5.2)
endif
50 continue
70 continue
c ouality output
i = 0
ii = 15
do 80 iii = 1,7
i = i + 1
ii = ii - 1
80 diag(iii) = (drsv(i) + drsv(ii))/2.0
call sort(diag, kdiag, 7)
call sort(drsv, kdrsv, 14)
do 90 i = 1,7
j = kdiag(i)
idrc(i) = idirec(j)
90 continue
if((noaxp .eq. 0) .and. (iprn .ge. 0)) write(punt,100) idrc, diag
100 format(/25x,18hquality evaluation//,
* 34h diagonals in order of strength ,7(4x,a2)/,
* 19h ave. of end points,15x,7f6.2//)
sum(1) = 0.0
do 110 i=1,14
c do not use upper point in average if limited by surface
if((z.lt.abs(test(6))*1.732).and.(kdrsv(i).eq. 6)) goto 110
sum(1) = sum(1) + drsv(i)
110 continue
avdrms = sum(1)/14.
if(z .lt. abs(test(6))*1.732) avdrms = sum(1)/13.
drmin = drsv(1)
icmin = kdrsv(1)
if((z .ge. abs(test(6))*1.732).or.(kdrsv(1) .ne. 6)) goto 115
drmin = drsv(2)
icmin = kdrsv(2)
115 jab = 4
if((nrwt.ge. 4.).and.(rmssv.le. 0.4).and.(avdrms.ge. 0.5)) jab = 3
if((nrwt.ge.5.).and.(rmssv.le. 0.4).and.(drmin .ge. 0.15)) jab=2
if((nrwt.ge.6.).and.(rmssv.le. 0.2).and.(drmin .ge. 0.30)) jab=1
iqa = iclass(jab)
if((noaxp.eq.0) .and. (iprn .ge. 0))
* write(punt,120)nrwt,rmssv,drmin,avdrms,iqa
120 format(10x,50h number rms min drms ave drms quality/
* 10x,i10,3f10.2,9x,a1/)
c
c restore variable values
latep = latsv
lonep = lonsv
z = zsv
org = orgsv
inst = instsv
if (drmin .gt. -.03 .or. test(6) .gt. 0.) then
c restore all values, such as azimuth, angle of incidence, and
c standard error to those of final solution for npunch.
130 iprnsv = iprn
iprn = -2
inst = 9
freor = .false.
savla = latsv
savlo = lonsv
savez = zsv
savor = orgsv
call quakes
c call output(0)
c call output(1)
iprn = iprnsv
inst = instsv
return
else
kgoon = kgoon + 1
c the first time an event reaches here, kgoon will equal 1.
c the 2nd time kgoon will equal 2, so the event will not be rurun again.
if((inst .ne. 9) .and. (inst .ne. 8)) return
if(kgoon .ge. 2) then
c do the same as above, starting with statement 130
c goto 130
iprnsv = iprn
iprn = -2
inst = 9
freor = .false.
savla = latsv
savlo = lonsv
savez = zsv
savor = orgsv
call quakes
c call output(0)
c call output(1)
iprn = iprnsv
inst = instsv
return
endif
call back(t6*aala(icmin),t6*aalo(icmin),savla,savlo,latsv,lonsv)
savez = zsv + aalz(icmin)*dez
if(z .lt. 0.0) z = 0.0
write(punt,1000)
1000 format(///' *** run again starting at best nearby point ***',/)
call quakes
goto 17
endif
end
c end boxau
c cosh.for []
function cosh(x)
cosh = (exp(x) + exp(-x))/2.
return
end
c end cosh
c coshi.for []
function coshi(x)
p = x + sqrt(x*x - 1.)
coshi = alog(p)
return
end
c end coshi
c critic.for []
subroutine critic(delta,az,w,iuse,nrp,nr,ldx)
c select critical stations
c 1) closest 4 with p-phase readings.
c 2) additional with p- or s-phase readings if they reduce the a gap > 72 deg
c by 5 deg or more.
c 3) s is used at critical stations. if there are no s phases selected,
c then s is used from the closest non-critical station.
include 'params.inc'
logical swtchp, swtchs
dimension dtemp(npa),key(npa)
dimension delta(npa),az(npa),w(npa),iuse(npa),ldx(npa)
integer punt
common /punt/ punt
write(punt,20)
20 format(' determining which stations are critical.')
do 40 j = 1, nr
iuse(j) = 0
if (w(j) .eq. 0.) iuse(j) = 1
40 continue
do 50 j = 1, nrp
dtemp(j) = delta(j)
50 continue
call sort(dtemp,key,nrp)
ns = 0
nwt = 0
c always use first 4 with p-wt .gt. 0.0
do 100 i = 1, nrp
j = key(i)
c check for p
if (w(j) .gt. 0.0) then
nwt = nwt + 1
iuse(j) = 1
c add s only if p is selected
if (ldx(j) .ne. 0) then
k = ldx(j)
if (w(k) .gt. 0.0) then
iuse(k) = 1
ns = ns + 1
endif
endif
endif
if (nwt .ge. 5) goto 125
100 continue
c fall through if there are 4 or fewer p phases that have wt <> 0
125 if (nwt .eq. 0) then
sge72 = 400.
nge72 = 0.
else
call sumgap(iuse, az, nr, sge72, nge72, w)
endif
c try one at a time for stations that reduce gap.
do 200 i = 1, nrp
j = key(i)
k = ldx(j)
if (iuse(j) .eq. 1) then
c p already in use
if (k .gt. 0) then
if (iuse(k) .eq. 1) then
c s also in use, so skip on to next phase
goto 200
endif
else
c no s
goto 200
endif
endif
c either p and/or s are not now being used - check weights
if (w(j) .eq. 0.) then
if (k .gt. 0) then
c there is an s phase
if (w(k) .eq. 0.) then
c p and s weights are 0, so skip to next phase
goto 200
endif
endif
endif
c either p and/or s has non zero weight and is not yet being used
c try adding the p or s to test gap reduction
swtchp = .false.
swtchs = .false.
if ((iuse(j) .eq. 0) .and. (w(j) .gt. 0.)) then
swtchp = .true.
iuse(j) = 1
else
swtchs = .true.
iuse(k) = 1
endif
call sumgap(iuse, az, nr, sgtr, ngtr, w)
c check if the number of gaps larger than 72 has been increased
if (ngtr .gt. nge72) then
nge72 = ngtr
c check if the reduction in sum of gaps > 72 has been reduced by 5 or more deg
else if (sgtr .le. (sge72 - 5.)) then
nge72 = ngtr
sge72 = sgtr
else
c there was no or not enough improvement
if(swtchp) iuse(j) = 0
if(swtchs) iuse(k) = 0
goto 200
endif
c this phase reduced the gap, so use it!
if (swtchp) then
nwt = nwt + 1
c add s if available
if (k .ne. 0) then
if (w(k) .gt. 0.0) then
iuse(k) = 1
ns = ns + 1
endif
endif
endif
if (swtchs) then
ns = ns + 1
endif
200 continue
400 if (ns .gt. 1) goto 500
c use closest s if none found at any critical station
do 350 i = 1,nrp
j = key(i)
if (ldx(j) .eq. 0) goto 350
k = ldx(j)
if (w(k) .le. 0.) goto 350
iuse(k) = 1
iuse(i) = 1
goto 500
350 continue
500 return
end
c end critic
c cvrtop.for []
subroutine cvrtop(x, y, r, theta)
c
c-------- bruce julian
c
c-------- cvrtop - convert from rectangular to polar coordinates
c
c (output - may overlay x, y)
c
c-------- standard fortran funct. required: atan2
c-------- funct. required: hypot
c
r = hypot(x, y)
theta = 0.
if ((y .ne. 0.) .or. (x .ne. 0.)) theta = atan2(y, x)
return
end
c
c end cvrtop
c cyldly.for []
subroutine cyldly
* (kno, alat, alon, lat1, lon1, x0, y0, z, dly, sdly, iprn,
* modset, test8)
save cyldy, cylrd, cylrd1, cylup, cylup1, cyldn, cyldn1,
* xc, yc, ncyl, setcyl
real lat1, lon1
character*1 ins, iew
include 'params.inc'
parameter (ndly = 11)
parameter (pi = 3.14159265)
parameter (rad = pi/180.)
parameter (mxcyl = 50)
logical setcyl
real cyld
c cyld horizontal distance to inner cylinder
real cyldin(mxcyl)
c cyldin(i) distance to cylinder center for eq within
c inner radius
integer cyldy(mxcyl)
c cyldy(i) delay assigned to cylinder i
integer cylmd(mxcyl)
c cylmd(i) velocity model assigned to cylinder i
real cylrd(mxcyl)
c cylrd(i) inner radius of cylinder number i
real cylrd1(mxcyl)
c cylrd1(i) outer radius of cylinder number i
integer cyldm(mxcyl)
c cyldm(i) delay model for transition table entry i
integer cyldmin(mxcyl)
c cyldmin(i) delay model for inner table entry i
real cylup(mxcyl)
c cylup(i) upper limit of inner cylinder i
real cylup1(mxcyl)
c cylup1(i) upper limit of outer cylinder i
real cyldn(mxcyl)
c cyldn(i) lower limit of inner cylinder i
real cyldn1(mxcyl)
c cyldn1(i) lower limit of outer cylinder i
real cylwt(mxcyl)
c cylwt(i) weight derived for transition zone i
real dly(ndly,nsn)
c dly(i, j) p-delay for model i, station j
integer ntrans
c ntrans number of entries in transition table
real sdly(ndly,nsn)
c sdly(i, j) s-delay for model i, station j
real xc(mxcyl)
c xc(i) x coordinate of center of cylinder i
real yc(mxcyl)
c yc(i) y coordinate of center of cylinder i
real z
c z current trial eq depth
real zmarg
c zmarg vertical distance from event to endcap
c (negative for lower endcap)
data setcyl /.false./
modset = 0
if(iprn .ge. 5) then
print *, ' '
print *, 'begin sub. cyldly to select delay and velocity'
print *, 'model(s) based on cylindrical domains.'
endif
if (.not. setcyl) then
cd print *, 'read in cylinder definitions with cylget'
call cylget
* ( mxcyl, cyldy, cylmd, cylrd, cylrd1, cylup, cylup1, cyldn,
* cyldn1, xc, yc, ncyl, lat1, lon1, test8)
setcyl = .true.
endif
c compute the x,y coordinates of the current trial epicenter
call delaz(lat1, lon1, delt, deldeg, azz, alat, alon)
x0 = delt*sin(azz*rad)
y0 = delt*cos(azz*rad)
if(iprn .ge. 5) then
call unfold2(alat, alon, la, ins, ala, lo, iew, alo)
print *, 'current epicenter = ', la, ins, ala, lo, iew, alo
print *, 'location of epicenter wrt reference station is'
print *, 'azimuth (deg) = ', azz
print *, 'distance (km) = ', delt
print *, 'x, y = ', x0, y0
print *, 'loop through ', ncyl, ' regions. regions must be in '
print *, 'order of preference in cases of overlap.'
endif
ntrans = 0
ninner = 0
do 20 i = 1, ncyl
cyld = sqrt((x0 - xc(i))**2 + (y0 - yc(i))**2)
if(iprn .ge. 5) then
print *, 'x, y of cylinder ', i, ' is ', xc(i), yc(i)
print *, 'dist to this cylinder ',
* 'which uses delay model ', cyldy(i), ' is ', cyld
print *, 'z = ', z, ' cylup,dn(i) = ', cylup(i), cyldn(i)
endif
if ((cyld .le. cylrd(i)) .and. (z .ge. cylup(i)) .and.
* (z .le. cyldn(i))) then
if(iprn .ge. 5)
* print *, 'location is within inner cylinder'
ninner = ninner + 1
cyldin(ninner) = cyld
cyldmin(ninner) = i
else if ((cyld .le. cylrd1(i)) .and. (z .ge. cylup1(1)) .and.
* (z .le. cyldn1(i))) then
c
cd print *, 'location is within transition zone, so compute '
cd print *, 'weight and add to list'
ntrans = ntrans + 1
cyldm(ntrans) = cyldy(i)
cd print *, 'zmarg is distance into upper or lower cap'
zmarg = 0.0
if ((z .lt. cylup(i)) .and. (z .gt. cylup1(i)))
* zmarg = cylup(i) - z
if ((z .gt. cyldn(i)) .and. (z .lt. cyldn1(i)))
* zmarg = cyldn(i) - z
cd print *, 'note that zmarg will be negative for lower ',
cd * 'cap region'
if (zmarg .eq. 0.0) then
c
cd print *, 'adjacent to cylinder '
cd print *, 'cyld - cylrd(i) ', cyld - cylrd(i)
cd print *, 'pi ', pi
cd print *, 'cylrd1(i) - cylrd(i) ', cylrd1(i) - cylrd(i)
cylwt(ntrans) = 0.5 + 0.5*cos( pi*(cyld - cylrd(i)) /
* (cylrd1(i) - cylrd(i)) )
else if (cyld .le. cylrd(i)) then
c
cd print *, 'within end caps'
if (zmarg .gt. 0) then
cd print *, 'within upper cap'
cylwt(ntrans) = 0.5 + 0.5*cos( pi*zmarg /
* (cylup(i) - cylup1(i)) )
else
c
cd print *, 'within lower cap'
cylwt(ntrans) = 0.5 + 0.5*cos( pi*zmarg /
* (cyldn(i) - cyldn1(i)) )
endif
else
c
cd print *, 'within corner zone'
if (zmarg .gt. 0) then
cd print *, 'within upper corner zone'
cylwt(ntrans) = 0.5 + 0.5*cos(pi*sqrt
* ( ((cyld - cylrd(i)) /
* (cylrd1(i) - cylrd(i)))**2 +
* (zmarg /
* (cylup(i) - cylup1(i)))**2 ) )
else
cd print *, 'within lower corner zone'
cylwt(ntrans) = 0.5 + 0.5*cos(pi*sqrt
* ( ((cyld - cylrd(i)) /
* (cylrd1(i) - cylrd(i)))**2 +
* (zmarg /
* (cyldn(i) - cyldn1(i)))**2 ) )
endif
endif
endif
20 continue
c use the parameters for the cylinder for which the eq is
c closest to the center
if (ninner .ne. 0) then
if(ninner .eq. 1) then
kno = cyldy(cyldmin(ninner))
modset = cylmd(cyldmin(ninner))
return
else
ntouse = 1
do 22 i = 2, ninner
if(cyldin(i) .lt. cyldin(ntouse)) then
ntouse = i
endif
22 continue
kno = cyldy(cyldmin(ntouse))
modset = cylmd(cyldmin(ntouse))
return
endif
endif
if (ntrans .eq. 0) then
cd print *, 'not within any cylinders, so use default model (1)'
kno = 1
return
else
cd print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd print *, 'models: ', (cyldm(i), i = 1, ntrans)
c reduce list if there are 2 or more entries
if (ntrans .gt. 1) call cylred(cyldm, cylwt, ntrans, sumwt)
kno = ndly
c
if (ntrans .eq. 1) then
if (cylwt(1) .ge. 1.0) then
kno = cyldm(1)
return
else
c fill in with default model
ntrans = 2
cylwt(2) = 1.0 - cylwt(1)
cyldm(2) = 1
if (iprn .ge. 5) then
print *, 'weights: ', (cylwt(i), i = 1, ntrans)
print *, 'models: ', (cyldm(i), i = 1, ntrans)
endif
endif
c
else if (ntrans .eq. 3) then
cylwt(1) = cylwt(1)/sumwt
cylwt(2) = cylwt(2)/sumwt
cylwt(3) = cylwt(3)/sumwt
sumwt = 1.0
cd print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd print *, 'models: ', (cyldm(i), i = 1, ntrans)
else
c
c ntrans = 2
if (sumwt .gt. 1.0) then
cylwt(1) = cylwt(1)/sumwt
cylwt(2) = cylwt(2)/sumwt
else
c fill in with default model
ntrans = 3
cylwt(3) = 1. - sumwt
sumwt = 1.0
cyldm(3) = 1
cd print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd print *, 'models: ', (cyldm(i), i = 1, ntrans)
endif
endif
endif
c compute combined delays
do 30 i = 1, nsn
dly(ndly,i) = 0.0
sdly(ndly,i) = 0.0