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

Revision 1669, 14.3 KB checked in by friberg, 15 years ago (diff)

First commit of rayloc_ew in EW-CENTRAL CVS

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