Changeset 3168
 Timestamp:
 12/05/07 11:13:22 (12 years ago)
 Location:
 trunk/src/seismic_processing/rayloc_ew
 Files:

 6 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/seismic_processing/rayloc_ew/PROGRAMMING_NOTES
r1669 r3168 1 1 2 Files rayloc_message_rw.h and rayloc_message_rw.c contain 2 3 API for reading and writing TYPEi_RAYLOC messages. … … 14 15 Ilya Dricker (i.dricker@isti.com) 15 16 and ISTI team 17 18 20071205 update from Mitch Withers. Synced rayloc_ew with ray_loc 19 from Hydra. Involved changes to ellip.f and robust_util.f. Also 20 needs a new prm file that I had previously allowed to reside in 21 the glass model dir. Added that dir and its contents along with 22 the new res_ac.prm file to rayloc_ew's directory. There are also 23 problems with the hydra ew that should be fixed. 24 25 A known feature is rayloc will complain with a less than illuminating 26 error when a station is missing from the station file. Something like: 27 28 20071205_UTC_15:12:50 lib_rayloc: Failed to write Input Structure to the file /g 29 aia/home/mwithers/Projects/EWSupport/Working/run_smeagol/model/RayLocInputa.txt 30 20071205_UTC_15:12:50 rayloc_ew: lib_rayloc FAILED: retVal = 1 31 20071205_UTC_15:12:50 rayloc_ew: Output Message size = NULL 32 20071205_UTC_15:12:50 ERROR: rayloc_MessageToRaylocHeader: RAYLOC HEADER MESSAGE 33 is NULL 34 20071205_UTC_15:12:50 Failed to convert RAYLOC message into structure 35 36 You have been warned. 
trunk/src/seismic_processing/rayloc_ew/ellip.f
r1669 r3168 1 1 real function ellip(phase,edist,edepth,slat,bazim) 2 2 3 c IGD 06/01/04 Added dirname common for portability 4 3 c MMW 20071129 merged Hydra and ew version 4 5 c IGD 06/01/04 Added dirname common for portability 6 5 7 C========================================================================== 6 8 C … … 43 45 c Changed lue to merr and allocated lut with freeunit. 44 46 C========================================================================= 45 common /my_dir/dirname 46 character*1001 dirname 47 47 48 common /my_dir/dirname 49 character*1001 dirname 50 48 51 character *(*) phase 49 52 50 53 real edist,edepth,ecolat,bazim,slat,azim,degrad, 51 54 ^ sc0,sc1,sc2,tcor, … … 65 68 66 69 c.. 67 68 70 parameter(NUMPH=57) 69 71 … … 112 114 113 115 save INIT,ECHO1,ECHO2,/tau1/,/tau2/,/tau3/,Nd,deldst,degrad 114 115 116 116 117 c.. fixed … … 167 168 168 169 c Make sure we use dirname if one is provided 169 if(dirname(1:2).eq.'. ') then 170 open(lut,file='tau.table',form='formatted', 171 > iostat=ios,err=998,status='old') 172 else 173 open(lut, file=dirname(1:ntrbl(dirname))//'/'//'tau.table',form='formatted', 174 > iostat=ios,err=998,status='old') 175 endif 176 rewind(lut) 170 if(dirname(1:2).eq.'. ') then 171 open(lut,file='tau.table',form='formatted', 172 > iostat=ios,err=998,status='old') 173 else 174 open(lut, file=dirname(1:ntrbl(dirname))//'/'//'tau.table',form='formatted', 175 > iostat=ios,err=998,status='old') 176 endif 177 178 c open(lut,file='tau.table',form='formatted', 179 c > iostat=ios,err=998,status='old',readonly,shared) 180 rewind(lut) 181 177 182 ip=0 178 183 20 continue … … 207 212 208 213 21 continue 209 210 c IGD 06/18/04 Closed statemet included 211 close(unit=lut) 214 212 215 if(ECHO1) then 213 216 write(merr,*) 'Number of phases: ',ip1 … … 262 265 call phase_alias(phase,edist,ip) 263 266 endif 264 267 265 268 if(ip.gt.0) then 266 269 if(ECHO1) … … 275 278 return 276 279 endif 277 280 278 281 c... distance index 279 282 idist = 1 + int( (edistdi1(ip))/ deldst ) 280 283 if(edist.lt.di1(ip)) idist =1 281 284 if(edist.gt.di2(ip)) idist= np(ip)1 282 285 283 286 c... depth index 284 c IGD 07/29/04 Next line to prevent crashes if jdepth is undefined285 c IGD 07/29/04 It is not clear to me how much we screw up processing286 jdepth = Nd287 287 do 25 j = 1,Nd1 288 288 if ((dpth(j).le.edepth).and.(dpth(j+1).ge.edepth))then … … 291 291 endif 292 292 25 continue 293 c Extrapolate below 700 km (& make behavior deterministic). Added by 294 c R. Buland on 18 March 2004. 295 jdepth=Nd1 293 296 26 continue 294 297 295 298 if(ECHO2) write(merr,*) 'idist, jdepth;',idist,jdepth 296 299 297 300 c... Compute tauvalues and ellip 298 301 c need to allow for zero entries (where phase 299 302 c description strongly depth dependent) 300 303 301 304 c tau0 302 305 a0 = t0(ip,idist,jdepth) … … 356 359 357 360 ellip=tcor 361 358 362 return 359 363 … … 403 407 character*(*) phase 404 408 405 if(phase(1:3).eq.'Pn ') then 409 if(phase(1:3).eq.'Pg '.or.phase(1:3).eq.'Pb '.or. 410 1 phase(1:3).eq.'Pn ') then 406 411 c phase='P ' 407 412 ip=2 408 else if(phase(1:3).eq.'Sn ') then 413 else if(phase(1:3).eq.'Sg '.or.phase(1:3).eq.'Sb '.or. 414 1 phase(1:3).eq.'Sn ') then 409 415 c phase='S ' 410 416 ip=33 411 else if(phase(1:4).eq.'pPn ') then 417 else if(phase(1:4).eq.'pPg '.or.phase(1:4).eq.'pPb '.or. 418 1 phase(1:4).eq.'pPn ') then 412 419 c phase='pP ' 413 420 ip=8 … … 415 422 c phase='pP ' 416 423 ip=8 417 else if(phase(1:5).eq.'pwPn ') then 424 else if(phase(1:5).eq.'pwPg '.or.phase(1:5).eq.'pwPb '.or. 425 1 phase(1:5).eq.'pwPn ') then 418 426 c phase='pP ' 419 427 ip=8 420 else if(phase(1:4).eq.'sPn ') then 428 else if(phase(1:4).eq.'sPg '.or.phase(1:4).eq.'sPb '.or. 429 1 phase(1:4).eq.'sPn ') then 421 430 c phase='sP ' 422 431 ip=13 423 else if(phase(1:4).eq.'pSn ') then 432 else if(phase(1:4).eq.'pSg '.or.phase(1:4).eq.'pSb '.or. 433 1 phase(1:4).eq.'pSn ') then 424 434 c phase='pS ' 425 435 ip=37 426 else if(phase(1:4).eq.'sSn ') then 436 else if(phase(1:4).eq.'sSg '.or.phase(1:4).eq.'sSb '.or. 437 1 phase(1:4).eq.'sSn ') then 427 438 c phase='sS ' 428 439 ip=40 429 else if(phase(1:4).eq.'SPn ') then 440 else if(phase(1:4).eq.'SPg '.or.phase(1:4).eq.'SPb '.or. 441 1 phase(1:4).eq.'SPn ') then 430 442 c phase='SP ' 431 443 ip=55 432 else if(phase(1:4).eq.'SnP ') then 444 else if(phase(1:4).eq.'SgP '.or.phase(1:4).eq.'SbP '.or. 445 1 phase(1:4).eq.'SnP ') then 433 446 c phase='SP ' 434 447 ip=55 435 else if(phase(1:4).eq.'PSn ') then 448 else if(phase(1:4).eq.'PSg '.or.phase(1:4).eq.'PSb '.or. 449 1 phase(1:4).eq.'PSn ') then 436 450 c phase='PS ' 437 451 ip=56 438 else if(phase(1:5).eq.'PnPn ') then 452 else if(phase(1:5).eq.'PgPg '.or.phase(1:5).eq.'PbPb '.or. 453 1 phase(1:5).eq.'PnPn ') then 439 454 c phase='PP ' 440 455 ip=30 441 else if(phase(1:5).eq.'SnSn ') then 456 else if(phase(1:5).eq.'SgSg '.or.phase(1:5).eq.'SbSb '.or. 457 1 phase(1:5).eq.'SnSn ') then 442 458 c phase='SS ' 443 459 ip=53 … … 463 479 c phase="S'S' ' 464 480 ip=54 465 else if(delta.le.100.0.and.phase.eq.'pPdiff ') then 481 else if(phase.eq.'Pdif') then 482 c phase='Pdiff' 483 ip=3 484 else if(phase.eq.'Sdif') then 485 c phase='Sdiff' 486 ip=34 487 else if(delta.le.100.0.and.phase.eq.'pPdif ') then 466 488 c phase='pP ' 467 489 ip=8 468 else if(delta.le.100.0.and.phase.eq.'pwPdif f') then490 else if(delta.le.100.0.and.phase.eq.'pwPdif ') then 469 491 c phase='pP ' 470 492 ip=8 471 else if(delta.le.100.0.and.phase.eq.'sPdif f') then493 else if(delta.le.100.0.and.phase.eq.'sPdif ') then 472 494 c phase='sP ' 473 495 ip=13 474 else if(delta.le.100.0.and.phase.eq.'pSdif f') then496 else if(delta.le.100.0.and.phase.eq.'pSdif ') then 475 497 c phase='pS ' 476 498 ip=37 477 else if(delta.le.100.0.and.phase.eq.'sSdif f') then499 else if(delta.le.100.0.and.phase.eq.'sSdif ') then 478 500 c phase='sS ' 479 501 ip=40 480 else if(delta.le.165.0.and.phase.eq.'PKPdif f') then502 else if(delta.le.165.0.and.phase.eq.'PKPdif ') then 481 503 c phase='PKPbc ' 482 504 ip=5 483 else if(delta.le.165.0.and.phase.eq.'pPKPdif f') then505 else if(delta.le.165.0.and.phase.eq.'pPKPdif ') then 484 506 c phase='pPKPbc ' 485 507 ip=10 486 else if(delta.le.165.0.and.phase.eq.'sPKPdif f') then508 else if(delta.le.165.0.and.phase.eq.'sPKPdif ') then 487 509 c phase='sPKPbc 488 510 ip=15 
trunk/src/seismic_processing/rayloc_ew/input.f
r1669 r3168 79 79 c 80 80 c Reset to generic phase codes so they will be used and reidentified. 81 82 83 84 81 c if(phase.eq.'Pg'.or.phase.eq.'Pn'.or.phase.eq.'Pb'.or. 82 c 1 phase.eq.'P*') phase='P' 83 c if(phase.eq.'Sg'.or.phase.eq.'Sn'.or.phase.eq.'Sb'.or. 84 c 1 phase.eq.'S*') phase='S' 85 85 c 86 86 c Check to see if you should use the phase in the earthquake location 
trunk/src/seismic_processing/rayloc_ew/libtau.f
r1669 r3168 887 887 if(ln.le.0) ln=index(phcd(jb),' ')1 888 888 if(ln.le.0) ln=len(phcd(jb)) 889 phnm(n)=phcd(jb)(1:ln)//'dif f'889 phnm(n)=phcd(jb)(1:ln)//'dif' 890 890 10 continue 891 891 return 
trunk/src/seismic_processing/rayloc_ew/makefile.sol
r2055 r3168 7 7 # Revision history: 8 8 # $Log$ 9 # Revision 1.3 2007/12/05 16:13:22 withers 10 # Synced with hydra version of rayloc 11 # 9 12 # Revision 1.2 2006/01/19 19:04:55 friberg 10 13 # rayloc upgraded to use a library version of rayloc_message_rw.c … … 41 44 42 45 CFLAGS = ${GLOBALFLAGS} g I ${I} DUSE_LOGIT 43 FFLAGS = g 46 FFLAGS = g f77=tab 44 47 B = $(EW_HOME)/$(EW_VERSION)/bin 45 48 L = $(EW_HOME)/$(EW_VERSION)/lib … … 81 84 82 85 rayloc_ew: $(RAYLOC_EW_OBJ) $(F_OBJS) $(C_OBJS) $(F_STUB) 83 ${F77} o rayloc_ew $(RAYLOC_EW_OBJ) $(F_OBJS) $(C_OBJS) $(F_STUB) lm lposix4 86 ${F77} o rayloc_ew $(RAYLOC_EW_OBJ) $(F_OBJS) $(C_OBJS) $(F_STUB) lm lposix4 84 87 85 88 $(TARGET1): $(F_OBJS) $(F_MAIN) 
trunk/src/seismic_processing/rayloc_ew/robust_util.f
r1669 r3168 206 206 c 207 207 include 'locprm.inc' 208 character*2 1flgs,cstor208 character*22 flgs,cstor 209 209 c Logical unit numbers. 210 210 common /unit/ min,mout,merr … … 806 806 c 807 807 include 'locprm.inc' 808 character*8 phcd1(maxtt),ophase, ph 808 logical phasid,log 809 809 c integer stdin, stdout 810 real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),811 1 ellip, elcor812 810 double precision a,b,c,x,s,avdr 813 811 c Physical constants. … … 821 819 c Stationreading parameters. 822 820 include 'pickcm.inc' 823 data dtds/19.96/824 821 c 825 822 c stdin = 5 … … 849 846 1 ierr=normeq(navsum,ihold) 850 847 c 851 c Initialize the depth852 c853 2 if (odep .le. 0.0) odep = 1.0854 zs = odep855 call depset(zs,usrc)856 c857 848 c Loop over all stationreadings. 858 849 c 859 850 2 do 3 j = 1, nphase 860 851 c 861 852 c Fetch data from station parameter storage. … … 872 863 1 delta,azim) 873 864 c 874 n = 0 875 call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1) 876 c 877 c Find the computed traveltimes and partials for phases other than P and S 878 c 879 do 111 i = 1, n 880 ph = phcd1(i) 881 ilen = ntrbl(phase) 882 if(phase(1:ilen) .eq. ph(1:ilen)) then 883 travt = tt1(i) 884 dtdd = dtdd1(i) 885 dtdh = dtdh1(i) 886 ophase = phase 887 phase = ph 888 go to 100 889 endif 890 111 continue 891 c 892 c Phase not found for that distance and depth 893 c 894 pflg = ' ' 895 res=0. 896 go to 5 897 c 898 100 continue 899 c 900 c This was the original call to get the ellipticity and elevation correction, 901 c which is replaced by subroutine used by Engdahl. Correction for the 902 c elevation was taken out of subroutine datum 903 c call datum(elev,dtdd,olat,delta,azim*deg,tcor,slatc) 904 c 905 azimuth = azim*deg 906 tcor=ellip(phase,delta,zs,glat,azimuth) 907 if(dtdd.gt.dtds) dtdd=dtds 908 elcor=(elev/5.57)*sqrt (1.0(5.57*dtdd/(radius*rad))**2) 909 res=trvtimosectravt(tcor+elcor) 910 c 911 c 912 c write(merr,333)j,sta,ophase,ph,pflg,delta,azimuth, 913 c 2 trvtim,res,dtdd,dtdh 914 c333 format('nclose: ',i5,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x, 915 c 1 f8.2,1x,f6.2,1x,f8.4,1x,f8.4) 865 log=phasid() 916 866 c 917 867 c All done. Push data back out to station parameter storage. 918 868 c 919 5call drdwr(2,j,1,13)869 call drdwr(2,j,1,13) 920 870 3 continue 921 871 c … … 973 923 c by E. R. Engdahl. 974 924 c 975 double precision a,b,c,x,s,avdr,eps 976 c double fac 925 double precision a,b,c,x,s,avdr,eps,fac 977 926 c Physical constants. 978 927 common /circle/ pi,rad,deg,radius … … 1337 1286 c 1338 1287 include 'locprm.inc' 1339 character*8 ophase, ph,phcd1(maxtt) 1288 logical phasid 1340 1289 c integer stdin, stdout 1341 real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),1342 1 ellip, elcor1343 1290 common /unit/ min, mout, merr 1344 1291 c Physical constants. … … 1352 1299 c Scores. 1353 1300 common /sccom/ scgen(maxph) 1354 data dtds/19.96/1355 1301 c 1356 1302 c stdout = 6 … … 1360 1306 c 1361 1307 navsum=0 1362 c1363 c Initialize depth for the calculation of the traveltimes and partial1364 c derivatives1365 c1366 if (odep .le. 0.0) odep = 1.01367 zs = odep1368 call depset(zs,usrc)1369 c1370 c The first time that ellip is called, it initializes the routine and then1371 c tcor is calculated on subsequent calls. See further down in the1372 c subroutine.1373 c1374 tcor=ellip( 'P ',100.0,zs,glat,45.0)1375 1308 c 1376 1309 c Loop over all stationreadings. … … 1401 1334 enddo 1402 1335 c 1403 6 call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1) 1404 c 1405 c This part of the code is stupid. Since some of the phases read in might not 1406 c be identified by branch (e.g PKPdf), we find the first phase that is closest 1407 c to read in phase. For example, if the phase read in is "P", but the first 1408 c phase in time is "Pn" than the phase is rename to "Pn". 1409 c 1410 do 111 i = 1, n 1411 ph = phcd1(i) 1412 ilen = ntrbl(phase) 1413 if(phase(1:ilen) .eq. ph(1:ilen)) then 1414 travt = tt1(i) 1415 dtdd = dtdd1(i) 1416 dtdh = dtdh1(i) 1417 ophase = phase 1418 phase = ph 1419 go to 100 1336 c ID the phases and compute residuals. 1337 if(phasid()) then 1338 c 1339 c Apply the RRES command. 1340 if(rres.eq.'r'.and.(res.lt.reslim(1).or.res.gt.reslim(2))) then 1341 pflg=' ' 1342 go to 8 1420 1343 endif 1421 111 continue 1422 c 1423 c Phase not found for that distance and depth 1424 c 1425 pflg = ' ' 1426 c write(merr,334)sta,phase 1427 c334 format('setup: ','Sta= ',a5,' no phase computed for ',a8) 1428 go to 8 1429 c 1430 100 continue 1431 c 1432 c Calculate the ellipticity correction for different phases. 1433 c 1434 azimuth = azim*deg 1435 tcor=ellip(phase,delta,zs,glat,azimuth) 1436 if(dtdd.gt.dtds) dtdd=dtds 1437 elcor=(elev/5.57)*sqrt (1.0(5.57*dtdd/(radius*rad))**2) 1438 res=trvtimosectravt(tcor+elcor) 1439 c res0 = res 1440 c 1441 c Apply the RRES command. 1442 if(rres.eq.'r'.and.(res.lt.reslim(1).or.res.gt.reslim(2))) then 1443 pflg=' ' 1444 go to 8 1445 endif 1446 c 1447 ierr=lpush(navsum,j,res) 1448 c write(merr,333)j,navsum,sta,ophase,ph,pflg,delta,azimuth, 1449 c 2 trvtim,res,dtdd,dtdh 1450 c333 format('setup: ',i5,i5,1x,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x, 1451 c 1 f8.2,1x,f6.2,1x,f8.4,1x,f8.4) 1344 c Push the residual for statistical purposes. 1345 ierr=lpush(navsum,j,res) 1346 endif 1452 1347 c 1453 1348 c Restore data to station parameter storage. … … 1481 1376 c 1482 1377 include 'locprm.inc' 1483 character*8 phcd1(maxtt),ph, ophase 1378 logical phasid 1484 1379 c integer stdout,stdin 1485 real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),1486 1 ellip, elcor1487 1380 c Physical constants. 1488 1381 common /circle/ pi,rad,deg,radius … … 1493 1386 c Stationreading parameters. 1494 1387 include 'pickcm.inc' 1495 data dtds/19.96/1496 1388 c 1497 1389 c stdout = 6 … … 1502 1394 navsum=0 1503 1395 c 1504 c Initialize the depth1505 c1506 if (odep .le. 0.0) odep = 1.01507 zs = odep1508 call depset(zs,usrc)1509 1396 c Loop over all stationreadings. 1510 1397 do 1 j=1,nphase … … 1517 1404 1 delta,azim) 1518 1405 c 1519 call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1) 1520 c 1521 c This part of the code is tricky. Since some of the phases read in might not 1522 c be identified by branch (e.g PKPdf), we find the first phase that is closest 1523 c to read in phase. For example, if the phase read in is "P", but the first 1524 c phase in time is "Pn" than the phase is rename to "Pn". 1525 c 1526 do 111 i = 1, n 1527 ph = phcd1(i) 1528 ilen = ntrbl(phase) 1529 if(phase(1:ilen) .eq. ph(1:ilen)) then 1530 travt = tt1(i) 1531 dtdd = dtdd1(i) 1532 dtdh = dtdh1(i) 1533 ophase = phase 1534 phase = ph 1535 go to 100 1536 endif 1537 111 continue 1538 c 1539 c Phase not found for that distance and depth 1540 c 1541 pflg = ' ' 1542 go to 8 1543 c 1544 100 continue 1545 c 1546 c 1547 azimuth = azim*deg 1548 if(dtdd.gt.dtds) dtdd=dtds 1549 elcor=(elev/5.57)*sqrt (1.0(5.57*dtdd/(radius*rad))**2) 1550 tcor=ellip(phase,delta,zs,glat,azim) 1551 res=trvtimosectravt(tcor+elcor) 1552 c 1553 c write(merr,333)j,sta,ophase,ph,pflg,delta,azimuth, 1554 c 2 trvtim,res,dtdd,dtdh 1555 c333 format('stats: ',i5,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x, 1556 c 1 f8.2,1x,f6.2,1x,f8.4,1x,f8.4) 1406 c ID the phase and compute the residual. 1407 if(phasid()) then 1557 1408 c 1558 1409 c Update the statistical sums. 1559 c 1560 ierr=lpush(navsum,j,res) 1410 ierr=lpush(navsum,j,res) 1411 endif 1561 1412 c 1562 1413 c Save computed data in station parameter storage. 1563 1414 c 1564 8call drdwr(2,j,1,5)1415 call drdwr(2,j,1,5) 1565 1416 1 continue 1566 1417 c … … 1654 1505 return 1655 1506 end 1507 1508 c Withers modified declaration of phasid function to prefent f90 1509 c from thinking its main. (gotta love .for) 1510 c logical function phasid 1511 logical function phasid() 1512 1513 c 1514 c Reidentify the current phase and return its traveltime residual. 1515 c 1516 implicit none 1517 real depmin,cn 1518 parameter(depmin=1.,cn=5.) 1519 c Miscellaneous constants. 1520 include 'locprm.inc' 1521 character*8 phcd1(maxtt), tmpcod 1522 character*20 prmmod 1523 logical first 1524 integer*4 n,i,j,ntrbl,mprm,ln 1525 c integer*4 min,mout,merr 1526 real*4 zs,tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt), 1527 1 ellip,elcor,pi,rad,deg,radius,tcor,ecor,zs0,usrc(2), 1528 2 travt(maxtt),win(maxtt),spd,amp 1529 c Physical constants. 1530 common /circle/ pi,rad,deg,radius 1531 c Logical units. 1532 c common /unit/ min,mout,merr 1533 c Hypocentral parameters. 1534 include 'hypocm.inc' 1535 c Stationreading parameters. 1536 include 'pickcm.inc' 1537 data first,zs0/.true.,1./ 1538 c 1539 c The first time through initialize statistical parameters and the 1540 c ellipticity routine. 1541 if(first) then 1542 first=.false. 1543 call freeunit(mprm) 1544 call getprm(mprm,prmmod) 1545 close(mprm) 1546 tcor=ellip('P ',100.0,zs,glat,45.0) 1547 c print *,'Phasid: ellip set' 1548 endif 1549 c 1550 c Initialize depth for the calculation of the traveltimes and partial 1551 c derivatives 1552 zs = amax1(odep,depmin) 1553 if(zs.ne.zs0) then 1554 zs0=zs 1555 call depset(zs,usrc) 1556 c print *,'Phasid: depset  zs glat =',zs,glat 1557 endif 1558 c 1559 c Compute travel times of candidate phases. 1560 c print *,'Phasid: scode depth delta azim elev = ',sta,zs,delta, 1561 c 1 azim*deg,elev 1562 call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1) 1563 c 1564 c Compute and apply traveltime corrections. 1565 do i = 1, n 1566 c print *,'Trtm: phcd tt dtdd dtdh = ',phcd1(i),tt1(i),dtdd1(i), 1567 c 1 dtdh1(i) 1568 ecor=elcor(phcd1(i),elev,dtdd1(i)) 1569 if((phcd1(i).eq.'Pg'.or.phcd1(i).eq.'Pb'.or.phcd1(i).eq.'Pn'.or. 1570 1 phcd1(i).eq.'P').and.dtdh1(i).gt.0.) then 1571 tcor=ellip('Pup',delta,zs,glat,azim*deg) 1572 c print *,'Pup' 1573 else if((phcd1(i).eq.'Sg'.or.phcd1(i).eq.'Sb'.or. 1574 1 phcd1(i).eq.'Sn'.or.phcd1(i).eq.'S').and.dtdh1(i).gt.0.) then 1575 tcor=ellip('Sup',delta,zs,glat,azim*deg) 1576 c print *,'Sup' 1577 else 1578 tcor=ellip(phcd1(i),delta,zs,glat,azim*deg) 1579 endif 1580 travt(i)=tt1(i)+(tcor+ecor)+osec 1581 c Base the association window on half the spread. 1582 call getstt(phcd1(i),delta,dtdd1(i),spd,amp) 1583 win(i)=cn*spd 1584 c print *,'Phasid: phase tcor ecor travt win = ', 1585 c 1 phcd1(i),tcor,ecor,travt(i),win(i) 1586 enddo 1587 c pause 1588 c 1589 c Identify the phase using the Chicxulub algorithm. 1590 c 1591 c Handle a generic P (phase class 'p'). 1592 if(phase.eq.'P') then 1593 i=0 1594 res=1e30 1595 do j=1,n 1596 if(phcd1(j).eq.'Pg'.or.phcd1(j).eq.'Pb'.or.phcd1(j).eq. 1597 1 'Pn'.or.phcd1(j).eq.'P'.or.phcd1(j).eq.'Pdif'.or. 1598 2 phcd1(j)(1:3).eq.'PKP'.or.phcd1(j).eq.'PKiKP') then 1599 if(abs(trvtimtravt(j)).lt.res) then 1600 i=j 1601 res=abs(trvtimtravt(j)) 1602 endif 1603 endif 1604 enddo 1605 c Handle a generic S (phase class 's'). 1606 else if(phase.eq.'S') then 1607 i=0 1608 res=1e30 1609 do j=1,n 1610 if(phcd1(j).eq.'Sg'.or.phcd1(j).eq.'Sb'.or.phcd1(j).eq. 1611 1 'Sn'.or.phcd1(j).eq.'S') then 1612 if(abs(trvtimtravt(j)).lt.res) then 1613 i=j 1614 res=abs(trvtimtravt(j)) 1615 endif 1616 endif 1617 enddo 1618 c Handle a generic PKP (phase class 'k'). 1619 else if(phase.eq.'PKP') then 1620 i=0 1621 res=1e30 1622 do j=1,n 1623 if(phcd1(j)(1:3).eq.'PKP'.or.phcd1(j).eq.'PKiKP') then 1624 if(abs(trvtimtravt(j)).lt.res) then 1625 i=j 1626 res=abs(trvtimtravt(j)) 1627 endif 1628 endif 1629 enddo 1630 c Handle a generic PP (phase class '2'). 1631 else if(phase.eq.'PP') then 1632 i=0 1633 res=1e30 1634 do j=1,n 1635 if(phcd1(j).eq.'PgPg'.or.phcd1(j).eq.'PbPb'.or.phcd1(j).eq. 1636 1 'PnPn'.or.phcd1(j).eq.'PP') then 1637 if(abs(trvtimtravt(j)).lt.res) then 1638 i=j 1639 res=abs(trvtimtravt(j)) 1640 endif 1641 endif 1642 enddo 1643 c Handle a generic core phase (phase class 'c'). 1644 else if(phase.eq.'pPKP'.or.phase.eq.'sPKP'.or.phase.eq.'SKP'.or. 1645 1 phase.eq.'PKKP'.or.phase.eq.'P''P'''.or.phase.eq.'SKS') then 1646 ln=ntrbl(phase) 1647 i=0 1648 res=1e30 1649 do j=1,n 1650 if(phcd1(j)(1:ln).eq.phase) then 1651 if(abs(trvtimtravt(j)).lt.res) then 1652 i=j 1653 res=abs(trvtimtravt(j)) 1654 endif 1655 endif 1656 enddo 1657 c Handle an unidentified phase (phase class ' '). 1658 else if(phase.eq.' ') then 1659 i=0 1660 res=1e30 1661 do j=1,n 1662 c if(abs(trvtimtravt(j)).lt.10.) print *,'phase res win = ', 1663 c 1 phcd1(j),trvtimtravt(j),win(j),delta 1664 if(abs(trvtimtravt(j)).le.win(j)) then 1665 if(abs(trvtimtravt(j)).lt.res) then 1666 i=j 1667 res=abs(trvtimtravt(j)) 1668 endif 1669 endif 1670 enddo 1671 c if(i.gt.0) print *,'Select ',phcd1(i) 1672 c Otherwise, look for an exact match (phase class 'e'). 1673 else 1674 do i=1,n 1675 if(phase.eq.phcd1(i)) go to 1 1676 enddo 1677 i=0 1678 endif 1679 c 1680 c Take care of special cases. 1681 if(i.eq.0.and.phase.ne.' ') then 1682 c If an explicit P crustal phase doesn't exist, let it float. 1683 if(phase.eq.'Pg'.or.phase.eq.'Pb'.or.phase.eq.'Pn') then 1684 i=0 1685 res=1e30 1686 do j=1,n 1687 if(phcd1(j).eq.'Pg'.or.phcd1(j).eq.'Pb'.or.phcd1(j).eq. 1688 1 'Pn') then 1689 if(abs(trvtimtravt(j)).lt.res.and. 1690 1 abs(trvtimtravt(j)).lt.win(j)) then 1691 i=j 1692 res=abs(trvtimtravt(j)) 1693 endif 1694 endif 1695 enddo 1696 c If an explicit S crustal phase doesn't exist, let it float. 1697 else if(phase.eq.'Sg'.or.phase.eq.'Sb'.or.phase.eq.'Sn') then 1698 i=0 1699 res=1e30 1700 do j=1,n 1701 if(phcd1(j).eq.'Sg'.or.phcd1(j).eq.'Sb'.or.phcd1(j).eq. 1702 1 'Sn') then 1703 if(abs(trvtimtravt(j)).lt.res.and. 1704 1 abs(trvtimtravt(j)).lt.win(j)) then 1705 i=j 1706 res=abs(trvtimtravt(j)) 1707 endif 1708 endif 1709 enddo 1710 c Allow natural extensions of P and PKP phases. 1711 else 1712 tmpcod=' ' 1713 if(phase.eq.'Pg') then 1714 tmpcod='Pb' 1715 else if(phase.eq.'Pn') then 1716 tmpcod='P' 1717 else if(phase.eq.'Sg') then 1718 tmpcod='Sb' 1719 else if(phase.eq.'Sn') then 1720 tmpcod='S' 1721 else if(phase.eq.'PKPdf') then 1722 tmpcod='PKiKP' 1723 else if(phase.eq.'PKiKP') then 1724 tmpcod='PKPdf' 1725 else if(phase.eq.'PKPbc') then 1726 tmpcod='PKPdif' 1727 else if(phase.eq.'PKPdif') then 1728 tmpcod='PKPbc' 1729 else if(phase.eq.'Pdif') then 1730 tmpcod='P' 1731 endif 1732 if(tmpcod.ne.' ') then 1733 do i=1,n 1734 if(tmpcod.eq.phcd1(i)) go to 1 1735 enddo 1736 endif 1737 endif 1738 endif 1739 c 1740 c Got it. Compute the residual. 1741 c 1742 1 if(i.gt.0) then 1743 c print *,'Phasid: phase phcd res = ',phase,phcd1(i), 1744 c 1 trvtimtravt(i) 1745 phase=phcd1(i) 1746 res=trvtimtravt(i) 1747 dtdd=dtdd1(i) 1748 dtdh=dtdh1(i) 1749 phasid=.true. 1750 else 1751 c 1752 c Phase not found for that distance and depth 1753 c 1754 pflg = ' ' 1755 res=0. 1756 phasid=.false. 1757 endif 1758 return 1759 end 1760 1761 real*4 function elcor(phase,elev,dtdd) 1762 c 1763 c Return the elevation correction. 1764 c 1765 implicit none 1766 real vp,vs,kmpd 1767 parameter(vp=5.80,vs=3.46,kmpd=111.19493) 1768 character*(*) phase 1769 integer*4 ntrbl,j 1770 real*4 elev,dtdd 1771 c 1772 c Figure out the arriving wave type. 1773 do j=ntrbl(phase),1,1 1774 if(phase(j:j).eq.'P') then 1775 c This phase arrives as a P. Do the correction using Vp. 1776 elcor=(elev/vp)*sqrt(abs(1. 1777 1 amin1((vp*dtdd/kmpd)**2,1.))) 1778 c print *,'Elcor: P corr = ',phase,elcor 1779 return 1780 endif 1781 if(phase(j:j).eq.'S') then 1782 c This phase arrives as an S. Do the correction using Vs. 1783 elcor=(elev/vs)*sqrt(abs(1. 1784 1 amin1((vs*dtdd/kmpd)**2,1.))) 1785 c print *,'Elcor: S corr = ',phase,elcor 1786 return 1787 endif 1788 enddo 1789 return 1790 end 1791 1792 subroutine getprm(nprm,modnam) 1793 c 1794 c Read the parameter file. 1795 c 1796 implicit none 1797 character*(*) modnam 1798 character*20 flnm 1799 integer*4 nprm,i,j,k,ntrbl 1800 include 'resprm.inc' 1801 data flnm/'res_ac.prm'/ 1802 character*1001 dirname 1803 common /my_dir/dirname 1804 c 1805 c Withers modified to be as other open statements in rayloc_ew 1806 c open(nprm,file=flnm,access='sequential',form='formatted', 1807 c 1 status='old',readonly,shared) 1808 c 1809 c Note that this param file is new to this version of rayloc_ew 1810 c (as of 20071205). It should live in the same place as the model 1811 c (e.g. ak135.*). Withers 1812 c 1813 c Make sure we use dirname if one is provided 1814 if(dirname(1:2).eq.'. ') then 1815 open(nprm,file=flnm,access='sequential',form='formatted', 1816 > status='old') 1817 else 1818 open(nprm,file=dirname(1:ntrbl(dirname))//'/'//flnm, 1819 > access='sequential',form='formatted',status='old') 1820 endif 1821 1822 read(nprm,101)modnam 1823 101 format(1x,a) 1824 c print *,'Getprm: modnam = ',modnam 1825 c 1826 do i=1,mxbr 1827 read(nprm,102,end=1)phnm(i),ms(i),me(i),xm(i),cut(i) 1828 102 format(1x,a8,2i8,f8.1,f8.0) 1829 c write(*,102)phnm(i),ms(i),me(i),xm(i),cut(i) 1830 do k=1,3 1831 read(nprm,103)mo(k,i),(prm(j,k,i),j=1,mo(k,i)) 1832 103 format(6x,i5,1p4e15.7:/(11x,1p4e15.7)) 1833 c write(*,103)mo(k,i),(prm(j,k,i),j=1,mo(k,i)) 1834 c pause 1835 enddo 1836 enddo 1837 i=mxbr+1 1838 c 1839 1 nbr=i1 1840 if(nbr.lt.mxbr) phnm(nbr+1)=' ' 1841 c print *,'Getprm: modnam nbr = ',modnam(1:ntrbl(modnam)),nbr 1842 close(nprm) 1843 return 1844 end 1845 1846 subroutine getstt(phcd,delta,dtdd,spd,amp) 1847 c 1848 c Given a phase code, phcd, sourcereceiver distance, delta (in deg), 1849 c and dt/d(delta), dtdd (in s/deg), Getstt returns the observed 1850 c spread, spd (sec), and relative amplitude, amp of the phase at that 1851 c distance. Note that the spread is the robust equivalent of one 1852 c standard deviation. If statistics aren't available, spd and amp are 1853 c set to zero. 1854 c 1855 implicit none 1856 character*(*) phcd 1857 integer*4 j,iupcor,margin 1858 real*4 delta,dtdd,spd,amp,xcor,tcor,x0,evlply 1859 include 'resprm.inc' 1860 c 1861 c Find the desired phase in the statistics common. 1862 do j=1,nbr 1863 if(phcd.eq.phnm(j)) then 1864 c Figure the correction to surface focus. 1865 if(iupcor(phcd(1:1),dtdd,xcor,tcor).le.0) then 1866 c print *,'Iupcor failed on phase ',phcd 1867 xcor=0. 1868 endif 1869 c print *,'Getstt: phcd delta xcor = ',phcd,delta,xcor 1870 c Correct the phase to surface focus. 1871 if(phcd(1:1).eq.'p'.or.phcd(1:1).eq.'s') then 1872 c For depth phases, the deeper, the closer. 1873 x0=deltaxcor 1874 else 1875 c For other phases, the deeper, the farther. 1876 x0=delta+xcor 1877 endif 1878 c If statistics are defined, do the interpolation. 1879 if(phnm(j).ne.phnm(min0(j+1,mxbr)).or.j.ge.nbr) then 1880 margin=15 1881 else 1882 margin=0 1883 endif 1884 if(x0.ge.ms(j)16.and.x0.le.me(j)+margin) then 1885 x0=x0xm(j) 1886 spd=amin1(amax1(evlply(x0,mo(3,j),prm(1,3,j)),.7),10.) 1887 amp=amax1(evlply(x0,mo(1,j),prm(1,1,j)),0.) 1888 c print *,'Getstt: x0 spd amp =',x0,spd,amp,ms(j),me(j) 1889 c pause 1890 return 1891 endif 1892 endif 1893 enddo 1894 c If we don't find anything, return zero. 1895 spd=0. 1896 amp=0. 1897 return 1898 end 1899 1900 function evlply(x,m,b) 1901 c 1902 c Evlply evaluates the polynomial. 1903 c 1904 dimension b(1) 1905 c 1906 xx=1. 1907 evlply=b(1) 1908 do j=2,m 1909 xx=xx*x 1910 evlply=evlply+xx*b(j) 1911 enddo 1912 return 1913 end
Note: See TracChangeset
for help on using the changeset viewer.