Changeset 3168


Ignore:
Timestamp:
12/05/07 11:13:22 (11 years ago)
Author:
withers
Message:

Synced with hydra version of rayloc

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 
    12Files rayloc_message_rw.h and rayloc_message_rw.c contain 
    23API for reading and writing TYPEi_RAYLOC messages. 
     
    1415Ilya Dricker (i.dricker@isti.com) 
    1516and ISTI team 
     17 
     1820071205 update from Mitch Withers.  Synced rayloc_ew with ray_loc 
     19from Hydra.  Involved changes to ellip.f and robust_util.f.  Also 
     20needs a new prm file that I had previously allowed to reside in 
     21the glass model dir.  Added that dir and its contents along with 
     22the new res_ac.prm file to rayloc_ew's directory.  There are also 
     23problems with the hydra ew that should be fixed. 
     24 
     25A known feature is rayloc will complain with a less than illuminating 
     26error when a station is missing from the station file.  Something like: 
     27 
     2820071205_UTC_15:12:50 lib_rayloc: Failed to write Input Structure to the file /g 
     29aia/home/mwithers/Projects/EWSupport/Working/run_smeagol/model/RayLocInputa.txt 
     3020071205_UTC_15:12:50 rayloc_ew: lib_rayloc FAILED: retVal = -1 
     3120071205_UTC_15:12:50 rayloc_ew: Output Message size = NULL 
     3220071205_UTC_15:12:50 ERROR: rayloc_MessageToRaylocHeader: RAYLOC HEADER MESSAGE 
     33 is NULL 
     3420071205_UTC_15:12:50 Failed to convert RAYLOC message into structure 
     35 
     36You have been warned. 
  • trunk/src/seismic_processing/rayloc_ew/ellip.f

    r1669 r3168  
    11      real function ellip(phase,edist,edepth,slat,bazim) 
    22 
    3 c IGD 06/01/04 Added dirname common for portability  
    4  
     3c MMW 20071129 merged Hydra and ew version 
     4 
     5c IGD 06/01/04 Added dirname common for portability 
     6  
    57C========================================================================== 
    68C 
     
    4345c   Changed lue to merr and allocated lut with freeunit. 
    4446C========================================================================= 
    45         common /my_dir/dirname 
    46         character*1001 dirname 
    47   
     47  
     48      common /my_dir/dirname 
     49      character*1001 dirname 
     50 
    4851      character *(*) phase 
    49  
     52  
    5053      real edist,edepth,ecolat,bazim,slat,azim,degrad, 
    5154     ^     sc0,sc1,sc2,tcor, 
     
    6568  
    6669c.. 
    67  
    6870      parameter(NUMPH=57) 
    6971  
     
    112114  
    113115      save INIT,ECHO1,ECHO2,/tau1/,/tau2/,/tau3/,Nd,deldst,degrad 
    114  
    115116  
    116117c..   fixed 
     
    167168 
    168169c       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 
     178c       open(lut,file='tau.table',form='formatted', 
     179c    >           iostat=ios,err=998,status='old',readonly,shared) 
     180        rewind(lut) 
     181  
    177182        ip=0 
    17818320      continue 
     
    207212  
    20821321      continue 
    209  
    210 c       IGD 06/18/04 Closed statemet included          
    211         close(unit=lut) 
     214  
    212215        if(ECHO1) then 
    213216          write(merr,*) 'Number of phases: ',ip-1 
     
    262265        call phase_alias(phase,edist,ip) 
    263266      endif 
    264  
     267  
    265268      if(ip.gt.0) then 
    266269        if(ECHO1) 
     
    275278        return 
    276279      endif 
    277  
     280  
    278281c...  distance index 
    279282      idist = 1 + int( (edist-di1(ip))/ deldst ) 
    280283      if(edist.lt.di1(ip)) idist =1 
    281284      if(edist.gt.di2(ip)) idist= np(ip)-1 
    282  
     285  
    283286c...  depth index 
    284 c     IGD 07/29/04 Next line to prevent crashes if jdepth is undefined 
    285 c     IGD 07/29/04 It is not clear to me how much we screw up processing 
    286       jdepth = Nd 
    287287      do 25 j = 1,Nd-1 
    288288        if ((dpth(j).le.edepth).and.(dpth(j+1).ge.edepth))then 
     
    291291        endif 
    292292 25   continue 
     293c   Extrapolate below 700 km (& make behavior deterministic).  Added by  
     294c   R. Buland on 18 March 2004. 
     295      jdepth=Nd-1 
    293296 26   continue 
    294  
     297  
    295298      if(ECHO2) write(merr,*) 'idist, jdepth;',idist,jdepth 
    296  
     299  
    297300c...  Compute tau-values and ellip 
    298301c     need to allow for zero entries (where phase 
    299302c     description strongly depth dependent) 
    300  
     303  
    301304c tau0 
    302305         a0 = t0(ip,idist,jdepth) 
     
    356359  
    357360      ellip=tcor 
     361  
    358362      return 
    359363  
     
    403407      character*(*) phase 
    404408  
    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 
    406411c       phase='P       ' 
    407412        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 
    409415c       phase='S       ' 
    410416        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 
    412419c       phase='pP      ' 
    413420        ip=8 
     
    415422c       phase='pP      ' 
    416423        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 
    418426c       phase='pP      ' 
    419427        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 
    421430c       phase='sP      ' 
    422431        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 
    424434c       phase='pS      ' 
    425435        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 
    427438c       phase='sS      ' 
    428439        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 
    430442c       phase='SP      ' 
    431443        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 
    433446c       phase='SP      ' 
    434447        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 
    436450c       phase='PS      ' 
    437451        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 
    439454c       phase='PP      ' 
    440455        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 
    442458c       phase='SS      ' 
    443459        ip=53 
     
    463479c       phase="S'S'    ' 
    464480        ip=54 
    465       else if(delta.le.100.0.and.phase.eq.'pPdiff  ') then 
     481        else if(phase.eq.'Pdif') then 
     482c         phase='Pdiff' 
     483          ip=3 
     484        else if(phase.eq.'Sdif') then 
     485c         phase='Sdiff' 
     486          ip=34 
     487      else if(delta.le.100.0.and.phase.eq.'pPdif   ') then 
    466488c       phase='pP      ' 
    467489        ip=8 
    468       else if(delta.le.100.0.and.phase.eq.'pwPdiff ') then 
     490      else if(delta.le.100.0.and.phase.eq.'pwPdif  ') then 
    469491c       phase='pP      ' 
    470492        ip=8 
    471       else if(delta.le.100.0.and.phase.eq.'sPdiff  ') then 
     493      else if(delta.le.100.0.and.phase.eq.'sPdif   ') then 
    472494c       phase='sP      ' 
    473495        ip=13 
    474       else if(delta.le.100.0.and.phase.eq.'pSdiff  ') then 
     496      else if(delta.le.100.0.and.phase.eq.'pSdif   ') then 
    475497c       phase='pS      ' 
    476498        ip=37 
    477       else if(delta.le.100.0.and.phase.eq.'sSdiff  ') then 
     499      else if(delta.le.100.0.and.phase.eq.'sSdif   ') then 
    478500c       phase='sS      ' 
    479501        ip=40 
    480       else if(delta.le.165.0.and.phase.eq.'PKPdiff ') then 
     502      else if(delta.le.165.0.and.phase.eq.'PKPdif  ') then 
    481503c       phase='PKPbc ' 
    482504        ip=5 
    483       else if(delta.le.165.0.and.phase.eq.'pPKPdiff') then 
     505      else if(delta.le.165.0.and.phase.eq.'pPKPdif ') then 
    484506c       phase='pPKPbc ' 
    485507        ip=10 
    486       else if(delta.le.165.0.and.phase.eq.'sPKPdiff') then 
     508      else if(delta.le.165.0.and.phase.eq.'sPKPdif ') then 
    487509c       phase='sPKPbc 
    488510        ip=15 
  • trunk/src/seismic_processing/rayloc_ew/input.f

    r1669 r3168  
    7979c 
    8080c   Reset to generic phase codes so they will be used and reidentified. 
    81         if(phase.eq.'Pg'.or.phase.eq.'Pn'.or.phase.eq.'Pb'.or. 
    82         1       phase.eq.'P*') phase='P' 
    83         if(phase.eq.'Sg'.or.phase.eq.'Sn'.or.phase.eq.'Sb'.or. 
    84         1       phase.eq.'S*') phase='S' 
     81c       if(phase.eq.'Pg'.or.phase.eq.'Pn'.or.phase.eq.'Pb'.or. 
     82c       1       phase.eq.'P*') phase='P' 
     83c       if(phase.eq.'Sg'.or.phase.eq.'Sn'.or.phase.eq.'Sb'.or. 
     84c       1       phase.eq.'S*') phase='S' 
    8585c 
    8686c Check to see if you should use the phase in the earthquake location 
  • trunk/src/seismic_processing/rayloc_ew/libtau.f

    r1669 r3168  
    887887      if(ln.le.0) ln=index(phcd(jb),' ')-1 
    888888      if(ln.le.0) ln=len(phcd(jb)) 
    889       phnm(n)=phcd(jb)(1:ln)//'diff' 
     889      phnm(n)=phcd(jb)(1:ln)//'dif' 
    890890 10   continue 
    891891      return 
  • trunk/src/seismic_processing/rayloc_ew/makefile.sol

    r2055 r3168  
    77#    Revision history: 
    88#     $Log$ 
     9#     Revision 1.3  2007/12/05 16:13:22  withers 
     10#     Synced with hydra version of rayloc 
     11# 
    912#     Revision 1.2  2006/01/19 19:04:55  friberg 
    1013#     rayloc upgraded to use a library version of rayloc_message_rw.c 
     
    4144 
    4245CFLAGS = ${GLOBALFLAGS} -g -I ${I} -DUSE_LOGIT 
    43 FFLAGS = -g 
     46FFLAGS = -g -f77=tab 
    4447B = $(EW_HOME)/$(EW_VERSION)/bin 
    4548L = $(EW_HOME)/$(EW_VERSION)/lib 
     
    8184 
    8285rayloc_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  
    8487 
    8588$(TARGET1):         $(F_OBJS) $(F_MAIN) 
  • trunk/src/seismic_processing/rayloc_ew/robust_util.f

    r1669 r3168  
    206206c 
    207207      include 'locprm.inc' 
    208         character*21 flgs,cstor 
     208        character*22 flgs,cstor 
    209209c   Logical unit numbers. 
    210210      common /unit/ min,mout,merr 
     
    806806c 
    807807      include 'locprm.inc' 
    808       character*8 phcd1(maxtt),ophase, ph 
     808        logical phasid,log 
    809809c       integer stdin, stdout 
    810       real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2), 
    811      1 ellip, elcor 
    812810      double precision a,b,c,x,s,avdr 
    813811c   Physical constants. 
     
    821819c   Station-reading parameters. 
    822820      include 'pickcm.inc' 
    823         data dtds/19.96/ 
    824821c 
    825822c       stdin = 5 
     
    8498461     ierr=normeq(navsum,ihold) 
    850847c 
    851 c Initialize the depth 
    852 c 
    853 2       if (odep .le. 0.0) odep = 1.0 
    854         zs = odep 
    855       call depset(zs,usrc) 
    856 c 
    857848c   Loop over all station-readings. 
    858849c 
    859       do 3 j = 1, nphase 
     8502     do 3 j = 1, nphase 
    860851c 
    861852c   Fetch data from station parameter storage. 
     
    872863        1                       delta,azim) 
    873864c 
    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=trvtim-osec-travt-(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() 
    916866c 
    917867c   All done.  Push data back out to station parameter storage. 
    918868c 
    919 5     call drdwr(2,j,1,13) 
     869      call drdwr(2,j,1,13) 
    9208703     continue 
    921871c 
     
    973923c   by E. R. Engdahl. 
    974924c 
    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 
    977926c   Physical constants. 
    978927      common /circle/ pi,rad,deg,radius 
     
    13371286c 
    13381287      include 'locprm.inc' 
    1339       character*8 ophase, ph,phcd1(maxtt) 
     1288        logical phasid 
    13401289c       integer stdin, stdout 
    1341       real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2), 
    1342      1 ellip, elcor 
    13431290        common /unit/ min, mout, merr 
    13441291c   Physical constants. 
     
    13521299c   Scores. 
    13531300      common /sccom/ scgen(maxph) 
    1354         data dtds/19.96/ 
    13551301c 
    13561302c       stdout = 6 
     
    13601306c 
    13611307      navsum=0 
    1362 c 
    1363 c   Initialize depth for the calculation of the traveltimes and partial 
    1364 c      derivatives 
    1365 c 
    1366       if (odep .le. 0.0) odep = 1.0 
    1367         zs = odep 
    1368       call depset(zs,usrc) 
    1369 c 
    1370 c The first time that ellip is called, it initializes the routine and then 
    1371 c   tcor is calculated on subsequent calls.  See further down in the 
    1372 c   subroutine. 
    1373 c 
    1374         tcor=ellip( 'P       ',100.0,zs,glat,45.0) 
    13751308c 
    13761309c   Loop over all station-readings. 
     
    14011334      enddo 
    14021335c 
    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 
     1336c   ID the phases and compute residuals. 
     1337      if(phasid()) then 
     1338c 
     1339c   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 
    14201343                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=trvtim-osec-travt-(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) 
     1344c   Push the residual for statistical purposes. 
     1345                ierr=lpush(navsum,j,res) 
     1346        endif 
    14521347c 
    14531348c   Restore data to station parameter storage. 
     
    14811376c 
    14821377      include 'locprm.inc' 
    1483       character*8 phcd1(maxtt),ph, ophase 
     1378        logical phasid 
    14841379c       integer stdout,stdin 
    1485       real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2), 
    1486      1 ellip, elcor 
    14871380c   Physical constants. 
    14881381      common /circle/ pi,rad,deg,radius 
     
    14931386c   Station-reading parameters. 
    14941387      include 'pickcm.inc' 
    1495         data dtds/19.96/ 
    14961388c 
    14971389c       stdout = 6 
     
    15021394      navsum=0 
    15031395c 
    1504 c Initialize the depth 
    1505 c 
    1506         if (odep .le. 0.0) odep = 1.0 
    1507         zs = odep 
    1508       call depset(zs,usrc) 
    15091396c   Loop over all station-readings. 
    15101397      do 1 j=1,nphase 
     
    15171404        1                       delta,azim) 
    15181405c 
    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=trvtim-osec-travt-(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) 
     1406c   ID the phase and compute the residual. 
     1407        if(phasid()) then 
    15571408c 
    15581409c   Update the statistical sums. 
    1559 c 
    1560       ierr=lpush(navsum,j,res) 
     1410                ierr=lpush(navsum,j,res) 
     1411        endif 
    15611412c 
    15621413c   Save computed data in station parameter storage. 
    15631414c 
    1564 8     call drdwr(2,j,1,5) 
     1415      call drdwr(2,j,1,5) 
    156514161       continue 
    15661417c 
     
    16541505      return 
    16551506      end 
     1507 
     1508c   Withers modified declaration of phasid function to prefent f90 
     1509c   from thinking its main.  (gotta love .for) 
     1510c     logical function phasid 
     1511      logical function phasid() 
     1512 
     1513c 
     1514c   Re-identify the current phase and return its travel-time residual. 
     1515c 
     1516      implicit none 
     1517        real depmin,cn 
     1518        parameter(depmin=1.,cn=5.) 
     1519c   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 
     1525c       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 
     1529c   Physical constants. 
     1530      common /circle/ pi,rad,deg,radius 
     1531c   Logical units. 
     1532c     common /unit/ min,mout,merr 
     1533c   Hypocentral parameters. 
     1534      include 'hypocm.inc' 
     1535c   Station-reading parameters. 
     1536      include 'pickcm.inc' 
     1537        data first,zs0/.true.,-1./ 
     1538c 
     1539c   The first time through initialize statistical parameters and the  
     1540c   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) 
     1547c               print *,'Phasid: ellip set' 
     1548        endif 
     1549c 
     1550c   Initialize depth for the calculation of the traveltimes and partial 
     1551c      derivatives 
     1552        zs = amax1(odep,depmin) 
     1553        if(zs.ne.zs0) then 
     1554                zs0=zs 
     1555                call depset(zs,usrc) 
     1556c               print *,'Phasid: depset - zs glat =',zs,glat 
     1557        endif 
     1558c 
     1559c   Compute travel times of candidate phases. 
     1560c       print *,'Phasid: scode depth delta azim elev = ',sta,zs,delta, 
     1561c       1 azim*deg,elev 
     1562      call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1) 
     1563c 
     1564c   Compute and apply travel-time corrections. 
     1565        do i = 1, n 
     1566c               print *,'Trtm: phcd tt dtdd dtdh = ',phcd1(i),tt1(i),dtdd1(i), 
     1567c       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) 
     1572c                       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) 
     1576c                       print *,'Sup' 
     1577                else 
     1578                        tcor=ellip(phcd1(i),delta,zs,glat,azim*deg) 
     1579                endif 
     1580                travt(i)=tt1(i)+(tcor+ecor)+osec 
     1581c   Base the association window on half the spread. 
     1582          call getstt(phcd1(i),delta,dtdd1(i),spd,amp) 
     1583          win(i)=cn*spd 
     1584c               print *,'Phasid: phase tcor ecor travt win = ', 
     1585c       1        phcd1(i),tcor,ecor,travt(i),win(i) 
     1586        enddo 
     1587c       pause 
     1588c 
     1589c   Identify the phase using the Chicxulub algorithm. 
     1590c 
     1591c   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(trvtim-travt(j)).lt.res) then 
     1600              i=j 
     1601              res=abs(trvtim-travt(j)) 
     1602            endif 
     1603          endif 
     1604        enddo 
     1605c   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(trvtim-travt(j)).lt.res) then 
     1613              i=j 
     1614              res=abs(trvtim-travt(j)) 
     1615            endif 
     1616          endif 
     1617        enddo 
     1618c   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(trvtim-travt(j)).lt.res) then 
     1625              i=j 
     1626              res=abs(trvtim-travt(j)) 
     1627            endif 
     1628          endif 
     1629        enddo 
     1630c   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(trvtim-travt(j)).lt.res) then 
     1638              i=j 
     1639              res=abs(trvtim-travt(j)) 
     1640            endif 
     1641          endif 
     1642        enddo 
     1643c   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(trvtim-travt(j)).lt.res) then 
     1652              i=j 
     1653              res=abs(trvtim-travt(j)) 
     1654            endif 
     1655          endif 
     1656        enddo 
     1657c   Handle an unidentified phase (phase class ' '). 
     1658      else if(phase.eq.' ') then 
     1659        i=0 
     1660        res=1e30 
     1661        do j=1,n 
     1662c         if(abs(trvtim-travt(j)).lt.10.) print *,'phase res win = ', 
     1663c    1     phcd1(j),trvtim-travt(j),win(j),delta 
     1664          if(abs(trvtim-travt(j)).le.win(j)) then 
     1665            if(abs(trvtim-travt(j)).lt.res) then 
     1666              i=j 
     1667              res=abs(trvtim-travt(j)) 
     1668            endif 
     1669          endif 
     1670        enddo 
     1671c       if(i.gt.0) print *,'Select ',phcd1(i) 
     1672c   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 
     1679c 
     1680c   Take care of special cases. 
     1681      if(i.eq.0.and.phase.ne.' ') then 
     1682c   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(trvtim-travt(j)).lt.res.and. 
     1690        1                abs(trvtim-travt(j)).lt.win(j)) then 
     1691                i=j 
     1692                res=abs(trvtim-travt(j)) 
     1693              endif 
     1694            endif 
     1695          enddo 
     1696c   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(trvtim-travt(j)).lt.res.and. 
     1704        1                abs(trvtim-travt(j)).lt.win(j)) then 
     1705                i=j 
     1706                res=abs(trvtim-travt(j)) 
     1707              endif 
     1708            endif 
     1709          enddo 
     1710c   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 
     1739c 
     1740c   Got it.  Compute the residual. 
     1741c 
     17421       if(i.gt.0) then 
     1743c               print *,'Phasid: phase phcd res = ',phase,phcd1(i), 
     1744c       1        trvtim-travt(i) 
     1745                phase=phcd1(i) 
     1746                res=trvtim-travt(i) 
     1747                dtdd=dtdd1(i) 
     1748                dtdh=dtdh1(i) 
     1749                phasid=.true. 
     1750        else 
     1751c 
     1752c Phase not found for that distance and depth 
     1753c 
     1754                pflg = ' ' 
     1755                res=0. 
     1756                phasid=.false. 
     1757        endif 
     1758        return 
     1759        end 
     1760 
     1761        real*4 function elcor(phase,elev,dtdd) 
     1762c 
     1763c   Return the elevation correction. 
     1764c 
     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 
     1771c 
     1772c   Figure out the arriving wave type. 
     1773        do j=ntrbl(phase),1,-1 
     1774                if(phase(j:j).eq.'P') then 
     1775c   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.))) 
     1778c                       print *,'Elcor: P corr = ',phase,elcor 
     1779                        return 
     1780                endif 
     1781                if(phase(j:j).eq.'S') then 
     1782c   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.))) 
     1785c                       print *,'Elcor: S corr = ',phase,elcor 
     1786                        return 
     1787                endif 
     1788        enddo 
     1789      return 
     1790        end 
     1791 
     1792      subroutine getprm(nprm,modnam) 
     1793c 
     1794c   Read the parameter file. 
     1795c 
     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 
     1804c 
     1805c    Withers modified to be as other open statements in rayloc_ew 
     1806c     open(nprm,file=flnm,access='sequential',form='formatted', 
     1807c    1 status='old',readonly,shared) 
     1808c 
     1809c    Note that this param file is new to this version of rayloc_ew 
     1810c    (as of 20071205).  It should live in the same place as the model 
     1811c    (e.g. ak135.*).  Withers 
     1812c 
     1813c       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) 
     1824c     print *,'Getprm:  modnam = ',modnam 
     1825c 
     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) 
     1829c        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)) 
     1833c           write(*,103)mo(k,i),(prm(j,k,i),j=1,mo(k,i)) 
     1834c             pause 
     1835         enddo 
     1836      enddo 
     1837      i=mxbr+1 
     1838c 
     1839 1    nbr=i-1 
     1840      if(nbr.lt.mxbr) phnm(nbr+1)=' ' 
     1841c     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) 
     1847c 
     1848c   Given a phase code, phcd, source-receiver distance, delta (in deg),  
     1849c   and dt/d(delta), dtdd (in s/deg), Getstt returns the observed  
     1850c   spread, spd (sec), and relative amplitude, amp of the phase at that  
     1851c   distance.  Note that the spread is the robust equivalent of one  
     1852c   standard deviation.  If statistics aren't available, spd and amp are  
     1853c   set to zero. 
     1854c 
     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' 
     1860c 
     1861c   Find the desired phase in the statistics common. 
     1862      do j=1,nbr 
     1863        if(phcd.eq.phnm(j)) then 
     1864c   Figure the correction to surface focus. 
     1865          if(iupcor(phcd(1:1),dtdd,xcor,tcor).le.0) then 
     1866c           print *,'Iupcor failed on phase ',phcd 
     1867            xcor=0. 
     1868          endif 
     1869c         print *,'Getstt:  phcd delta xcor = ',phcd,delta,xcor 
     1870c   Correct the phase to surface focus. 
     1871          if(phcd(1:1).eq.'p'.or.phcd(1:1).eq.'s') then 
     1872c   For depth phases, the deeper, the closer. 
     1873            x0=delta-xcor 
     1874          else 
     1875c   For other phases, the deeper, the farther. 
     1876            x0=delta+xcor 
     1877          endif 
     1878c   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=x0-xm(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.) 
     1888c           print *,'Getstt:  x0 spd amp =',x0,spd,amp,ms(j),me(j) 
     1889c           pause 
     1890            return 
     1891          endif 
     1892        endif 
     1893      enddo 
     1894c   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) 
     1901c 
     1902c   Evlply evaluates the polynomial. 
     1903c 
     1904      dimension b(1) 
     1905c 
     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.