program extracticeman

c     Code returns the P- and S-wave speed for the mantle beneath Iceland
c     given r (radius), theta (lat), and phi (lon).  
c     The velocities are extracted from the ICEMAN models, for more info 
c     see Allen et al, Imaging the mantle beneath Iceland using integrated
c     seismological techniques, JGR in press 2002.

c     Richard M Allen (rallen@geology.wisc.edu) - May 2002

c     INPUT
c
c     You will need the following files: 
c     extracticeman.f (this file) ICEMAN-HP ICEMAN-LP ICEMAN-S model.cube.iceman
c
c     Code opens the file "icein" which should specify the velocity models you 
c     would like to use followed by a list of r, theta, phi:
c     1st line: ICEMAN-HP or ICEMAN-LP
c     2rd line: plumehead or null
c     3rd line: absolute or anomaly or percent
c     4th line: r, theta, phi
c     5th line: r, theta, phi
c     etc
c
c     The first line specifies if you want to use a mantly P-velocity model:
c         ICEMAN-HP	- P-velocity model using ~1 Hz compressional phase traveltimes
c         ICEMAN-LP	- P-velocity model using 10-30 sec compressional phases 
c     The second line specifies if you want to add the plumehead portion of the model
c         null 		- no plumehead will be added
c         plumehead 	- the horizontal low velocity anomaly in the upper 200 km
c                         of the mantle beneath Iceland will be added TO THE S-VEL ONLY. 
c     The third line specifies how your want the velocity anomaly returned
c         absolute	- absolute velocity (= ICEMAN anomaly + IASPI91) - km/s
c         anomaly	- the ICEMAN velocity anomaly - km/s
c         percent	- the ICEMAN velocity anomaly as a percentage (= ICEMAN/IASPI91)

c     OUTPUT
c     The returned velocity values are written to a file iceout.  The file contains:
c     r, theta, phi, P-vel, S-vel

c     EXAMPLE - input file: icein
c	 ICEMAN-LP
c	 null
c	 percent
c	 5920.00   64.70  -17.50
c	 5970.00   64.70  -17.50
c	 6020.00   64.70  -17.50
c	 6070.00   64.70  -17.50
c	 6120.00   64.70  -17.50
c	 6170.00   64.70  -17.50
c	 6220.00   64.70  -17.50
c	 6270.00   64.70  -17.50
c	 6320.00   64.70  -17.50
c	 6370.00   64.70  -17.50
c
c     EXAMPLE - output file: iceout
c	 r theta phi Vp Vs
c	  5920.  64.7 -17.5  0.  0.
c	  5970.  64.7 -17.5  0.  0.
c	  6020.  64.7 -17.5 -1.79484766 -2.67779002
c	  6070.  64.7 -17.5 -1.57987272 -3.22103154
c	  6120.  64.7 -17.5 -1.21828759 -3.26796325
c	  6170.  64.7 -17.5 -0.794649483 -3.12111122
c	  6220.  64.7 -17.5 -0.00886660355 -2.30398237
c	  6270.  64.7 -17.5  1.29382737 -0.0660613834
c	  6320.  64.7 -17.5  2.0659201  1.70112513
c	  6370.  64.7 -17.5  0.  0.


      implicit none
      logical infile,pmod,smod,cmod,plumehead,abs,pct
      character*40 manp,mans,addph,mc,absanom
      double precision r,th,ph,rr,rth,rph,pvs,pvp,prho,torad
      double precision dvelicav,velicav
      real realrad,iaspivp,iaspivs,iaspirho
      integer l

      write(*,*)
      write(*,*)'EXTRACTICEMAN: '
      write(*,*)'  will extract Vp and Vs from the ICEMAN vel models'
      write(*,*)

ccc   initialize
      mc='model.cube.iceman'
      mans='ICEMAN-S'
      torad=3.141592654/180.

c
ccc   check all necessary files are present
c
      infile=.FALSE.
      inquire(file='icein',exist=infile)
      if (infile) then
        write(*,*)'found: icein'
      else
        write(*,*)'ERROR: can not find: icein'
        stop
      endif
      infile=.FALSE.
      inquire(file='ICEMAN-HP',exist=infile)
      if (infile) then
        write(*,*)'found: ICEMAN-HP'
      else
        write(*,*)'ERROR: can not find: ICEMAN-HP'
        stop
      endif
      infile=.FALSE.
      inquire(file='ICEMAN-LP',exist=infile)
      if (infile) then
        write(*,*)'found: ICEMAN-LP'
      else
        write(*,*)'ERROR: can not find: ICEMAN-LP'
        stop
      endif
      infile=.FALSE.
      inquire(file='ICEMAN-S',exist=infile)
      if (infile) then
        write(*,*)'found: ICEMAN-S'
      else
        write(*,*)'ERROR: can not find: ICEMAN-S'
        stop
      endif
      infile=.FALSE.
      inquire(file='model.cube.iceman',exist=infile)
      if (infile) then
        write(*,*)'found: model.cube.iceman'
      else
        write(*,*)'ERROR: can not find: model.cube.iceman'
        stop
      endif
      write(*,*)'found all necessary input files, continuing...'
      write(*,*)

c
ccc   open input file, read chosen models
c
      open(unit=1, file='icein',status='old')
      read(1,*)manp
      read(1,*)addph
      read(1,*)absanom
      if (manp .eq. 'ICEMAN-HP' .OR. manp .eq. 'ICEMAN-LP') then
        call read_plume_model(manp,mans,mc)
        write(*,*)'read in all the model files, continuing...'
        write(*,*)
      else
        write(*,*)'ERROR: 1st line of icein'
        write(*,*)'  must be: ICEMAN-HP or ICEMAN-LP'
        stop
      endif

      if (addph .eq. 'plumehead') then
        write(*,*)'will add the horizontal low velocity anomaly '
        write(*,*)'  to the S-VELOCITY MODEL ONLY'
        plumehead=.TRUE.
      else if (addph .eq. 'null') then
        write(*,*)'will not add plumehead anomaly'
        plumehead=.FALSE.
      else
        write(*,*)'ERROR: 2nd line of icein'
        write(*,*)'  must be: plumehead or null'
        stop
      endif

      if (absanom .eq. 'absolute') then
        write(*,*)'will return absolute velocities (iceman+IASPI)'
        abs=.TRUE.
        pct=.FALSE.
      else if (absanom .eq. 'anomaly') then
        write(*,*)'will return velocity anomalies (km/s) '
        abs=.FALSE.
        pct=.FALSE.
      else if (absanom .eq. 'percent') then
        write(*,*)'will return percent velocity anomalies '
        abs=.FALSE.
        pct=.TRUE.
      else
        write(*,*)'ERROR: 3rd line of icein'
        write(*,*)'  must be: absolute or anomaly'
        stop
      endif
      write(*,*)

c
ccc   loop over all input r, theta, phi
c
      write(*,*)'looping through r, theta, phi in icein'
      open(2,file='iceout')
      write(2,*)'r theta phi Vp Vs'
      write(*,*)'r theta phi Vp Vs'

100   read(1,*,end=200)r,th,ph
      rth = th*torad
      rph = ph*torad
      call get_plume(r,rth,rph,pvs,pvp,prho,abs)

      if (plumehead .AND. r.ge.6161. .AND. r.le.6321.) then
        velicav = 3.97 + (0.0026*(6371.-r))
        dvelicav = velicav-4.516
        pvs=pvs+dvelicav
      endif

      if (pct) then
        realrad=sngl(r)
        l=0
        call mdiaspiice(realrad,l,iaspivp,iaspivs,iaspirho)
        pvp = 100.*pvp/dreal(iaspivp)
        pvs = 100.*pvs/dreal(iaspivs)
      endif

      write(2,*) r,th,ph,pvp,pvs
      write(*,*) r,th,ph,pvp,pvs
      goto 100

200   continue
      close(1)
      close(2)
      write(*,*)'end of icein file'
      write(*,*)
      write(*,*)'output vel is file iceout: r theta phi Vp Vs'
      write(*,*)

      end


c--------------------------------------------------------------

      subroutine read_plume_model(argvp,argvs,argv)

c     reads plume velocity models into common block for later use

      implicit none

      integer MAXPMOD
      parameter (MAXPMOD=100000)
      integer npvp,npvs,i
      character*40 argv,argvp,argvs
      double precision pvp(MAXPMOD),pvs(MAXPMOD),torad
      integer nsq,nxy,nz,ny,nx
      double precision ola,olo,oazi,zmin,zmax,ymin,ymax,xmin,xmax
      double precision dz,dy,dx
      common/model_plume/npvp,pvp,npvs,pvs
      common/grid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &                  zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx

      torad=3.141592654/180.

c     read Vp
      open(unit=10,file=argvp,status='old')
      read(10,*) npvp
      if (npvs.gt.MAXPMOD) then
        stop 'ERROR:read_plume_model: increase MAXPMOD'
      endif
      read(10,*) (pvp(i),i=1,npvp)
      close(10)
      write(*,*)' read: ', argvp

c     read Vs
      write(*,*)' reading: ', argvs
      open(unit=10,file=argvs,status='old')
      read(10,*) npvs
      if (npvs.gt.MAXPMOD) then
        stop 'ERROR:read_plume_model: increase MAXPMOD'
      endif
      read(10,*) (pvs(i),i=1,npvs)
      close(10)
      write(*,*)' read: ', argvs

c     read model cube
      open(unit=10,file=argv,status='old')
      read(10,*) ola,olo,oazi,nz,ny,nx,
     &                  zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
      close(10)
      nxy = nx*ny
      nsq = nxy*nz
      ola = ola*torad
      olo = olo*torad
      oazi = oazi*torad
      write(*,*)' read: ', argv

      if (npvp .lt. nsq) then
        write(*,*) 'ERROR: npvp solution doesnt have enough points!'
      endif
      if (npvs .lt. nsq) then
        write(*,*) 'ERROR: npvs solution doesnt have enough points!'
      endif


      return
      end subroutine read_plume_model

c--------------------------------------------------------------

      subroutine get_plume(radius,theta,phi,ppvs,ppvp,pprho,abs)

ccc   returns velocity anomlaies for given r,theta(lat),phi(lon) point
ccc   prho scales with pvs
ccc   if abs=TRUE then absolute velocity is returned: ICEMAN+iasp91

      implicit none

c     include "constants.h"
      double precision SCALE_ALPHA,SCALE_RHO
      parameter(SCALE_ALPHA=0.55d0,SCALE_RHO=0.4d0)

      double precision radius,theta,phi,ppvs,ppvp,pprho
      double precision x,y,z,pq(100),ptheta,pphi
      integer nq,iq(100),i,l
      logical out,man2surface
      real realrad,iaspivp,iaspivs,iaspirho

      integer MAXPMOD
      parameter (MAXPMOD=100000)
      integer npvp,npvs
      double precision pvp(MAXPMOD),pvs(MAXPMOD)
      logical abs
      common/model_plume/npvp,pvp,npvs,pvs

c     if man2surface=TRUE then vel at 50km is continued to the surface
      man2surface=.FALSE.

      ppvp = 0.
      ppvs = 0.
      pprho = 0.

      call rotate(0,theta,phi,ptheta,pphi)
c     write(*,*)'th,ph,th,ph: ',theta,phi,ptheta,pphi
      call rtfxyz(0,ptheta,pphi,radius,x,y,z)
c     write(*,*)'x,y,z: ',x,y,z
      if (z .lt. 5971.) return
      if (man2surface) then
        if (z.ge.6321.) z=6321.
      endif
      call getc(x,y,z,nq,iq,pq,out)
      if (out) then
c       write(*,*) 'outside box'
        return
      endif
c     write(*,*)'iq: ',(iq(i),i=1,nq)
c     write(*,*)'pq: ',(pq(i),i=1,nq)
      do i=1,nq
        ppvp = ppvp + pq(i)*pvp(iq(i))
        ppvs = ppvs + pq(i)*pvs(iq(i))
      enddo

      pprho=SCALE_RHO*ppvs

      if (abs) then
        realrad=sngl(radius)
        l=0
        call mdiaspiice(realrad,l,iaspivp,iaspivs,iaspirho)
        ppvp = ppvp + dreal(iaspivp)
        ppvs = ppvs + dreal(iaspivs)
        pprho = pprho + dreal(iaspirho)
      endif


c     write(*,*) radius,ppvp,ppvs,pprho

      return
      end subroutine get_plume

c----------------------------------------------------------------

      subroutine rotate(mode,tetain,fiin,tetaout,fiout)

c     Rotates between convensional (on the map) lat lon and the 'new' coords
c     used in the code, ie relative to a defined origin.
c     Origin of 'new' coords (ola olo) passed through common block
c     mode 0, rotate from conventional to new
c     mode 1, rotate from new to conventional

      integer mode
      double precision tetain,fiin,tetaout,fiout

      double precision ang1,ang2,ct,st,cf,sf,ca,sa,co,so,som
      double precision fdiff,sfo,cfo,teller,rnoemer,azimuth

      integer nsq,nxy,nz,ny,nx
      double precision ola,olo,oazi,zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
      common/grid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &			zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx


cinclude 'defines'

      if (mode == 0) then
         ang1 = olo
         ang2 = oazi
      else if (mode == 1) then
         ang2 = olo
         ang1 = oazi
      else 
        tetaout = tetain
        fiout = fiin
        mode = -1
        return
      endif

      ct = cos(tetain)
      st = sin(tetain)
      cf = cos(fiin)
      sf = sin(fiin)
      ca = cos(ola)
      sa = sin(ola)
      co = cos(ang1)
      so = sin(ang1)
      som = ct*ca*(cf*co + sf*so) + st*sa
      tetaout = asin(som) 
      fdiff = fiin - ang1
      sfo = sin(fdiff)
      cfo = cos(fdiff)
      teller = ct*sfo
      rnoemer = st*ca - ct*sa*cfo
      azimuth = atan2(teller,rnoemer)
      fiout = ang2 - azimuth

      return
      end

c----------------------------------------------------------------

      subroutine rtfxyz(mode,teta,fi,r,x,y,z)

c     Code to convert from the rotated ('new') coords ie relative to
c     defined origin, to x,y,z

c     teta = latitude, fi = longitude, relative to new (rotated)
c     northpole and equator. r = distance from Earth's center (0,0,0).
c     if mode=0 r,teta,fi --> x,y,z - if r km, x,y,z will be in km
c     if mode=1 x,y,z --> r,teta,fi

      implicit none

      integer mode
      double precision teta,fi,r,x,y,z
      double precision ct,x2,y2,z2


      if (mode == 0) then
         ct = cos(teta)
         x = r*ct*cos(fi)
         y = r*ct*sin(fi)
         z = r*sin(teta)
      elseif (mode == 1) then
         x2 = x*x
         y2 = y*y
         z2 = z*z
         r = sqrt(x2 + y2 + z2)
         fi = atan2(y,x)
         teta = atan2(z,sqrt(x2 + y2))
      else 
         x=0;y=0;z=0
         r=0;teta=0;fi=0
         mode=-1
      endif

      return
      end subroutine rtfxyz

c----------------------------------------------------------------

      subroutine getc(xp,yp,zp,np,i,f,out)

c     routine to associate grid points around some point

c     passed from above:
c     xp,yp,zp - coords of point of interest
c     obtained here:
c     np - number of grid points associated
c     i(np) - grid indices of these grid points
c     f(np) - weights of grid points, based on distance from point
c       ...total number of points depends on smoothing
c     out - logical - if returned true then poind outside grid

c     xyz grid counts from top layer, node 1 is at xmin,ymax and
c     node nx is at xmax,ymax. node nx+1 is at xmin,ymax-1, etc. till
c     node nx*ny. Next node is at z-1, xmin,ymax, etc.

      implicit none

      integer maxpt
      parameter (maxpt=8)

      integer np,l,i(maxpt)
      double precision f(maxpt),xp,yp,zp,x0,y0,z0
      logical out

      integer nsq,nxy,nz,ny,nx
      double precision ola,olo,oazi,zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
      common/grid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &			zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx

      np = 0
      do l = 1,maxpt 
         i(l) = 0; f(l) = 0. 
      enddo

      x0 = xp
      y0 = yp
      z0 = zp

c     check if x,y,z, is inside square grid: if not return
      out = .false.
      if (x0 .lt. xmin .or. x0 .gt. xmax .or. y0 .lt. ymin .or. 
     &		y0 .gt. ymax .or. z0 .lt. zmin .or. z0 .gt. zmax) then
         out = .true. 
         return
      endif

      call trilin(x0,y0,z0,np,i,f)
c      write(*,*)'i: ',(i(l),l=1,8)
c      write(*,*)'f: ',(f(l),l=1,8)

      return
      end subroutine getc

c-------------------------------------------------------------------

      subroutine trilin(xq,yq,zq,n8,i8,f8)

      implicit none

      integer maxpt
      parameter (maxpt=8)
      integer n8,i8(maxpt),ixmn,ixmx,iymn,iymx,izmn,izmx,lub
      double precision f8(maxpt),xq,yq,zq
      double precision xdist,ydist,zdist,x,y,z,xy,xz,yz,xyz

      integer nsq,nxy,nz,ny,nx
      double precision ola,olo,oazi,zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
      common/grid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &			zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx

      n8 = 8
      
      xdist = (xq - xmin)/dx
      ydist = (ymax - yq)/dy
      zdist = (zmax - zq)/dz

c     coordinates of upper left back corner:
      ixmn = int(xdist) + 1
      ixmx = ixmn + 1
      if (ixmn == nx) then
         ixmx = ixmn; ixmn = ixmx - 1  
      endif
      iymn = int(ydist) + 1
      iymx = iymn + 1
      if (iymn == ny) then
         iymx = iymn; iymn = iymx - 1 
      endif
      izmn = int(zdist) + 1
      izmx = izmn + 1
      if (izmn == nz) then
         izmx = izmn; izmn = izmx - 1 
      endif

      lub = (izmn-1)*nxy + (iymn-1)*nx + ixmn 
      i8(1) = lub
      i8(2) = lub + 1
      i8(3) = lub + nx
      i8(4) = i8(3) + 1
      i8(5) = lub + nxy
      i8(6) = i8(5) + 1
      i8(7) = i8(5) + nx
      i8(8) = i8(7) + 1

c     relative distance of point to 3 sides of cell:
      x = xdist - (ixmn-1)
      y = (iymx-1) - ydist
      z = (izmx-1) - zdist
      
      xy = x*y
      xz = x*z
      yz = y*z
      xyz = xy*z

      f8(1) = yz - xyz
      f8(2) = xyz
      f8(3) = z - yz - xz + xyz
      f8(4) = xz - xyz
      f8(5) = y - xy - yz + xyz
      f8(6) = xy - xyz
      f8(7) = 1 - x - y - z + xy + yz + xz - xyz
      f8(8) = x - xy - xz + xyz

c     # 1 cube:
c     #
c     #        1________2 
c     #        /       /|
c     #      3/______4/ |
c     #       | |_____|_|6 
c     #       |/5     |/ 
c     #      7/_______/8 
c     #
c     #
c     # the index of a node:  i = (iz-1)*nx*ny + (iy-1)*nx + ix

      return
      end

c-------------------------------------------------------------------


      subroutine mdiaspiice(r,l,vp,vs,rho)
c     gives vp and vs at a level of r km in iaspi91 earth model
c     input: r - radius in km
c            l - layer # (top layer is 1, centre 9) l may be 0.
c     output: vp - p velocity - km/s
c             vs - s velocity - km/s
c             l - layer number where r is located (only if l=0 upon input)
c     if l=0 and r exactly at boundary, the lower layer is chosen.
c     if l>0 and r not in the specified layer, vp and vs are extrapolated!
c	rho is a dummy arg right now 

c     RMA May 2000: modified to use iceland bg model for crust
c                   iaspi91 from core to 120km
c                   120 to 30km linear from 4.49 (isasp) to 4.1 (ice)
c                   30km to surface uses ice model:
c                        6371 -  0     2.25
c                        6366 -  5     3.75
c                        6341 - 30     3.85

c     Vp is calculated from Vs for the iceland model:
c     possible range of vpvs is 1.73-1.97 (avg: 1.85)
      vpvs=1.85
c      write(*,*)'mdiaspiice: vpvs: ',vpvs

        x=r/6371.
	rup=r

c	IF(L.GT.0) GO TO (97,95,90,80,70,60,50,40,30,20,10),L	
	if (rup.ge.0..and.rup.lt.1217.1) then 
10	    layer=11
            vp=11.24094-4.09689*x**2
   	    vs=3.56454-3.45241*x**2
	    go to 100
	
	else if (rup.ge.1217.1.and.rup.lt.3482.) then
20   	    layer=10
            vp=10.03904+3.75665*x-13.67046*x**2
            vs=0.0
	    go to 100
         
        else if (rup.ge.3482 .and.rup.lt.3631.) then
30            layer=9
             vp=14.49470-1.47089*x
             vs=8.16616-1.58206*x
	     go to 100
        
        else if (rup.ge.3631. .and. rup.lt.5611.) then
40           layer=8
            vp=25.1486-41.1538*x+51.9932*x**2-26.6083*x**3
            vs=12.9303-21.2590*x+27.8988*x**2-14.1080*x**3
	    go to 100
        
        else if (rup.ge.5611. .and. rup.lt.5711.) then
50            layer=7
            vp=25.96984-16.93412*x
            vs=20.76890-16.53147*x
	    go to 100
        
        else if (rup.ge.5711. .and. rup.lt.5961.) then
60            layer=6
            vp=29.38896-21.40656*x
            vs=17.70732-13.50652*x
	    go to 100
        
        else if (rup.ge.5961..and.rup.lt.6161.) then
70            layer=5
            vp=30.78765-23.25415*x
            vs=15.24213-11.08552*x
	    go to 100
        
        else if (rup.ge.6161. .and. rup.lt.6251.) then
80            layer=4
            vp=25.41389-17.69722*x
            vs=5.75020-1.27420*x
	    go to 100
        
c        else if (rup.ge.6251. .and. rup.lt.6336.) then
c90            layer=3
c            vp=8.78541-0.74953*x
c            vs=6.706231-2.248585*x
c	    go to 100
c        
c        else if (rup.ge.6336. .and. rup.lt.6351.) then 
c95            layer=2
c            vp=6.50
c            vs=3.75
c	    go to 100
c        
c        else if (rup.ge.6351. .and. rup.le.6371.) then 
c97            layer=1
c            vp=5.80
c            vs=3.36
c	go to 100
c        endif 

c       this is the iceland model....
        else if (rup.ge.6251. .and. rup.lt.6341.) then
             vs = 4.1 + (6341.-rup)*0.39/90.
             vp = vs * vpvs
             goto 100
        else if (rup.ge.6341. .and. rup.lt.6366.) then
             vs = 3.75 + (6366.-rup)*0.10/25.
             vp = vs * vpvs
             goto 100
        else if (rup.ge.6366. .and. rup.le.6371.) then
             vs = 2.25 + (6371.-rup)*1.50/5.
             vp = vs * vpvs
             goto 100
        endif 

100	if (l.eq.0) then
		l=layer
	end if
        return
        end