      program tseries
c     example program to extract the loding time series for STJO
c     for the period 1993-1995, in the CM reference frame

c     BEFORE ALTERING THE PROGRAM...
c     compile and run this program and compare the output with the
c     data in STJO.ntol.orig 

c     the program will need to be altered if different time periods 
c     and/or different sites are desired 

c     user input 
c       reference frame=iref
c       rlt=latitude of point grid data should be interpolated to 
c       rln=longitude of point grid data should be interpolated to 
c     note: the program is currently set up to interpolate the data to stjo
c     STJO  47.60  307.32  

c     output variables: 
c       ns=north displacement 
c       ew=east displacement 
c       up=radial displacement 
c     all dislpacements in mm

c     program assumes directory structure given on the website, i.e. top directories 
c     cf_ggfc/  ce_ggfc/ cm_ggfc/
c     with the year directories inside

      implicit none
      integer ilon,ilat,nlon,nlat
      integer iyr,idy,ihr,nday,iref
      integer i1000,i100,i10,i1,kcnt
      character dig(0:9)*1,chyr*3,cdoy*3,cref(3)*2,
     +      chhr(2)*2,chyrw*4,sta*4
      character*40 infl,outfl
      data dig/'0','1','2','3','4','5','6','7','8','9'/
      data chhr/'06','18'/
      data cref/'ce','cm','cf'/

      parameter(nlon=144,nlat=73)

      real*4 dr(nlon,nlat),dvt(nlon,nlat),
     +    dvl(nlon,nlat)
      real*4 dri(nlon,nlat),dvti(nlon,nlat),
     +    dvli(nlon,nlat)
      real*4 rlt,rln,corlt,erln,time,thr
      real*4 rlon(nlon),rlat(nlat) 
      real*4 ns,ew,up

c     user set reference frame
c     if(iref .eq. 1) ref = ce
c     if(iref .eq. 2) ref = cm
c     if(iref .eq. 3) ref = cf
      iref=2

c     user input point to which data are to be interpolated
      sta='STJO'
      rln=307.32
      rlt=47.60

c     convert latitude to co-latitude
      corlt=90.-rlt
c     if necessary, convert longitude to east longitude
      erln=rln
      if(rln .lt. 0) erln=rln+360.
c     open output file with the station name
      outfl=sta//'.ntol'
      open(unit=21,status='unknown',file=outfl)

      do ilat=1,nlat
        rlat(ilat)=(ilat-1)*2.5
      end do
      do ilon=1,nlon
         rlon(ilon)=(ilon-1)*2.5
      end do

c     read loop over all years

      do iyr=1993,1995
        i1000=iyr/1000
        i100=(iyr-i1000*1000)/100
        i10=(iyr-100*i100-i1000*1000)/10
        i1=(iyr-100*i100-10*i10-i1000*1000)
        chyrw=dig(i1000)//dig(i100)//
     +         dig(i10)//dig(i1)
        nday=365
        if(mod(iyr,4) .eq. 0) nday=366

c       read loop over all days

        do idy=1,nday
         i100=idy/100
         i10=(idy-100*i100)/10
         i1=(idy-100*i100-10*i10)
         cdoy=dig(i100)//dig(i10)//dig(i1)

c        generate file names from iyr and idy

         infl=cref(iref)//'_ggfc/'//chyrw//'/'//chyrw//'_'//cdoy
         open(unit=11,status='unknown',file=infl)
c        read loop over all epochs
         do ihr=1,2
           if(ihr .eq. 2) thr=6.
           if(ihr .eq. 4) thr=18.
           read (11,1002) 
1002       format(a4,1x,a3,1x,a2,1x,a2)
           do ilon=1,nlon
             do ilat=1,nlat
               read(11,1003) dvt(ilon,ilat),
     +         dvl(ilon,ilat),dr(ilon,ilat)
1003           format(3(f10.4,2x))
             end do
          end do

c         set up arrays for the bicubic spline interpolation

c          PLEASE NOTE: the routine splie2 only needs to be called
c          once for each grid; BUT the routine splin2 should be called
c          for each point you would like to interpolate to

           call splie2(rlon,rlat,dvt,nlon,nlat,dvti)
           call splie2(rlon,rlat,dvl,nlon,nlat,dvli)
           call splie2(rlon,rlat,dr,nlon,nlat,dri)


c          perform interpolation
           call splin2(rlon,rlat,dvt,dvti,nlon,nlat,erln,corlt
     +           ,ns)
           call splin2(rlon,rlat,dvl,dvli,nlon,nlat,erln,corlt
     +           ,ew)
           call splin2(rlon,rlat,dr,dri,nlon,nlat,erln,corlt
     +           ,up)
           time=float(iyr)+(float(idy-1)+thr/24.)/float(nday)
           write(21,1004) iyr,idy,ihr,time,ns,ew,up
1004       format(3(2x,i4),2x,f14.7,3(2x,f10.4))
          end do
        end do
      end do
      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
