source: trunk/src/seismic_processing/rayloc_ew/ellip.f @ 3168

Revision 3168, 15.1 KB checked in by withers, 11 years ago (diff)

Synced with hydra version of rayloc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1      real function ellip(phase,edist,edepth,slat,bazim)
2
3c MMW 20071129 merged Hydra and ew version
4
5c IGD 06/01/04 Added dirname common for portability
6 
7C==========================================================================
8C
9C    Ellipticity correction for any given phase using
10C    Dziewonski & Gilbert representation
11C    The ellipticity corrections are found by linear interpolation
12C    in terms of values calculated for the ak135 model for a wide
13C    range of phases to match the output of the iasp software
14C
15C    Input Parameters:
16C    character
17C          phase : a  string specifying the PHASE,   -e.g P, ScP etc.
18C                  'phase' should have at least 8 characters length
19C
20C    real
21C          edist  :  epicentral distance to station (in degrees)
22C          edepth :  depth of event (km)
23C          slat   :  epicentral latitude of source (in degrees)
24C          bazim  :  azimuth from source to station (in degrees)
25C
26C    Output:
27C    real
28C          ellip  :  time correction for path to allow for ellipticity
29C                    (ellip has to be added to the AK135 travel times)
30C                     returns ellip=0. for phases unknown to this software)
31C
32C    Usage:-one call: tcor=ellip(phase,edist,edepth,slat,bazim)
33C           to initialize
34C          -every next call computes ellip from input parameters
35C
36C==========================================================================
37C
38C   B.L.N. Kennett RSES,ANU        May 1995
39C   (based on earlier routine by D.J. Brown)
40 
41C  (modified for processing of large data sets
42C   by W. Spakman, Earth Sciences, Utrecht University, June 1996)
43
44c   Modified logical units for compatibility with RayLocator, July 2003.
45c   Changed lue to merr and allocated lut with freeunit.
46C=========================================================================
47 
48      common /my_dir/dirname
49      character*1001 dirname
50
51      character *(*) phase
52 
53      real edist,edepth,ecolat,bazim,slat,azim,degrad,
54     ^     sc0,sc1,sc2,tcor,
55     ^     tau0, a0,b0,h0,d0,e0,f0,g0,
56     ^     tau1, a1,b1,h1,d1,e1,f1,g1,
57     ^     tau2, a2,b2,h2,d2,e2,f2,g2
58      integer Nd,ios,lut
59 
60c..
61      logical INIT,ECHO1,ECHO2
62 
63c     no screen output if ECHO1=.false.
64c     Only put ECHO1=.true. to familiarize yourself
65c     Only put ECHO2=.true. to get some additional
66c     numbers on output (which only make sense with
67c     the code at hand)
68 
69c..
70      parameter(NUMPH=57)
71 
72c     NUMPH: number of phases supported by the tau-tables.
73c     NUMPH should exactly match the number of phases in the
74c     tau-tables
75c     To add a phase to the software proceed as follows:
76c       . increase NUMPH by 1
77c       . add the phase code at the end of the phcod data statement
78c       . add the proper tau-table entries at the bottom of
79c         file 'tau.tables'
80c       The only syntax check on the tau.table file is that
81c       the order in which the phases appear in the file should be
82c       the same as the order of the phases in the phcod data statement
83 
84      character*8 phdum,phcod(NUMPH)
85      integer phnch(NUMPH),np(NUMPH)
86 
87      parameter(MDEL=50)
88      real dpth(6),delta(NUMPH,MDEL),di1(NUMPH),di2(NUMPH)
89      real t0(NUMPH,MDEL,6),t1(NUMPH,MDEL,6),t2(NUMPH,MDEL,6)
90 
91      common /tau1/ phcod
92      common /tau2/ phnch,np
93      common /tau3/ dpth,delta,di1,di2,t0,t1,t2
94
95        integer minr,mout,merr
96        common /unit/ minr,mout,merr
97 
98      data phcod/
99     &"Pup     ","P       ","Pdiff   ","PKPab   ","PKPbc   ","PKPdf   ",
100     &"PKiKP   ","pP      ","pPKPab  ","pPKPbc  ","pPKPdf  ","pPKiKP  ",
101     &"sP      ","sPKPab  ","sPKPbc  ","sPKPdf  ","sPKiKP  ","PcP     ",
102     &"ScP     ","SKPab   ","SKPbc   ","SKPdf   ","SKiKP   ","PKKPab  ",
103     &"PKKPbc  ","PKKPdf  ","SKKPab  ","SKKPbc  ","SKKPdf  ","PP      ",
104     &"P'P'    ","Sup     ","S       ","Sdiff   ","SKSac   ","SKSdf   ",
105     &"pS      ","pSKSac  ","pSKSdf  ","sS      ","sSKSac  ","sSKSdf  ",
106     &"ScS     ","PcS     ","PKSab   ","PKSbc   ","PKSdf   ","PKKSab  ",
107     &"PKKSbc  ","PKKSdf  ","SKKSac  ","SKKSdf  ","SS      ","S'S'    ",
108     &"SP      ","PS      ","PnS     "/
109 
110c     In addition to the phase names listed above a number of phase aliases
111c     are available in the routine phase_alias, e.g. Pn --> P etc
112c     The input phase code is first checked against the phcod array
113c     and next against the phase aliases.
114 
115      save INIT,ECHO1,ECHO2,/tau1/,/tau2/,/tau3/,Nd,deldst,degrad
116 
117c..   fixed
118      data dpth   / 0.0, 100.0, 200.0, 300.0, 500.0, 700.0 /
119      data degrad / 0.01745329/
120      data Nd     / 6 /
121      data deldst / 5.0 /
122      data INIT   /.true./
123 
124c..   changeable
125      data ECHO1  /.false./
126      data ECHO2  /.false./
127c     data lut    / 10 /
128 
129c...
130 
131      ellip=0.
132 
133      if(INIT) then
134c...   Initialize at first call and return
135        INIT=.false.
136 
137c       check on the length of phase
138        l=len(phase)
139        if(l.lt.8) then
140          write(merr,*)
141     >    'character variable `phase` should have at least length 8'
142          call exit(104)
143        endif
144 
145c       initialize arrays
146        do i=1,NUMPH
147          phnch(i)=lnblk(phcod(i))
148          np(i)     =0
149          di1(i)    =0.
150          di2(i)    =0.
151          do j=1,MDEL
152            delta(i,j)=0.
153            do k=1,6
154              t0(i,j,k)=0.
155              t1(i,j,k)=0.
156              t2(i,j,k)=0.
157            enddo
158          enddo
159        enddo
160 
161c       Read tau.table from unit 'lut'
162 
163        ios=0
164        if(ECHO1)
165     >  write(merr,*) 'ellip: open and read ellipticity tau-table'
166 
167        call freeunit(lut)
168
169c       Make sure we use dirname if one is provided
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 
182        ip=0
18320      continue
184          ip=ip+1
185          if(ip.gt.NUMPH) then
186c           stop reading NUMPH: limit is reached'
187            goto 21
188          endif
189c         phase_code, number_of_epicentral_distance_entries, min_max_dist.
190          read(lut,'(a8,i2,2f10.1)',end=21,iostat=ios)
191     >         phdum,np(ip),di1(ip),di2(ip)
192          if(ios.ne.0) then
193            write(merr,*) 'ellip: read error on phase record err= ',ios,
194     >                 ' on tau-table'
195            call exit(104)
196          endif
197          if(ECHO1) write(merr,'(''reading: '',i3,1x,a8,1x,i2,2f10.1)')
198     >             ip,phdum,np(ip),di1(ip),di2(ip)
199          if(phdum(1:8).ne.phcod(ip)) then
200            write(merr,*)
201     >     'syntax of tau-table does not conform to ',
202     >     'syntax of phase data statement'
203            call exit(104)
204          endif
205          do i=1,np(ip)
206           read(lut,*,iostat=ios,err=999) delta(ip,i)
207           read(lut,*,iostat=ios,err=999) (t0(ip,i,m),m=1,6)
208           read(lut,*,iostat=ios,err=999) (t1(ip,i,m),m=1,6)
209           read(lut,*,iostat=ios,err=999) (t2(ip,i,m),m=1,6)
210          enddo
211        go to 20
212 
21321      continue
214 
215        if(ECHO1) then
216          write(merr,*) 'Number of phases: ',ip-1
217          write(merr,*) 'Phase codes     : '
218 
219c         print out of available phases
220          i=NUMPH/6
221          do j=1,i
222            i1=1+(j-1)*6
223            i2=i1-1+6
224            write(merr,'(7a9)') (phcod(k),k=i1,i2)
225          enddo
226          i=i*6+1
227          j=NUMPH-i
228          if(j.gt.0) then
229            write(merr,'(7a9)') (phcod(k),k=i,NUMPH)
230          endif
231          write(merr,*) ' '
232          write(merr,*)
233     >       'See also the phase aliases in routine phase_alias'
234          write(merr,*) ' '
235        endif
236c.....  End initialization
237        return
238      endif
239 
240c...  set up source dependent constants
241      azim   = bazim*degrad
242      ecolat = (90.0-slat)*degrad
243      call ellref(ecolat,sc0,sc1,sc2)
244 
245      if(ECHO2) then
246        write(merr,*) 'phase,edist,edepth,ecolat,azim'
247        write(merr,*)  phase(1:8),edist,edepth,ecolat,azim
248        write(merr,*) 'sc0,sc1,sc2: ',sc0,sc1,sc2
249      endif
250 
251c...  select phase
252      ip = -1
253      nc=min(lnblk(phase),8)
254      do 10 i=1,NUMPH
255        if(nc.ne.phnch(i)) goto 10
256        if (phase(1:nc) .eq. phcod(i)(1:nc)) then
257          ip = i
258          go to 11
259        endif
260 10   continue
261 11   continue
262 
263      if(ip.eq.-1) then
264c       check phase aliases
265        call phase_alias(phase,edist,ip)
266      endif
267 
268      if(ip.gt.0) then
269        if(ECHO1)
270     >  write(merr,'(''ip:'',i3,'' Selected phase: '',a8,
271     >              '' table distance range: '',2f8.1)')
272     >              ip,phcod(ip),di1(ip),di2(ip)
273      else
274        if(ECHO1)  write(merr,
275     >    '('' Selected phase '',a8,'' is not available'')')
276     >     phase(1:8)
277        ellip=0.
278        return
279      endif
280 
281c...  distance index
282      idist = 1 + int( (edist-di1(ip))/ deldst )
283      if(edist.lt.di1(ip)) idist =1
284      if(edist.gt.di2(ip)) idist= np(ip)-1
285 
286c...  depth index
287      do 25 j = 1,Nd-1
288        if ((dpth(j).le.edepth).and.(dpth(j+1).ge.edepth))then
289          jdepth = j
290          goto 26
291        endif
292 25   continue
293c   Extrapolate below 700 km (& make behavior deterministic).  Added by
294c   R. Buland on 18 March 2004.
295      jdepth=Nd-1
296 26   continue
297 
298      if(ECHO2) write(merr,*) 'idist, jdepth;',idist,jdepth
299 
300c...  Compute tau-values and ellip
301c     need to allow for zero entries (where phase
302c     description strongly depth dependent)
303 
304c tau0
305         a0 = t0(ip,idist,jdepth)
306         b0 = t0(ip,idist,jdepth+1)
307         h0 = t0(ip,idist+1,jdepth+1)
308         d0 = t0(ip,idist+1,jdepth)
309         e0 = a0 +
310     ^       (d0-a0)*(edist-delta(ip,idist))/
311     ^       (delta(ip,idist+1)-delta(ip,idist))
312         f0 = b0 +
313     ^       (h0-b0)*(edist-delta(ip,idist))/
314     ^       (delta(ip,idist+1)-delta(ip,idist))
315         g0 = e0 +
316     ^       (f0-e0)*(edepth-dpth(jdepth))/
317     ^       (dpth(jdepth+1)-dpth(jdepth))
318         tau0 = g0
319c tau1
320         a1 = t1(ip,idist,jdepth)
321         b1 = t1(ip,idist,jdepth+1)
322         h1 = t1(ip,idist+1,jdepth+1)
323         d1 = t1(ip,idist+1,jdepth)
324         e1 = a1 +
325     ^       (d1-a1)*(edist-delta(ip,idist))/
326     ^       (delta(ip,idist+1)-delta(ip,idist))
327         f1 = b1 +
328     ^       (h1-b1)*(edist-delta(ip,idist))/
329     ^       (delta(ip,idist+1)-delta(ip,idist))
330         g1 = e1 +
331     ^       (f1-e1)*(edepth-dpth(jdepth))/
332     ^       (dpth(jdepth+1)-dpth(jdepth))
333         tau1 = g1
334c tau2
335         a2 = t2(ip,idist,jdepth)
336         b2 = t2(ip,idist,jdepth+1)
337         h2 = t2(ip,idist+1,jdepth+1)
338         d2 = t2(ip,idist+1,jdepth)
339         e2 = a2 +
340     ^       (d2-a2)*(edist-delta(ip,idist))/
341     ^       (delta(ip,idist+1)-delta(ip,idist))
342         f2 = b2 +
343     ^       (h2-b2)*(edist-delta(ip,idist))/
344     ^       (delta(ip,idist+1)-delta(ip,idist))
345         g2 = e2 +
346     ^       (f2-e2)*(edepth-dpth(jdepth))/
347     ^       (dpth(jdepth+1)-dpth(jdepth))
348         tau2 = g2
349c
350         caz = cos(azim)
351         cbz = cos(2.0*azim)
352 
353         if(ECHO2) then
354           write(merr,*) 'tau0,tau1,tau2:',tau0,tau1,tau2
355           write(merr,*) 'azim,caz,cbz',azim,caz,cbz
356         endif
357 
358         tcor = sc0*tau0 + sc1*cos(azim)*tau1 + sc2*cos(2.0*azim)*tau2
359 
360      ellip=tcor
361 
362      return
363 
364999   write(merr,*) 'ellip: read error ',ios,' on tau-table'
365      call exit(104)
366998   write(merr,*) 'ellip: open error ',ios,' on tau-table'
367      call exit(104)
368      end
369 
370c--------
371 
372      subroutine ellref(ecolat,sc0,sc1,sc2)
373c...  s3 = sqrt(3.0)/2.0
374      real ecolat,s3,sc0,sc1,sc2
375 
376      data s3 /0.8660254/
377 
378      sc0 = 0.25*(1.0+3.0*cos(2.0*ecolat))
379      sc1 = s3*sin(2.0*ecolat)
380      sc2 = s3*sin(ecolat)*sin(ecolat)
381      return
382      end
383 
384c--------
385 
386      integer function lnblk(s)
387      character *(*) s
388      l = len(s)
389      do i=l,1,-1
390        if (s(i:i).gt.' ') then
391          lnblk=i
392          return
393        endif
394      end do
395      lnblk=0
396      return
397      end
398 
399c------
400 
401      subroutine phase_alias(phase,delta,ip)
402 
403c-    check for alternative phase names
404c     input phase, delta
405c     output ip (index of phcod)
406 
407      character*(*) phase
408 
409      if(phase(1:3).eq.'Pg '.or.phase(1:3).eq.'Pb '.or.
410        1 phase(1:3).eq.'Pn ') then
411c       phase='P       '
412        ip=2
413      else if(phase(1:3).eq.'Sg '.or.phase(1:3).eq.'Sb '.or.
414        1 phase(1:3).eq.'Sn ') then
415c       phase='S       '
416        ip=33
417      else if(phase(1:4).eq.'pPg '.or.phase(1:4).eq.'pPb '.or.
418        1 phase(1:4).eq.'pPn ') then
419c       phase='pP      '
420        ip=8
421      else if(phase(1:4).eq.'pwP ') then
422c       phase='pP      '
423        ip=8
424      else if(phase(1:5).eq.'pwPg '.or.phase(1:5).eq.'pwPb '.or.
425        1 phase(1:5).eq.'pwPn ') then
426c       phase='pP      '
427        ip=8
428      else if(phase(1:4).eq.'sPg '.or.phase(1:4).eq.'sPb '.or.
429        1 phase(1:4).eq.'sPn ') then
430c       phase='sP      '
431        ip=13
432      else if(phase(1:4).eq.'pSg '.or.phase(1:4).eq.'pSb '.or.
433        1 phase(1:4).eq.'pSn ') then
434c       phase='pS      '
435        ip=37
436      else if(phase(1:4).eq.'sSg '.or.phase(1:4).eq.'sSb '.or.
437        1 phase(1:4).eq.'sSn ') then
438c       phase='sS      '
439        ip=40
440      else if(phase(1:4).eq.'SPg '.or.phase(1:4).eq.'SPb '.or.
441        1 phase(1:4).eq.'SPn ') then
442c       phase='SP      '
443        ip=55
444      else if(phase(1:4).eq.'SgP '.or.phase(1:4).eq.'SbP '.or.
445        1 phase(1:4).eq.'SnP ') then
446c       phase='SP      '
447        ip=55
448      else if(phase(1:4).eq.'PSg '.or.phase(1:4).eq.'PSb '.or.
449        1 phase(1:4).eq.'PSn ') then
450c       phase='PS      '
451        ip=56
452      else if(phase(1:5).eq.'PgPg '.or.phase(1:5).eq.'PbPb '.or.
453        1 phase(1:5).eq.'PnPn ') then
454c       phase='PP      '
455        ip=30
456      else if(phase(1:5).eq.'SgSg '.or.phase(1:5).eq.'SbSb '.or.
457        1 phase(1:5).eq.'SnSn ') then
458c       phase='SS      '
459        ip=53
460      else if(phase(1:2).eq.'p ') then
461c       phase='Pup     '
462        ip=1
463      else if(phase(1:2).eq.'s ') then
464c       phase='Sup     '
465        ip=32
466      elseif(phase(1:6).eq."P'P'ab") then
467c       phase="P'P'    '
468        ip=31
469      else if(phase(1:6).eq."P'P'bc") then
470c       phase="P'P'    '
471        ip=31
472      else if(phase(1:6).eq."P'P'df") then
473c       phase="P'P'    '
474        ip=31
475      else if(phase(1:6).eq."S'S'ac") then
476c       phase="S'S'    '
477        ip=54
478      else if(phase(1:6).eq."S'S'df") then
479c       phase="S'S'    '
480        ip=54
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
488c       phase='pP      '
489        ip=8
490      else if(delta.le.100.0.and.phase.eq.'pwPdif  ') then
491c       phase='pP      '
492        ip=8
493      else if(delta.le.100.0.and.phase.eq.'sPdif   ') then
494c       phase='sP      '
495        ip=13
496      else if(delta.le.100.0.and.phase.eq.'pSdif   ') then
497c       phase='pS      '
498        ip=37
499      else if(delta.le.100.0.and.phase.eq.'sSdif   ') then
500c       phase='sS      '
501        ip=40
502      else if(delta.le.165.0.and.phase.eq.'PKPdif  ') then
503c       phase='PKPbc '
504        ip=5
505      else if(delta.le.165.0.and.phase.eq.'pPKPdif ') then
506c       phase='pPKPbc '
507        ip=10
508      else if(delta.le.165.0.and.phase.eq.'sPKPdif ') then
509c       phase='sPKPbc
510        ip=15
511      else
512        ip=-1
513      endif
514      return
515      end
Note: See TracBrowser for help on using the repository browser.