      program parameter_at_point
C-----------------------------------------------------------------
C     version 1.2
C     
C     Supporting information for the Geophysical Research Letters paper (2015):
C     
C     Inverting Interpolated Receiver Functions with Surface-Wave Dispersion and Gravity:
C          Application to the western US and adjacent Canada and Mexico
C        Chengping Chai, Charles J. Ammon, Monica Maceira, Robert B. Herrmann.
C      
C-----------------------------------------------------------------
C     Simple program to interpolate velocity/density depth profile at a given 
C     latitude, longitude (within the lat/lon range 25˚ to 55˚ N  and 
C     97˚ to 127˚ W ) by interpolating thev elocities between the node
C     points of the inversion. The model files can be found in the supplements of the
C     above paper.
C     Model files:
C        2015GL063733-ds01.txt (Initial model)
C        2015GL063733-ds02.txt (Inverted model) 
C
C     Rename the file to parameter_at_point.f, then compile the code with
C         gfortran parameter_at_point.f -o parameter_at_point
C     to run the program, use following syntax in the terminal
C         ./parameter_at_point [-h] -f FILENAME -t PARAMETER_TYPE -lat LAT -lon
C                             LON [-dep DEPTH]
C     Optional arguments:
C       -h, --help         show this help message and exit
C       -f FILENAME        Filename of the model
C       -t PARAMETER_TYPE  vs, vp, rho or all
C       -lat LAT           Latitude in decimal degrees, +N.
C       -lon LON           Longitude in decimal degrees, +E.
C       -dep DEPTH         Depth in km.
C     Filename of the model file (-f) is required.
C     
C     For example:
C     
C         ./parameter_at_point -f 2015GL063733-ds01.txt -t vp -lat 35.5 -lon -115
C     
C     Last update:
C              Chengping Chai 2015/01/25
C-----------------------------------------------------------------
      integer NLMAX,NPMAX,NTMAX,NWFMAX,NCUMAX,NGRAV
      integer NMAX,MAXCELLSX,MAXCELLSY,MODELSIZE,NMAXALLCELLS
      PARAMETER(NLMAX=100, NPMAX=150, NTMAX=2048, NWFMAX=6)
      PARAMETER(NCUMAX=4, NGRAV=1)
      PARAMETER(NMAX=NPMAX*NCUMAX+NTMAX*NWFMAX+2*NLMAX + NGRAV)
      PARAMETER(MAXCELLSX=32, MAXCELLSY=32)
      PARAMETER(MODELSIZE=MAXCELLSX*MAXCELLSY*NLMAX)   
      PARAMETER(NMAXALLCELLS = MAXCELLSX*MAXCELLSY*(1+NMAX))   

      real latmax,lonmin,dlat,dlon,clat,clon,celev
      real vs3d(MODELSIZE),rho3d(MODELSIZE),lyrthick3d(MODELSIZE)
      real vp3d(MODELSIZE),vsPrior(MODELSIZE),vs3d0(MODELSIZE)
      real smthwt3d(MODELSIZE),apmwt3d(MODELSIZE)
      real lyrtopdepth(MODELSIZE)
      real ap3d0(MODELSIZE), ae3d0(MODELSIZE),rho3d0(MODELSIZE)
      integer nx3d,ny3d,nz3d,nlyr3d,irmin,irmax,icmin,icmax
      character*2 geocodes(MAXCELLSX*MAXCELLSY)
C
      common/Model3D/nx3d,ny3d,latmax,lonmin,dlat,dlon,
     1    vs3d,rho3d,nlyr3d,lyrthick3d,vp3d,lyrtopdepth,
     2    vs3d0,vsPrior,apmwt3d,smthwt3d,ae3d0,ap3d0,
     3    geocodes,
     4    irmin,irmax,icmin,icmax

      character*256 fname
      integer ilenfname
      integer iarg
      integer NARG
      parameter (NARG=12)
      character*32 arg, arglist(NARG)
      character*4 parameter_type
      real point_data(NLMAX), vp11(NLMAX), vs11(NLMAX), rho11(NLMAX)
      real tops(NLMAX)
      logical flag
      integer ipalen
C-----------------------------------------------------------------
C     get argument list
      do i = 1, iargc()
        call getarg(i, arg)
        arglist(i) = arg
      end do
C      
      flag = .TRUE.
      if (arglist(1)(1:2) .EQ. '-h') then
        write(*,*) 'usage: parameter_at_point.py [-h] -f FILENAME -t',
     1   ' PARAMETER_TYPE -lat LAT -lon'
        write(*,*) 'LON [-dep DEPTH]'
        write(*,*) ''
        write(*,*) 'Interpolate velocity/density from the',
     1   'initial/final model.'
        write(*,*) ''
        write(*,*) 'optional arguments:'
        write(*,*) '  -h, --help         show this help message',
     1   ' and exit'
        write(*,*) '  -f FILENAME        Filename of the model'
        write(*,*) '  -t PARAMETER_TYPE  vs, vp, rho or all'
        write(*,*) '  -lat LAT           Latitude in decimal',
     1   ' degrees, +N.'
        write(*,*) '  -lon LON           Longitude in decimal',
     1  ' degrees, +E.'
        write(*,*) '  -dep DEPTH         Depth in km.'
      end if
C-----------------------------------------------------------------
C     set up the 3D grid
      nx3d = 30
      ny3d = 30
      nlyr3d = 99
      nz3d = nlyr3d
      latmax = 54.5
      lonmin = -126.5
      dlat = 1.0
      dlon = 1.0
      point_lat = 54.5
      point_lon = -124.5
      point_data = 0.
C-----------------------------------------------------------------
      if (arglist(1)(1:2) .EQ. '-f') then
        fname = arglist(2)
      else if (arglist(2)(1:2) .EQ. '-f') then
        fname = arglist(3)
      else
        flag = .FALSE.
        write(*,*) 'Filename of the model file (-f) is required.'
        stop
      end if
C-----------------------------------------------------------------
      if (flag) then
        if (arglist(3)(1:2) .EQ. '-t') then
          ipalen = lnblnk(arglist(4))
          parameter_type = arglist(4)(1:ipalen)
        else if (arglist(4)(1:2) .EQ. '-t') then
          ipalen = lnblnk(arglist(5))
          parameter_type = arglist(5)(1:ipalen)
        else
          flag = .FALSE.
          write(*,*) 'Parameter type (-t) is required.'
          stop
        end if
      end if
C-----------------------------------------------------------------
      if (flag) then
        if (arglist(5)(1:4) .EQ. '-lat') then
          read(arglist(6),*) point_lat
        else if (arglist(6)(1:4) .EQ. '-lat') then
          read(arglist(7),*) point_lat
        else
          flag = .FALSE.
          write(*,*) 'Latitude (-lat) is required.'
          stop
        end if
        if (arglist(7)(1:4) .EQ. '-lon') then
          read(arglist(8),*) point_lon
        else if (arglist(8)(1:4) .EQ. '-lon') then
          read(arglist(9),*) point_lon
        else
          flag = .FALSE.
          write(*,*) 'Longitude (-lon) is required.'
          stop
        end if
      end if
C-----------------------------------------------------------------
      if (flag) then
        if (arglist(9)(1:4) .EQ. '-dep') then
          read(arglist(10),*) point_depth
        else if (arglist(10)(1:4) .EQ. '-dep') then
          read(arglist(11),*) point_depth
        else 
          point_depth = -10000.0
        end if
      end if
C-----------------------------------------------------------------
      if (flag) then
        ilenfname = lnblnk(fname)
        call Read3DModelToTextfile(fname(1:ilenfname))
        tops(1:nlyr3d) = lyrtopdepth(1:nlyr3d)
        if (abs(point_depth + 10000.0) .LT. 1e-6) then
          if (parameter_type(1:2) .EQ. 'vs') then
            call vel_at_latlon(point_lat, point_lon, vs3d, vs11)
            write(*,*) 'layer top(km)  Vs(km/s)' 
            do ilay=1,nlyr3d
              write(*,'(2(f10.4,2x))') tops(ilay),vs11(ilay)
            end do
          else if (parameter_type(1:2) .EQ. 'vp') then
            call vel_at_latlon(point_lat, point_lon, vp3d, vp11)
            write(*,*) 'layer top(km)  Vp(km/s)' 
            do ilay=1,nlyr3d
              write(*,'(2(f10.4,2x))') tops(ilay),vp11(ilay)
            end do
          else if (parameter_type(1:3) .EQ. 'rho') then
            call vel_at_latlon(point_lat, point_lon, rho3d, rho11)
            write(*,*) 'layer top(km)  Rho(g/cm3)' 
            do ilay=1,nlyr3d
              write(*,'(2(f10.4,2x))') tops(ilay),rho11(ilay)
            end do
          else if (parameter_type(1:3) .EQ. 'all') then
            call vel_at_latlon(point_lat, point_lon, vp3d, vp11)
            call vel_at_latlon(point_lat, point_lon, vs3d, vs11)
            call vel_at_latlon(point_lat, point_lon, rho3d, rho11)
            write(*,*) 'layer top(km)  vp(km/s)    vs(km/s)   ',
     1    ' rho(km/s)'
            do ilay=1,nlyr3d
              write(*,'(4(f10.4,2x))') tops(ilay),vp11(ilay),
     1                             vs11(ilay),rho11(ilay)
            end do
          else
            write(*,*) 'Parameter type error, choose from',
     1    ' vs, vp, rho, all.'
          end if
        else
          if (parameter_type(1:2) .EQ. 'vs') then
            call vel_at_latlon(point_lat, point_lon, vs3d, vs11)
            call vel_at_depth(point_depth, tops, vs11, vs11_dep)
            write(*,*) 'layer top(km)  Vs(km/s)' 
            write(*,'(2(f10.4,2x))') point_depth, vs11_dep
          else if (parameter_type(1:2) .EQ. 'vp') then
            call vel_at_latlon(point_lat, point_lon, vp3d, vp11)
            call vel_at_depth(point_depth, tops, vp11, vp11_dep)
            write(*,*) 'layer top(km)  Vp(km/s)' 
            write(*,'(2(f10.4,2x))') point_depth, vp11_dep
          else if (parameter_type(1:3) .EQ. 'rho') then
            call vel_at_latlon(point_lat, point_lon, rho3d, rho11)
            call vel_at_depth(point_depth, tops, rho11, rho11_dep)
            write(*,*) 'layer top(km)  Rho(g/cm3)' 
            write(*,'(2(f10.4,2x))') point_depth, rho11_dep
          else if (parameter_type(1:3) .EQ. 'all') then
            call vel_at_latlon(point_lat, point_lon, vp3d, vp11)
            call vel_at_depth(point_depth, tops, vp11, vp11_dep)
            call vel_at_latlon(point_lat, point_lon, vs3d, vs11)
            call vel_at_depth(point_depth, tops, vs11, vs11_dep)
            call vel_at_latlon(point_lat, point_lon, rho3d, rho11)
            call vel_at_depth(point_depth, tops, rho11, rho11_dep)
            write(*,*) 'layer top(km)  vp(km/s)    vs(km/s)   ',
     1    ' rho(km/s)'
            write(*,'(4(f10.4,2x))') point_depth, vp11_dep,
     1                             vs11_dep,rho11_dep
          else
            write(*,*) 'Parameter type error, choose from',
     1    ' vs, vp, rho, all.'
          end if
        end if
      end if
      end program parameter_at_point
C ================================================================
      subroutine Read3DModelToTextFile(fname)
      character(*) fname
      character*128 tempstr,tempstr1
            integer NLMAX,NPMAX,NTMAX,NWFMAX,NCUMAX,NGRAV
      integer NMAX,MAXCELLSX,MAXCELLSY,MODELSIZE,NMAXALLCELLS
      PARAMETER(NLMAX=100, NPMAX=150, NTMAX=2048, NWFMAX=6)
      PARAMETER(NCUMAX=4, NGRAV=1)
      PARAMETER(NMAX=NPMAX*NCUMAX+NTMAX*NWFMAX+2*NLMAX + NGRAV)
      PARAMETER(MAXCELLSX=32, MAXCELLSY=32)
      PARAMETER(MODELSIZE=MAXCELLSX*MAXCELLSY*NLMAX)   
      PARAMETER(NMAXALLCELLS = MAXCELLSX*MAXCELLSY*(1+NMAX))   

      real latmax,lonmin,dlat,dlon,clat,clon,celev
      real vs3d(MODELSIZE),rho3d(MODELSIZE),lyrthick3d(MODELSIZE)
      real vp3d(MODELSIZE),vsPrior(MODELSIZE),vs3d0(MODELSIZE)
      real smthwt3d(MODELSIZE),apmwt3d(MODELSIZE)
      real lyrtopdepth(MODELSIZE)
      real ap3d0(MODELSIZE), ae3d0(MODELSIZE),rho3d0(MODELSIZE)
      integer nx3d,ny3d,nz3d,nlyr3d,irmin,irmax,icmin,icmax
      character*2 geocodes(MAXCELLSX*MAXCELLSY)
C
      common/Model3D/nx3d,ny3d,latmax,lonmin,dlat,dlon,
     1    vs3d,rho3d,nlyr3d,lyrthick3d,vp3d,lyrtopdepth,
     2    vs3d0,vsPrior,apmwt3d,smthwt3d,ae3d0,ap3d0,
     3    geocodes,
     4    irmin,irmax,icmin,icmax

      real temp
      
C     
C     Save the 3d model in one large text file
C     
      open (2,file=fname,status='old')
C     
      ncells = nx3d*ny3d
      do iycell=1, ny3d
         do ixcell=1, nx3d
            icell = (iycell-1)*nx3d + ixcell
            read(2,701) tempstr,irdcell,tempstr,icrow,tempstr,
     1                iccol,tempstr,clat,tempstr,clon,tempstr,
     2                geocodes(icell)
 701        format(a23,i3.3,a10,i3.3,a1,i3.3,
     1           a6,f7.3,a6,f8.3,a7,a2)

            icellModPntrOffset = ((iycell-1)*nx3d+ixcell-1)*nlyr3d
            read (2,*) nlyr3d
            do ilay=1,nlyr3d
               i3d = ilay + icellModPntrOffset
               read(2,*) irdlay, vp3d(i3d),vs3d(i3d), rho3d(i3d), 
     1              lyrthick3d(i3d), smthwt3d(i3d), vs3d0(i3d), 
     2              apmwt3d(i3d),ap3d0(i3d), ae3d0(i3d), 
     3              lyrtopdepth(i3d)
            end do
C     
         end do
      end do
      close (2)

C     end loop over cells

      return
      end subroutine Read3DModelToTextfile
C-----------------------------------------------------------------
      subroutine vel_at_latlon(point_lat,point_lon,prd_data,point_data)
            integer NLMAX,NPMAX,NTMAX,NWFMAX,NCUMAX,NGRAV
      integer NMAX,MAXCELLSX,MAXCELLSY,MODELSIZE,NMAXALLCELLS
      PARAMETER(NLMAX=100, NPMAX=150, NTMAX=2048, NWFMAX=6)
      PARAMETER(NCUMAX=4, NGRAV=1)
      PARAMETER(NMAX=NPMAX*NCUMAX+NTMAX*NWFMAX+2*NLMAX + NGRAV)
      PARAMETER(MAXCELLSX=32, MAXCELLSY=32)
      PARAMETER(MODELSIZE=MAXCELLSX*MAXCELLSY*NLMAX)   
      PARAMETER(NMAXALLCELLS = MAXCELLSX*MAXCELLSY*(1+NMAX))   

      real latmax,lonmin,dlat,dlon,clat,clon,celev
      real vs3d(MODELSIZE),rho3d(MODELSIZE),lyrthick3d(MODELSIZE)
      real vp3d(MODELSIZE),vsPrior(MODELSIZE),vs3d0(MODELSIZE)
      real smthwt3d(MODELSIZE),apmwt3d(MODELSIZE)
      real lyrtopdepth(MODELSIZE)
      real ap3d0(MODELSIZE), ae3d0(MODELSIZE),rho3d0(MODELSIZE)
      integer nx3d,ny3d,nz3d,nlyr3d,irmin,irmax,icmin,icmax
      character*2 geocodes(MAXCELLSX*MAXCELLSY)
!
      common/Model3D/nx3d,ny3d,latmax,lonmin,dlat,dlon,
     1    vs3d,rho3d,nlyr3d,lyrthick3d,vp3d,lyrtopdepth,
     2    vs3d0,vsPrior,apmwt3d,smthwt3d,ae3d0,ap3d0,
     3    geocodes,
     4    irmin,irmax,icmin,icmax

      logical flag
      real point_data(NLMAX)
      real data0(NLMAX), data2(NLMAX), data11(NLMAX)
      real data00(NLMAX), data02(NLMAX), data20(NLMAX), data22(NLMAX)
      real data01(NLMAX), data21(NLMAX)
      real prd_data(MODELSIZE)
      real lat0, lat2, lon0, lon2
      nlat = ny3d
      nlon = nx3d
      flag = .TRUE.
      data0 = 0.
      data2 = 0.
      data11 = 0.
      data00 = 0.
      data01 = 0.
      data02 = 0.
      data20 = 0.
      data21 = 0.
      data22 = 0.
      if ((point_lat .GT. latmax) .OR. 
     1          (point_lat .LT. latmax-nlat*dlat)) then
          write(*,*) 'Latitude out of range!'
          flag = .FALSE.
          stop
      end if
      if ((point_lon .LT. lonmin) .OR. 
     1          (point_lon .GT. lonmin+nlon*dlon)) then
          write(*,*) 'Longitude out of range!'
          flag = .FALSE.
          stop
      end if
      if (flag) then
          ilat0 = int((latmax-point_lat)/dlat)
          ilat2 = ilat0 + 1
          lat0 = latmax - real(ilat0) * dlat
          lat2 = latmax - real(ilat2) * dlat
          ilon0 = int((point_lon-lonmin)/dlon)
          ilon2 = ilon0 + 1
          lon0 = lonmin + real(ilon0) * dlon
          lon2 = lonmin + real(ilon2) * dlon
          if (abs(lat0-point_lat) .LT. 1e-7) then
            do ilay=1, nlyr3d
              icell = ilat0*nlon*nlyr3d + ilon0*nlyr3d + ilay
              data0(ilay) = prd_data(icell)
              icell = ilat0*nlon*nlyr3d + ilon2*nlyr3d + ilay
              data2(ilay) = prd_data(icell)
              data11(ilay) = (lon2-point_lon)/(lon2-lon0)*data0(ilay)+
     1                       (point_lon-lon0)/(lon2-lon0)*data2(ilay)
            end do
          else if (abs(lon0-point_lon) .LT. 1e-7) then
            do ilay=1, nlyr3d
              icell = ilat0*nlon*nlyr3d + ilon0*nlyr3d + ilay
              data0(ilay) = prd_data(icell)
              icell = ilat2*nlon*nlyr3d + ilon0*nlyr3d + ilay
              data2(ilay) = prd_data(icell)
              data11(ilay) = (point_lat-lat0)/(lat2-lat0)*data0(ilay)+
     1                       (lat2-point_lat)/(lat2-lat0)*data2(ilay)
            end do
          else
            do ilay=1, nlyr3d
              icell = ilat0*nlon*nlyr3d + ilon0*nlyr3d + ilay
              data00(ilay) = prd_data(icell)
              icell = ilat0*nlon*nlyr3d + ilon2*nlyr3d + ilay
              data02(ilay) = prd_data(icell)
              icell = ilat2*nlon*nlyr3d + ilon0*nlyr3d + ilay
              data20(ilay) = prd_data(icell)
              icell = ilat2*nlon*nlyr3d + ilon2*nlyr3d + ilay
              data22(ilay) = prd_data(icell)
              data01(ilay) = (lon2-point_lon)/(lon2-lon0)*data00(ilay)+
     1                       (point_lon-lon0)/(lon2-lon0)*data02(ilay)
              data21(ilay) = (lon2-point_lon)/(lon2-lon0)*data20(ilay)+
     1                       (point_lon-lon0)/(lon2-lon0)*data22(ilay)
              data11(ilay) = (point_lat-lat0)/(lat2-lat0)*data01(ilay)+
     1                       (lat2-point_lat)/(lat2-lat0)*data21(ilay)
            end do
          end if
          do ilay=1, nlyr3d
            point_data(ilay) = data11(ilay)
          end do
      end if
      end subroutine vel_at_latlon
C-----------------------------------------------------------------
      subroutine vel_at_depth(depth, tops, profile, depth_data)
            integer NLMAX,NPMAX,NTMAX,NWFMAX,NCUMAX,NGRAV
      integer NMAX,MAXCELLSX,MAXCELLSY,MODELSIZE,NMAXALLCELLS
      PARAMETER(NLMAX=100, NPMAX=150, NTMAX=2048, NWFMAX=6)
      PARAMETER(NCUMAX=4, NGRAV=1)
      PARAMETER(NMAX=NPMAX*NCUMAX+NTMAX*NWFMAX+2*NLMAX + NGRAV)
      PARAMETER(MAXCELLSX=32, MAXCELLSY=32)
      PARAMETER(MODELSIZE=MAXCELLSX*MAXCELLSY*NLMAX)   
      PARAMETER(NMAXALLCELLS = MAXCELLSX*MAXCELLSY*(1+NMAX))   

      real latmax,lonmin,dlat,dlon,clat,clon,celev
      real vs3d(MODELSIZE),rho3d(MODELSIZE),lyrthick3d(MODELSIZE)
      real vp3d(MODELSIZE),vsPrior(MODELSIZE),vs3d0(MODELSIZE)
      real smthwt3d(MODELSIZE),apmwt3d(MODELSIZE)
      real lyrtopdepth(MODELSIZE)
      real ap3d0(MODELSIZE), ae3d0(MODELSIZE),rho3d0(MODELSIZE)
      integer nx3d,ny3d,nz3d,nlyr3d,irmin,irmax,icmin,icmax
      character*2 geocodes(MAXCELLSX*MAXCELLSY)
C
      common/Model3D/nx3d,ny3d,latmax,lonmin,dlat,dlon,
     1    vs3d,rho3d,nlyr3d,lyrthick3d,vp3d,lyrtopdepth,
     2    vs3d0,vsPrior,apmwt3d,smthwt3d,ae3d0,ap3d0,
     3    geocodes,
     4    irmin,irmax,icmin,icmax

      integer depth_index,idep0, idep2
      real depth0, depth2
      real tops(NLMAX)
      real profile(NLMAX)
      depth_index = 0
      do ilay=1, nlyr3d
        if (abs(depth-tops(ilay)) .LT. 1e-6) then
          depth_index = ilay
        else
          depth_index = 0
        end if
      end do
      if (depth_index .NE. 0) then
        data1 = profile(depth_index)
      else if (depth > tops(1) .AND. depth < tops(nlyr3d)) then
        i = 1
        depth0 = tops(i)
        depth2 = tops(i+1)
        do while ((depth0-depth)*(depth2-depth) .GT. 0)
          i = i + 1
          depth0 = tops(i)
          depth2 = tops(i+1)
        end do
        do ilay=1, nlyr3d
          if (abs(depth0-tops(ilay)) .LT. 1e-6) then
            idep0 = ilay
          else if (abs(depth2-tops(ilay)) .LT. 1e-6) then
            idep2 = ilay
          end if
        end do
        data0 = profile(idep0)
        data2 = profile(idep2)
        data1 = (depth2-depth)/(depth2-depth0)*data0 + 
     1          (depth-depth0)/(depth2-depth0)*data2
      else
        write(*,*) 'Depth input error!'
        stop
      end if
      depth_data = data1
      end subroutine vel_at_depth