program extracticecrt

c     Code returns the S-velocity structure for the Icelandic crust 
c     given r (radius), theta (lat), and phi (lon).
c     The velocities are extracted from the ICECRTb S-vel model. 
c     The code also returns a P-velocity which is determined from the S-vel
c     using a Vp/Vs ratio defined: Vp/Vs = 1.78+0.004.depth.
c     For more info on the development and resolution of the S-vel model, and
c     on the origin of the Vp/Vs ratio  see Allen et al, Plume driven plumbing 
c     and crustal formation in Iceland, 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     extracticecrt.f (this file) sol.ICECRTb.abs model.cube.ICECRTb.abs
c
c     Code opens the file "icecrtin" containing a list of r, theta, phi:
c     1th line: r, theta, phi
c     2th line: r, theta, phi
c     3th line: r, theta, phi
c     etc

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

c     EXAMPLE - input file: icecrtin
c         6371.00 64.70 -20.00
c         6366.00 64.70 -20.00
c         6361.00 64.70 -20.00
c         6356.00 64.70 -20.00
c         6351.00 64.70 -20.00
c         6346.00 64.70 -20.00
c         6341.00 64.70 -20.00
c         6336.00 64.70 -20.00
c         6331.00 64.70 -20.00
c
c     EXAMPLE - output file: icecrtout
c         r theta phi Vp Vs
c          6371.  64.7 -20.  3.9540639  2.22138425
c          6366.  64.7 -20.  5.82956455  3.23864702
c          6361.  64.7 -20.  6.48264273  3.56189166
c          6356.  64.7 -20.  6.82676614  3.71019904
c          6351.  64.7 -20.  7.08961967  3.81162353
c          6346.  64.7 -20.  7.23998113  3.85105384
c          6341.  64.7 -20.  7.61496239  4.00787499
c          6336.  64.7 -20.  7.90348531  4.11639865
c          6331.  64.7 -20.  8.03936154  4.14400084

      implicit none
      character*40 mc,ms
      double precision torad,r,th,ph,rth,rph,pvs,pvp,prho
      logical infile

      write(*,*)
      write(*,*)'EXTRACTICECRT: '
      write(*,*)'  will extract Vs from the ICECRTb vel model'
      write(*,*)'  and return Vp using Vs and Vp/Vs = 1.78+0.004.depth'
      write(*,*)

ccc   initialize
      mc='model.cube.ICECRTb.abs'
      ms='sol.ICECRTb.abs'
      torad=3.141592654/180.

c
ccc   check all necessary files are present
c
      infile=.FALSE.
      inquire(file='icecrtin',exist=infile)
      if (infile) then
        write(*,*)'found: icecrtin'
      else
        write(*,*)'ERROR: can not find: icecrtin'
        stop
      endif
      infile=.FALSE.
      inquire(file=ms,exist=infile)
      if (infile) then
        write(*,*)'found: ',ms
      else
        write(*,*)'ERROR: can not find: ',ms
        stop
      endif
      infile=.FALSE.
      inquire(file=mc,exist=infile)
      if (infile) then
        write(*,*)'found: ',mc
      else
        write(*,*)'ERROR: can not find: ',mc
        stop
      endif
      write(*,*)'found all necessary input files, continuing...'
      write(*,*)

c
ccc   read in the velocity model
c
      call read_icecrust_model(ms,mc)
      write(*,*)'read in the model solution, continuing...'
      write(*,*)

c
ccc   open input file
c
      open(unit=1, file='icecrtin',status='old')

c
ccc   loop over all input r, theta, phi
c
      write(*,*)'looping through r, theta, phi in icecrtin'
      open(2,file='icecrtout')
      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_icecrust(r,rth,rph,pvs,pvp,prho)

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

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

      end





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

      subroutine read_icecrust_model(argv,argc)

c     reads Iceland crustal velocity models into common block for later use

      implicit none

      integer MAXCMOD
      parameter (MAXCMOD=1000000)
      integer ncvs,i
      character*40 argv,argc
      double precision cvs(MAXCMOD),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_icecrust/ncvs,cvs
      common/crtgrid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &			zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx

      torad=3.141592654/180.

c     read Vs
c     argv='sol.ICECRTb.abs'
      open(unit=10,file=argv,status='old')
      read(10,*) ncvs
      if (ncvs.gt.MAXCMOD) then
        stop 'ERROR:read_icecrust_model: increase MAXCMOD'
      endif
      read(10,*) (cvs(i),i=1,ncvs)
      close(10)

c     read model cube
c     argc='model.cube.ICECRTb.abs'
      open(unit=10,file=argc,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

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

      return
      end subroutine read_icecrust_model

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

      subroutine get_icecrust(radius,theta,phi,vs,vp,rho)

c     given r,theta(co-lat),phi, code obtains corresponding ICECRTb.abs Vs
c     if it is a crustal velocity then it is returned, otherwise 0. 
c     is returned for Vs.  NB sub returns ABSLOUTE Vs, Vp and rho.
c     Vp is obtained using Vp/Vs = 1.78+0.004.depth
c     rho is obtained from Vp using the relations from Darbyshire 2000
c     This is intended to overlay the crustal velocities over the ICEMAN
c     models

      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,vs,vp,rho
      double precision x,y,z,pq(100),ptheta,pphi,depth
      double precision prop
      real realrad,iaspivp,iaspivs,iaspirho
      integer nq,iq(100),i,l
      logical out,smooth

      integer MAXCMOD
      parameter (MAXCMOD=1000000)
      integer ncvs
      double precision cvs(MAXCMOD)
      common/model_icecrust/ncvs,cvs


      vp = 0.
      vs = 0.
      rho = 0.
      smooth=.FALSE.

      call rotatecrt(0,theta,phi,ptheta,pphi)
c     write(*,*)'th,ph,th,ph: ',theta,phi,ptheta,pphi
      call rtfxyzcrt(0,ptheta,pphi,radius,x,y,z)

c
ccc   get velocity
      call getccrt(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
        vs = vs + pq(i)*cvs(iq(i))
      enddo

c     write(*,*)'vs: ',vs

ccc     estimate vp 
        depth = 6371.-radius
        vp = vs*(1.78+(0.004*depth))
ccc     estiate rho (after Darbyshire 2000)
        if (vp.le.4.5) then
          rho = 1000*(-0.6997+(2.2302*vp)-(0.598*(vp**2.))
     &              +(0.07036*(vp**3.))-(0.0028311*(vp**4.)))
        else if (vp.gt.4.5 .AND. vp.le.6.6) then
          rho = 1000*(3.81-(6.0/vp))
        else if (vp.gt.6.6) then
          rho = 1000*(5.23-(15.38/vp))
        endif

        
c     write(*,*) radius,vp,vs,rho

      return
      end subroutine get_icecrust

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

      subroutine rotatecrt(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
      double precision dz,dy,dx
      common/crtgrid/ola,olo,oazi,nsq,nxy,nz,ny,nx,
     &			zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx


      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 rtfxyzcrt(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 rtfxyzcrt

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

      subroutine getccrt(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
      double precision dz,dy,dx
      common/crtgrid/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 trilincrt(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 getccrt

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

      subroutine trilincrt(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
      double precision dz,dy,dx
      common/crtgrid/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