       program grdintrp
c      program to interpolate 3-D deformation gridded results to a particular latitude and longitude
c      the program uses a bicubic spline subroutine from numerical recipes
c      and longitude
c
c      displacement grid is 1. deg x 1. deg
c      361 longitudes
c      181 latitudes
c
c      dr=radial displacements
c      vt=horizontal (north) displacement
c      vl=horizontal (east) displacement
c
c      the program takes as input the sine and cosine of the displacements (mm):
c      ss1,cs1,ss2,cs2
c
c      the user needs to specify iref(reference frame) and iout(output format)
c      within the program
c      and the program prompts the user for the station latitude and longitude
c
c      rlat and rlon are the latitudes and longitudes of displacement grid 
c
c      
c      program written: september, 2005 T. vanDam 

c      please note: i already know that i am not the most efficient or clear
c      programmer in the world. so you might have to work some things out on your own 

       implicit none
       integer nlon,nlat,i,j,iref,iout,ista
       parameter(nlon=360,nlat=181)
       real*4 ss1dr(nlon,nlat),cs1dr(nlon,nlat),ss2dr(nlon,nlat),
     +  cs2dr(nlon,nlat),
     +  ss1vt(nlon,nlat),cs1vt(nlon,nlat),ss2vt(nlon,nlat),
     +  cs2vt(nlon,nlat),
     +  ss1vl(nlon,nlat),cs1vl(nlon,nlat),ss2vl(nlon,nlat),
     +  cs2vl(nlon,nlat)
       real*4 ss1dri(nlon,nlat),cs1dri(nlon,nlat),ss2dri(nlon,nlat),
     +  cs2dri(nlon,nlat),
     +  ss1vti(nlon,nlat),cs1vti(nlon,nlat),ss2vti(nlon,nlat),
     +  cs2vti(nlon,nlat),
     +  ss1vli(nlon,nlat),cs1vli(nlon,nlat),ss2vli(nlon,nlat),
     +  cs2vli(nlon,nlat)
       real*4 utot(4),vttot(4),vltot(4)
       real*4 uamp1,uamp2,upha1,upha2,vtamp1,vtamp2,vtpha1,vtpha2,
     +  vlamp1,vlamp2,vlpha1,vlpha2
       real*4 rlat(nlat),rlon(nlon),rlt,rln,skp1,skp2,orlt,erln,dtr,
     + corlt
       parameter(dtr=6.2831853071795865/360.)
       character sta*4,filein*48
c      specify reference frame:
c      if reference frame=ce then iref=1 
c      if reference frame=cm then iref=2 
c      if reference frame=cf then iref=3 
       iref=1
       if(iref .eq. 1) filein='s1_s2_def_ce.dat'
       if(iref .eq. 2) filein='s1_s2_def_cm.dat'
       if(iref .eq. 3) filein='s1_s2_def_cf.dat'
       open(unit=11,status='old',file=filein)
       open(unit=12,status='old',file='in.grdintrp')
       open(unit=21,status='unknown',file='grdintrp.dat')
c      check the next line: iout=1 => output sines and cosines; iout=2 output amplitudes and phases
       iout=1
       write(21,*) '$$ Atmospheric Tidal loading displacememt'
       write(21,*) '$$ '
       write(21,*) '$$ calculated using grdinterp.f (t. van Dam)' 
       write(21,*) '$$ calculated using atmospheric tidal model Ray and'
       write(21,*) '$$ Ponte, Annales Geophysicae (2003) v21; 1897-1910'
       write(21,*) '$$ Farrell elastic Greens functions have been used' 
       write(21,*) '$$ '
       write(21,*) '$$ Displacement is defined positive up, north'
       write(21,*) '$$ and east directions'
       write(21,*) '$$ '
       if(iout .eq. 1) then
         write(21,*) '$$ Column order coss1 sins1 coss2 sins2'
         write(21,*) '$$ '
         write(21,*) '$$ Row order:'
         write(21,*) '$$ RADIAL'
         write(21,*) '$$ TANGENTIAL NS'
         write(21,*) '$$ TANGENTIAL EW'
         write(21,*) '$$ '
         write(21,*) '$$ '
         write(21,*) '$$ '
       else
         write(21,*) '$$ Phase lag relative to Greenwich and the lags'
         write(21,*) '$$ are positive'
         write(21,*) '$$ Column order: radial s1, radial s2, north s1,'
         write(21,*) '$$ north s2, east s1, east s2'
         write(21,*) '$$ '
         write(21,*) '$$ Row order:'
         write(21,*) '$$ Amplitude (mm)'
         write(21,*) '$$ Phase (deg)'
         write(21,*) '$$ '
         write(21,*) '$$ '
         write(21,*) '$$ '
       end if
       do i=1,nlat
         rlat(i)=(i-1)*1.
       end do
       rlat(nlat)=179.99
       do i=1,nlon
         rlon(i)=(i-1)*1.
       end do
       do i=1,nlon
         do j=1,nlat
           read(11,*) skp1,skp2,
     +                cs1dr(i,j),ss1dr(i,j),cs2dr(i,j),ss2dr(i,j),
     +                cs1vt(i,j),ss1vt(i,j),cs2vt(i,j),ss2vt(i,j),
     +                cs1vl(i,j),ss1vl(i,j),cs2vl(i,j),ss2vl(i,j)
         end do
       end do
       call splie2(rlon,rlat,ss1dr,nlon,nlat,ss1dri)
       call splie2(rlon,rlat,cs1dr,nlon,nlat,cs1dri)
       call splie2(rlon,rlat,ss2dr,nlon,nlat,ss2dri)
       call splie2(rlon,rlat,cs2dr,nlon,nlat,cs2dri)
       call splie2(rlon,rlat,ss1vt,nlon,nlat,ss1vti)
       call splie2(rlon,rlat,cs1vt,nlon,nlat,cs1vti)
       call splie2(rlon,rlat,ss2vt,nlon,nlat,ss2vti)
       call splie2(rlon,rlat,cs2vt,nlon,nlat,cs2vti)
       call splie2(rlon,rlat,ss1vl,nlon,nlat,ss1vli)
       call splie2(rlon,rlat,cs1vl,nlon,nlat,cs1vli)
       call splie2(rlon,rlat,ss2vl,nlon,nlat,ss2vli)
       call splie2(rlon,rlat,cs2vl,nlon,nlat,cs2vli)
       do ista=1,50000
         read(12,*,end=100) sta,rln,rlt
c        convert latitude to co-latitude
         corlt=90.-rlt
c        tf necessary, convert longitude to east longitude
         erln=rln
         if(rln .lt. 0) erln=rln+360.
         call splin2(rlon,rlat,cs1dr,cs1dri,nlon,nlat,erln,corlt
     +      ,utot(1))
         call splin2(rlon,rlat,ss1dr,ss1dri,nlon,nlat,erln,corlt
     +      ,utot(2))
         call splin2(rlon,rlat,cs2dr,cs2dri,nlon,nlat,erln,corlt
     +      ,utot(3))
         call splin2(rlon,rlat,ss2dr,ss2dri,nlon,nlat,erln,corlt
     +      ,utot(4))
         call splin2(rlon,rlat,cs1vt,cs1vti,nlon,nlat,erln,corlt
     +      ,vttot(1))
         call splin2(rlon,rlat,ss1vt,ss1vti,nlon,nlat,erln,corlt
     +      ,vttot(2))
         call splin2(rlon,rlat,cs2vt,cs2vti,nlon,nlat,erln,corlt
     +      ,vttot(3))
         call splin2(rlon,rlat,ss2vt,ss2vti,nlon,nlat,erln,corlt
     +      ,vttot(4))
         call splin2(rlon,rlat,cs1vl,cs1vli,nlon,nlat,erln,corlt
     +      ,vltot(1))
         call splin2(rlon,rlat,ss1vl,ss1vli,nlon,nlat,erln,corlt
     +      ,vltot(2))
         call splin2(rlon,rlat,cs2vl,cs2vli,nlon,nlat,erln,corlt
     +      ,vltot(3))
         call splin2(rlon,rlat,ss2vl,ss2vli,nlon,nlat,erln,corlt
     +      ,vltot(4))
         if(iout .eq. 1) then
           write(21,*) '$$ '
           write(21,*) '$$ station ',sta,'; coord.(long,lat)',rln,rlt
           write(21,*) '$$ '
           write (21,1000) (utot(i),i=1,4)
           write (21,1000) (vttot(i),i=1,4)
           write (21,1000) (vltot(i),i=1,4)
          else
           uamp1=sqrt(utot(1)**2+utot(2)**2)
           uamp2=sqrt(utot(3)**2+utot(4)**2)
           vtamp1=sqrt(vttot(1)**2+vttot(2)**2)
           vtamp2=sqrt(vttot(3)**2+vttot(4)**2)
           vlamp1=sqrt(vltot(1)**2+vltot(2)**2)
           vlamp2=sqrt(vltot(3)**2+vltot(4)**2)
           upha1=atan2(utot(2),utot(1))/dtr
           upha2=atan2(utot(4),utot(3))/dtr
           if(upha1 .lt. 0) upha1=upha1+360.
           if(upha2 .lt. 0) upha2=upha2+360.
           vtpha1=atan2(vttot(2),vttot(1))/dtr
           vtpha2=atan2(vttot(4),vttot(3))/dtr
           if(vtpha1 .lt. 0) vtpha1=vtpha1+360.
           if(vtpha2 .lt. 0) vtpha2=vtpha2+360.
           vlpha1=atan2(vltot(2),vltot(1))/dtr
           vlpha2=atan2(vltot(4),vltot(3))/dtr
           if(vlpha1 .lt. 0) vlpha1=vlpha1+360.
           if(vlpha2 .lt. 0) vlpha2=vlpha2+360.
           write(21,*) '$$ '
           write(21,*) '$$ station ',sta,'; coord.(long,lat)',rln,rlt
           write(21,*) '$$ '
           write(21,2000) uamp1,uamp2,vtamp1,vtamp2,vlamp1,vlamp2
           write(21,2000) upha1,upha2,vtpha1,vtpha2,vlpha1,vlpha2
          end if
1000      format(12f12.4)
2000      format(6f12.4)
        end do
100     continue
        stop
        end
      SUBROUTINE splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y)
      INTEGER m,n,NN
      REAL x1,x2,y,x1a(m),x2a(n),y2a(m,n),ya(m,n)
      PARAMETER (NN=500)
CU    USES spline,splint
      INTEGER j,k
      REAL y2tmp(NN),ytmp(NN),yytmp(NN)
      do 12 j=1,m
        do 11 k=1,n
          ytmp(k)=ya(j,k)
          y2tmp(k)=y2a(j,k)
11      continue
        call splint(x2a,ytmp,y2tmp,n,x2,yytmp(j))
12    continue
      call spline(x1a,yytmp,m,1.e30,1.e30,y2tmp)
      call splint(x1a,yytmp,y2tmp,m,x1,y)
      return
      END
      SUBROUTINE splie2(x1a,x2a,ya,m,n,y2a)
      INTEGER m,n,NN
      REAL x1a(m),x2a(n),y2a(m,n),ya(m,n)
      PARAMETER (NN=500)
CU    USES spline
      INTEGER j,k
      REAL y2tmp(NN),ytmp(NN)
      do 13 j=1,m
        do 11 k=1,n
          ytmp(k)=ya(j,k)
11      continue
        call spline(x2a,ytmp,n,1.e30,1.e30,y2tmp)
        do 12 k=1,n
          y2a(j,k)=y2tmp(k)
12      continue
13    continue
      return
      END
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
