c23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
c  This program takes a velocity file written by simulps9 and
c  computes a new velocity file on a new grid by interpolating
c  from the old file
c
c  Donna Eberhart-Phillips 09-feb-87
c
c  15-feb-1989  now allow different origins for old and new grids
c
      common/vmod3d/iuses,nx,ny,nz,xn(30),yn(30),zn(30),vel(30,30,40)
      common/vmodnw/ mx,my,mz,xm(30),ym(30),zm(30),vnew(30,30,40)
        common/arrays/iastns,iaeqs,ianod,ianodi,iandst,iasltn,iaobs
      common/step/istep
      character*10 day,tm
      parameter(zero=0.0,izero=0)
      ianod=12700                         !was 3200
c  input old origin of coordinates
      open(unit=02,status='old',form='formatted',read only)
      rewind (02)
      call setor1(2)
      close(02)
c  input new origin of coordinates
      open(unit=12,status='old',form='formatted',read only)
      rewind (12)
      call setor2(12)
      close (12)
c  input medium model
      open(unit=03,status='old',form='formatted',read only)
      rewind (03)
      open(unit=13,status='old',form='formatted',read only)
      rewind (13)
c  new grid
      call input3
      close(03)
      close(13)
c
c  Compute velocities on new grid by interpolation
c
      do 120 kv=1,iuses
      isp=kv-1
      do 100 k=1,mz
      kp=k+(kv-1)*mz
        do 90 j=1,my
          do 80 i=1,mx
          call latlon(xm(i),ym(j),lat,xlat,lon,xlon)
          dlat=float(lat)
          dlon=float(lon)
          call disto(dlat,xlat,dlon,xlon,x,y)
      write(16,1616) xm(i),ym(j),lat,xlat,lon,xlon
      write(16,1617) x,y,dlat,xlat,dlon,xlon
 1616 format(' latlon',2f8.2,i4,f6.2,i5,f6.2)
 1617 format(' disto ',2f8.2,f4.0,f6.2,f5.0,f6.2)
          call vel3(isp,x,y,zm(k),v)
          vnew(i,j,kp)=v
      if(isp.eq.1) then                                        !eh added 6 March 1994
			vnew(i,j,kp)=vnew(i,j,k)/v                           
            endif                                              !eh added end 
   80     continue
   90   continue
  100 continue
  120 continue
c
c   Output new velocity model to file 23
c
c  Write the final velocity model to file23 in a format that
c  can be used as input for another run.                     
      open(unit=23,carriagecontrol='list',status='new')
      call date(day)
      call time(tm)
      write(23,2302) istep,mx,my,mz,iuses,day,tm
 2302 format(f4.1,4i3,10x,'computed',2x,a9,1x,a8)
      write(23,2304) (xm(i),i=1,mx)
      write(23,2304) (ym(i),i=1,my)
      write(23,2304) (zm(i),i=1,mz)
 2304 format(30f6.1)
      write(23,2305) izero,izero,izero
 2305 format(3i3)
      do 38 kv=1,iuses
      do 36 k=1,mz
        k2=k+(kv-1)*mz
        write(23,2311) (vnew(i,1,k2),i=1,mx)
        do 34 j=2,my-1
        write(23,2311) (vnew(i,j,k2),i=1,mx)
   34 continue
        write(23,2311) (vnew(i,my,k2),i=1,mx)
   36 continue
   38 continue
c
 2311 format(30f5.2)
   50 close(unit=23)
c
      stop
c***** end of main program *****
      end
c
c
      subroutine bldmap
      integer ixloc,iyloc,izloc
      common/vmod3d/iuses,nx,ny,nz,xn(30),yn(30),zn(30),vel(30,30,40)
      common/locate/ xl,yl,zl,ixloc(6500),iyloc(6500),izloc(6500)
      common/step/istep
c     array size limits
      ixkms=6500
      iykms=6500
      izkms=6500
c
c     write(6,400)
c 400 format(' subroutine bldmap')
      xl=1.-xn(1)
      ixmax=xn(nx)+xl
      yl=1.-yn(1)
      iymax=yn(ny)+yl
      zl=1.-zn(1)
      izmax=zn(nz)+zl
c     write(6,402)ixmax,iymax,izmax
c 402 format(' array sizes: ',3i5)
c
c  Check for array size overflow
      if(ixmax.gt.ixkms.or.iymax.gt.iykms.or.izmax.gt.izkms)goto 330
      ix=1
      do 10 i=1,ixmax
c
      ix1=ix+1
c
      xnow=float(i)-xl
      if (xnow.ge.xn(ix1)) ix=ix1
c
      ixloc(i)=ix
   10 continue
c  Fill remainder of array with zeroes.
      do 12 i=ixmax,ixkms
   12 ixloc(i)=0
c
c
      iy=1
      do 15 i=1,iymax
c
      iy1=iy+1
c
      ynow=float(i)-yl
      if (ynow.ge.yn(iy1)) iy=iy1
c
      iyloc(i)=iy
   15 continue
c
c  Fill rest of array with zeroes.
      do 17 i=iymax,iykms
 17   iyloc(i)=0
c
      iz=1
      do 20 i=1,izmax
c
      iz1=iz+1
c
      znow=float(i)-zl
      if (znow.ge.zn(iz1)) iz=iz1
c
      izloc(i)=iz
   20 continue
c
c  Fill remainder of array with zeroes.
      do 22 i=izmax,izkms
  22  izloc(i)=0
      return
 330   continue
      write(16,331)ixkms,iykms,izkms
 331  format(' ***** error in array size in common/locate/',/,
     *' maximum map dimensions (km)=',/,' x=',i5,' y=',i5,' z=',i5)
      write(16,332)ixmax,iymax,izmax
  332 format(' Actual map size (km): ',/,' x=',i5,' y=',i5,' z=',i5)
      stop
c***** end of subroutine bldmap *****
      end
c
      subroutine disto(dlat,rlat,dlon,rlon,x,y)
c  this routine calculates distance of station or event
c  from given coordinate origin in terms of (possibly
c  rotated) cartesian coords x and y
c  uses short distance conversion factors from setorg
      common/short1/ xltkm1,xlnkm1,rota1,xlt1,xln1,snr1,csr1
c   calculate coords of station n w.r.t.
c  center of coords, doing rotation if any
      double precision drad,drlt
      parameter (drad=1.7453292d-02)
      parameter (drlt=9.9330647d-01)
c     plt=60.*ltds(n)+sltm(n)
c     pln=60.*lnds(n)+slnm(n)
      plt=60.*dlat+rlat
      pln=60.*dlon+rlon
c  now convert lat and lon differences to km
      x=pln-xln1
      y=plt-xlt1
      xlta=atan(drlt*tan(drad*(plt+xlt1)/120.))
      x=x*xlnkm1*cos(xlta)
      y=y*xltkm1
c  now do rotation
      if(rota1.eq.0.0) return
      ty=csr1*y+snr1*x
      x=csr1*x-snr1*y
      y=ty
      return
c***** end of subroutine disto *****
      end
c
c
      subroutine input3
c  this routine reads in the initial velocity model in the
c  form of velocity specified on a uniform but not
c  necessarily evenly spaced grid of points
c  (reads from file03 )
c  common block variables:
      common/vmod3d/iuses,nx,ny,nz,xn(30),yn(30),zn(30),vel(30,30,40)
      common/vmodnw/ mx,my,mz,xm(30),ym(30),zm(30),vnew(30,30,40)
        common/arrays/iastns,iaeqs,ianod,ianodi,iandst,iasltn,iaobs
      common/step/istep
      integer ixf(3200),iyf(3200),izf(3200)
      character*1 vtype(2)
      parameter(zero=0.0,izero=0)
      vtype(1)='P'
      vtype(2)='S'
c
c  for this version the gridpoints can be unevenly spaced
c  the origin of the coordinate system is at (x,y,z)=(0,0,0)
c  which will not in general correspond to the point
c  xn(1),yn(1),zn(1).
c
c input the number of gridpoints in x, y and z directions
      read(3,3023) istep,nx,ny,nz, iuses
        atemp=iuses*(nx-2)*(ny-2)*(nz-2)
        if(atemp.le.ianod)goto 40
        write(6,*)'atemp=',atemp,'ianod=',ianod
        write(16,42)
 42     format('0Too many nodes for program array sizes.')
        write(16,*)'atemp=',atemp,'ianod=',ianod
        stop
 40     continue
      write(16,3005) nx,ny,nz,iuses
c
c  input the x grid, y grid, and z grid
        read(3,3004) (xn(i),i=1,nx)
      write(16,3006) (xn(i),i=1,nx)
        read(3,3004) (yn(i),i=1,ny)
      write(16,3007) (yn(i),i=1,ny)
        read(3,3004) (zn(i),i=1,nz)
      write(16,3008) (zn(i),i=1,nz)
 3023 format(f4.1,4i3)
 3003 format(f4.2,4i3)
 3004 format(30f6.1)
c  read dummy line from fixed nodes 
      read(3,3003) xidum
c
 3005 format(//,' velocity grid size:',/,
     * ' nx =',i3,5x,'ny =',i3,5x,'nz =',i3,5x,'iuses=',i3)
c
 3006 format(/,' xgrid',/,3x,30f7.2)
 3007 format(/,' ygrid',/,3x,30f7.2)
 3008 format(/,' zgrid',/,3x,30f7.2,/)
c
c
c  now read in the velocity values
   65 write(16,3101)
      do 38 kv=1,iuses
        do 36 k=1,nz
        k2=k + (kv-1)*nz
        write(16,3015) k,vtype(kv)
        do 36 j=1,ny
        read(3,3011) (vel(i,j,k2),i=1,nx)
        write(16,3013) (vel(i,j,k2),i=1,nx)
   36   continue
   38 continue
 3013 format(30f7.2)
 3015 format(/,' layer',i3,5x,a1,' velocity')
 3011 format(30f5.2)
 3101 format(//,' velocity values on three-dimensional grid')
c
c
c  Set up an array which is used to point to node indices, for any x,y,z
      call bldmap
c
c  read in new velocity grid
c
c input the number of gridpoints in x, y and z directions
      read(13,3003) xdummy,mx,my,mz
        atemp=iuses*(mx-2)*(my-2)*(mz-2)
        if(atemp.le.ianod)goto 80
        write(16,42)
        write(16,*)'22 atemp=',atemp,'ianod=',ianod
      stop
c
c  input the x grid, y grid, and z grid
   80   read(13,3004) (xm(i),i=1,mx)
        read(13,3004) (ym(i),i=1,my)
        read(13,3004) (zm(i),i=1,mz)
c
 2005 format(//,' velocity grid size:',/,
     * ' nx =',i3,5x,'ny =',i3,5x,'nz =',i3)
c
      write(16,2006) (xm(i),i=1,mx)
 2006 format(/,' New xgrid',/,3x,30f7.2)
      write(16,2007) (ym(i),i=1,my)
 2007 format(/,' New ygrid',/,3x,30f7.2)
      write(16,2008) (zm(i),i=1,mz)
 2008 format(/,' New zgrid',/,3x,30f7.2,/)
      return
c
  980 continue
      write(16,1698) nparv,ianod
 1698 format(/,'  ****** STOP ******',/,i8,' velocity nodes, program',
     2 ' arrays only allow',i6)
      stop
  990 continue
      write(16,1699) npari,iandst
 1699 format(/,'  ****** STOP ******',/,i8,' parameters to invert for',
     2 ', program arrays only allow',i6)
      stop
c
c***** end of subroutine input3 *****
      end
c
      subroutine intmap(x,y,z,ip,jp,kp)
c  Modified by W. Prothero so a single call can get the indices
      integer ixloc,iyloc,izloc
      common/locate/xl,yl,zl,ixloc(6500),iyloc(6500),izloc(6500)
      ip=int(x+xl)
      ip=ixloc(ip)
      jp=int(yl+y)
      jp=iyloc(jp)
      kp=int(z+zl)
      kp=izloc(kp)
c  If an array element=0, the position is off the map.
      return
c***** end of subroutine intmap *****
      end
c
c
       subroutine latlon(x,y,lat,xlat,lon,xlon)
c  Subroutine to convert from Cartesian coordinates back to
c  latitude and longitude.
c
      common/short2/ xltkm2,xlnkm2,rota2,xlt2,xln2,snr2,csr2
c
      rad=1.7453292e-2
      rlt=9.9330647e-1
c
      fy=csr2*y-snr2*x
      fx=snr2*y+csr2*x
c
      fy=fy/xltkm2
      plt=xlt2+fy
c
      xlta=atan(rlt*tan(rad*(plt+xlt2)/120.))
      fx=fx/(xlnkm2*cos(xlta))
      pln=xln2+fx
c
      lat=plt/60.
      xlat=plt-lat*60.
      lon=pln/60.
      xlon=pln-lon*60.
c
      return
c
c***** end of subroutine latlon ******
      end
c
c
      subroutine setor1(ifil1)
c  this routine establishes the short distance conversion factors
c  given the origin of coordinates
c  of the old coordinate system
c  the rotation angle is converted to radians also
c  common block variables:
      common/short1/ xltkm1,xlnkm1,rota1,xlt1,xln1,snr1,csr1
c  local variables:
      double precision dlt1,dlt2,dxlt,drad,drlt
      parameter ( re=6378.163, ell=298.26)
      parameter (drad=1.7453292d-2)
      parameter (drlt=9.9330647d-1)
      write(16,2000)
 2000 format('/  origin :  latitude   longitude   rotation')
c  read in center of coordinates and rotation angle
      read(ifil1,*) ltdo,oltm,lndo,olnm,rota1
    5 write(16,2002) ltdo,oltm,lndo,olnm,rota1
 2002 format(12x,i3,1x,f5.2,2x,i4,1x,f5.2,3x,f7.2)
c  set up short-distance conversion factors, given center of coords
      rota1=rota1*drad       !convert from degrees to radians
c
      dxlt=dble(60.*ltdo+oltm)
      xln1=60.*lndo+olnm
      xlt1=sngl(dxlt)
c  conversion factor for latitude
      dlt1=datan(drlt*dtan(dxlt*drad/60.d0))
      dlt2=datan(drlt*dtan((dxlt+1.d0)*drad/60.d0))
      del=sngl(dlt2-dlt1)
      r=re*(1.0-sngl(dsin(dlt1)**2)/ell)
      xltkm1=r*del
c  conversion factor for longitude
      del=sngl(dacos(1.0d0-(1.0d0-dcos(drad/60.d0))*dcos(dlt1)**2))
      bc=r*del
      xlnkm1=bc/sngl(dcos(dlt1))
      write(16,3001) xltkm1,bc
 3001 format(5x,'short distance conversion factors',/,
     2       10x,'one min lat  ',f7.4,' km',/,
     3       10x,'one min lon  ',f7.4,' km',/)
c
c  convert coordinates with rotation cosines
      snr1=sin(rota1)
      csr1=cos(rota1)
      return
c***** end of subroutine setor1 *****
      end
c
      subroutine setor2(ifil2)
c  this routine establishes the short distance conversion factors
c  given the origin of coordinates
c  of the old coordinate system
c  the rotation angle is converted to radians also
c  common block variables:
      common/short2/ xltkm2,xlnkm2,rota2,xlt2,xln2,snr2,csr2
c  local variables:
      double precision dlt1,dlt2,dxlt,drad,drlt
      parameter ( re=6378.163, ell=298.26)
      parameter (drad=1.7453292d-2)
      parameter (drlt=9.9330647d-1)
      write(16,2000)
 2000 format('/  origin :  latitude   longitude   rotation')
c  read in center of coordinates and rotation angle
      read(ifil2,*) ltdo,oltm,lndo,olnm,rota2
    5 write(16,2002) ltdo,oltm,lndo,olnm,rota2
 2002 format(12x,i3,1x,f5.2,2x,i4,1x,f5.2,3x,f7.2)
c  set up short-distance conversion factors, given center of coords
      rota2=rota2*drad       !convert from degrees to radians
c
      dxlt=dble(60.*ltdo+oltm)
      xln2=60.*lndo+olnm
      xlt2=sngl(dxlt)
c  conversion factor for latitude
      dlt1=datan(drlt*dtan(dxlt*drad/60.d0))
      dlt2=datan(drlt*dtan((dxlt+1.d0)*drad/60.d0))
      del=sngl(dlt2-dlt1)
      r=re*(1.0-sngl(dsin(dlt1)**2)/ell)
      xltkm2=r*del
c  conversion factor for longitude
      del=sngl(dacos(1.0d0-(1.0d0-dcos(drad/60.d0))*dcos(dlt1)**2))
      bc=r*del
      xlnkm2=bc/sngl(dcos(dlt1))
      write(16,3001) xltkm2,bc
 3001 format(5x,'short distance conversion factors',/,
     2       10x,'one min lat  ',f7.4,' km',/,
     3       10x,'one min lon  ',f7.4,' km',/)
c
c  convert coordinates with rotation cosines
      snr2=sin(rota2)
      csr2=cos(rota2)
      return
c***** end of subroutine setor2 *****
      end
c
      subroutine vel3(isp,x,y,z,v)
c  This routine is Cliff Thurber's
      common/vmod3d/iuses,nx,ny,nz,xn(30),yn(30),zn(30),vel(30,30,40)
      common/locate/ xl,yl,zl,ixloc(6500),iyloc(6500),izloc(6500)
c  use Prothero's intmap here
      common/weight/ wv(8),ip,jp,kp,kpg
      call intmap(x,y,z,ip,jp,kp)
c
      ip1=ip+1
      jp1=jp+1
      kp1=kp+1
c
c    	write(16,100)x,xl,ip,y,yl,jp,z,zl,kp
c100 	format(3(2f11.3,i3))
      xf=(x-xn(ip))/(xn(ip1)-xn(ip))
      yf=(y-yn(jp))/(yn(jp1)-yn(jp))
      zf=(z-zn(kp))/(zn(kp1)-zn(kp))
      xf1=1.0-xf
      yf1=1.0-yf
      zf1=1.0-zf
c
      wv(1)=xf1*yf1*zf1
      wv(2)=xf*yf1*zf1
      wv(3)=xf1*yf*zf1
      wv(4)=xf*yf*zf1
      wv(5)=xf1*yf1*zf
      wv(6)=xf*yf1*zf
      wv(7)=xf1*yf*zf
      wv(8)=xf*yf*zf
c  calculate velocity
c  S-velocity is stored after P-velocity
      kpg=kp
      if(isp.eq.1) kp=kp+nz
      kp1=kp+1
      v=wv(1)*vel(ip,jp,kp)+wv(2)*vel(ip1,jp,kp)
     2 +wv(3)*vel(ip,jp1,kp)+wv(4)*vel(ip1,jp1,kp)
     * +wv(5)*vel(ip,jp,kp1)+wv(6)*vel(ip1,jp,kp1)
     * +wv(7)*vel(ip,jp1,kp1)+wv(8)*vel(ip1,jp1,kp1)
      return
c***** end of subroutine vel3 *****
      end
c
c