      SUBROUTINE linterp2d (LBx, UBx, LBy, UBy,
     &                      Xinp, Yinp, Finp,
     &                      Istr, Iend, Jstr, Jend,
     &                      Xout, Yout, Fout)
      implicit none
      integer*4  LLm,Lm,MMm,Mm,N, LLm0,MMm0
      parameter (LLm0=41,   MMm0=42,   N=32)
      parameter (LLm=LLm0,  MMm=MMm0)
      integer*4 Lmmpi,Mmmpi,iminmpi,imaxmpi,jminmpi,jmaxmpi
      common /comm_setup_mpi1/ Lmmpi,Mmmpi
      common /comm_setup_mpi2/ iminmpi,imaxmpi,jminmpi,jmaxmpi
      integer*4 NSUB_X, NSUB_E, NPP
      integer*4 NP_XI, NP_ETA, NNODES
      parameter (NP_XI=2,  NP_ETA=2,  NNODES=NP_XI*NP_ETA)
      parameter (NPP=1)
      parameter (NSUB_X=1, NSUB_E=1)
      integer*4 NWEIGHT
      parameter (NWEIGHT=1000)
      integer*4 stdout, Np, padd_X,padd_E
      common /stdout/stdout
      parameter (Np=N+1)
      parameter (Lm=(LLm+NP_XI-1)/NP_XI, Mm=(MMm+NP_ETA-1)/NP_ETA)
      parameter (padd_X=(Lm+2)/2-(Lm+1)/2)
      parameter (padd_E=(Mm+2)/2-(Mm+1)/2)
      integer*4 NSA, N2d,N3d, size_XI,size_ETA
      integer*4 se,sse, sz,ssz
      parameter (NSA=28)
      parameter (size_XI=7+(Lm+NSUB_X-1)/NSUB_X)
      parameter (size_ETA=7+(Mm+NSUB_E-1)/NSUB_E)
      parameter (sse=size_ETA/Np, ssz=Np/size_ETA)
      parameter (se=sse/(sse+ssz), sz=1-se)
      parameter (N2d=size_XI*(se*size_ETA+sz*Np))
      parameter (N3d=size_XI*size_ETA*Np)
      real Vtransform
      parameter (Vtransform=2)
      integer*4   NT, itemp
      integer*4   ntrc_salt, ntrc_pas, ntrc_bio, ntrc_sed
      parameter (itemp=1)
      parameter (ntrc_salt=1)
      parameter (ntrc_pas=0)
      parameter (ntrc_bio=0)
      parameter (ntrc_sed=0)
      parameter (NT=itemp+ntrc_salt+ntrc_pas+ntrc_bio+ntrc_sed)
      integer*4   ntrc_diats, ntrc_diauv, ntrc_diabio
      integer*4   ntrc_diavrt, ntrc_diaek, ntrc_diapv
      integer*4   ntrc_diaeddy, ntrc_surf
     &          , isalt
      parameter (isalt=itemp+1)
      parameter (ntrc_diabio=0)
      parameter (ntrc_diats=0)
      parameter (ntrc_diauv=0)
      parameter (ntrc_diavrt=0)
      parameter (ntrc_diaek=0)
      parameter (ntrc_diapv=0)
      parameter (ntrc_diaeddy=0)
      parameter (ntrc_surf=0)
      integer*4, intent(in) :: LBx, UBx, LBy, UBy
      integer*4, intent(in) :: Istr, Iend, Jstr, Jend
      real(kind=8), intent(in) :: Xinp(LBx:UBx)
      real(kind=8), intent(in) :: Yinp(LBy:UBy)
      real(kind=8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
      real(kind=8), intent(in) :: Xout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real(kind=8), intent(in) :: Yout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real(kind=8), intent(out) :: Fout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      integer*4 i, i1, i2, j, j1, j2, ii, jj
      real(kind=8) cff, x, x1, x2, y, y1, y2
      DO j=Jstr,Jend
        DO i=Istr,Iend
           DO ii=LBx,(UBx-1)
             if ((Xinp(ii).le.Xout(i,j)).and.
     &           (Xinp(ii+1).gt.Xout(i,j))) then
               i1=ii
               i2=ii+1
               goto 10
             endif
           enddo
           print*, 'Did not find i1 and i2'
           goto 100
 10        continue
           DO jj=LBy,(UBy-1)
             if ((Yinp(jj).le.Yout(i,j)).and.
     &           (Yinp(jj+1).gt.Yout(i,j))) then
               j1=jj
               j2=jj+1
               goto 20
             endif
           enddo
           print*, 'Did not find j1 and j2'
           goto 100
 20        continue
          IF (((LBx.le.i1).and.(i1.le.UBx)).and.
     &        ((LBy.le.j1).and.(j1.le.UBy))) THEN
            x1=Xinp(i1)
            x2=Xinp(i2)
            y1=Yinp(j1)
            y2=Yinp(j2)
            x=Xout(i,j)
            y=Yout(i,j)
            cff= Finp(i1,j1)*(x2-x )*(y2-y )
     &          +Finp(i2,j1)*(x -x1)*(y2-y )
     &          +Finp(i1,j2)*(x2-x )*(y -y1)
     &          +Finp(i2,j2)*(x -x1)*(y -y1)
            Fout(i,j)=cff/((x2-x1)*(y2-y1))
          END IF
        END DO
      END DO
      RETURN
 100  continue
      print*, 'error in linterp2d'
      END SUBROUTINE linterp2d
      SUBROUTINE cinterp2d (LBx, UBx, LBy, UBy,
     &                      Xinp, Yinp, Finp,
     &                      Istr, Iend, Jstr, Jend,
     &                      Xout, Yout, Fout)
      implicit none
      integer*4  LLm,Lm,MMm,Mm,N, LLm0,MMm0
      parameter (LLm0=41,   MMm0=42,   N=32)
      parameter (LLm=LLm0,  MMm=MMm0)
      integer*4 Lmmpi,Mmmpi,iminmpi,imaxmpi,jminmpi,jmaxmpi
      common /comm_setup_mpi1/ Lmmpi,Mmmpi
      common /comm_setup_mpi2/ iminmpi,imaxmpi,jminmpi,jmaxmpi
      integer*4 NSUB_X, NSUB_E, NPP
      integer*4 NP_XI, NP_ETA, NNODES
      parameter (NP_XI=2,  NP_ETA=2,  NNODES=NP_XI*NP_ETA)
      parameter (NPP=1)
      parameter (NSUB_X=1, NSUB_E=1)
      integer*4 NWEIGHT
      parameter (NWEIGHT=1000)
      integer*4 stdout, Np, padd_X,padd_E
      common /stdout/stdout
      parameter (Np=N+1)
      parameter (Lm=(LLm+NP_XI-1)/NP_XI, Mm=(MMm+NP_ETA-1)/NP_ETA)
      parameter (padd_X=(Lm+2)/2-(Lm+1)/2)
      parameter (padd_E=(Mm+2)/2-(Mm+1)/2)
      integer*4 NSA, N2d,N3d, size_XI,size_ETA
      integer*4 se,sse, sz,ssz
      parameter (NSA=28)
      parameter (size_XI=7+(Lm+NSUB_X-1)/NSUB_X)
      parameter (size_ETA=7+(Mm+NSUB_E-1)/NSUB_E)
      parameter (sse=size_ETA/Np, ssz=Np/size_ETA)
      parameter (se=sse/(sse+ssz), sz=1-se)
      parameter (N2d=size_XI*(se*size_ETA+sz*Np))
      parameter (N3d=size_XI*size_ETA*Np)
      real Vtransform
      parameter (Vtransform=2)
      integer*4   NT, itemp
      integer*4   ntrc_salt, ntrc_pas, ntrc_bio, ntrc_sed
      parameter (itemp=1)
      parameter (ntrc_salt=1)
      parameter (ntrc_pas=0)
      parameter (ntrc_bio=0)
      parameter (ntrc_sed=0)
      parameter (NT=itemp+ntrc_salt+ntrc_pas+ntrc_bio+ntrc_sed)
      integer*4   ntrc_diats, ntrc_diauv, ntrc_diabio
      integer*4   ntrc_diavrt, ntrc_diaek, ntrc_diapv
      integer*4   ntrc_diaeddy, ntrc_surf
     &          , isalt
      parameter (isalt=itemp+1)
      parameter (ntrc_diabio=0)
      parameter (ntrc_diats=0)
      parameter (ntrc_diauv=0)
      parameter (ntrc_diavrt=0)
      parameter (ntrc_diaek=0)
      parameter (ntrc_diapv=0)
      parameter (ntrc_diaeddy=0)
      parameter (ntrc_surf=0)
      real h(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real hinv(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real f(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real fomn(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /grid_h/h /grid_hinv/hinv /grid_f/f /grid_fomn/fomn
      real angler(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /grid_angler/angler
      real latr(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real lonr(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real latu(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real lonu(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real latv(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real lonv(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /grid_latr/latr /grid_lonr/lonr
      common /grid_latu/latu /grid_lonu/lonu
      common /grid_latv/latv /grid_lonv/lonv
      real pm(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pn(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real om_r(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real on_r(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real om_u(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real on_u(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real om_v(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real on_v(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real om_p(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real on_p(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pn_u(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pm_v(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pm_u(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pn_v(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /metrics_pm/pm    /metrics_pn/pn
      common /metrics_omr/om_r /metrics_on_r/on_r
      common /metrics_omu/om_u /metrics_on_u/on_u
      common /metrics_omv/om_v /metrics_on_v/on_v
      common /metrics_omp/om_p /metrics_on_p/on_p
      common /metrics_pnu/pn_u /metrics_pmv/pm_v
      common /metrics_pmu/pm_u /metrics_pnv/pn_v
      real dmde(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real dndx(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /metrics_dmde/dmde    /metrics_dndx/dndx
      real pmon_p(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pmon_r(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pmon_u(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pnom_p(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pnom_r(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pnom_v(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real grdscl(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /metrics_pmon_p/pmon_p /metrics_pnom_p/pnom_p
      common /metrics_pmon_r/pmon_r /metrics_pnom_r/pnom_r
      common /metrics_pmon_u/pmon_u /metrics_pnom_v/pnom_v
      common /metrics_grdscl/grdscl
      real rmask(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pmask(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real umask(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real vmask(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real pmask2(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /mask_r/rmask
      common /mask_p/pmask
      common /mask_u/umask
      common /mask_v/vmask
      common /mask_p2/pmask2
      real zob(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      common /Z0B_VAR/zob
      integer*4, intent(in) :: LBx, UBx, LBy, UBy
      integer*4, intent(in) :: Istr, Iend, Jstr, Jend
      real(kind=8), intent(in) :: Xinp(LBx:UBx)
      real(kind=8), intent(in) :: Yinp(LBy:UBy)
      real(kind=8), intent(in) :: Finp(LBx:UBx,LBy:UBy)
      real(kind=8), intent(in) :: Xout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real(kind=8), intent(in) :: Yout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      real(kind=8), intent(out) :: Fout(-1:Lm+2+padd_X,-1:Mm+2+padd_E)
      integer*4 i, ic, iter, i1, i2, j, jc, j1, j2, ii, jj
      real(kind=8) :: a11, a12, a21, a22
      real(kind=8) :: e11, e12, e21, e22
      real(kind=8) :: cff, d1, d2, dfc, dx, dy, eta, xi, xy, yx
      real(kind=8) :: f0, fx, fxx, fxxx, fxxy, fxy, fxyy, fy, fyy, fyyy
      real(kind=8), parameter :: C01 = 1.0D0/48.0D0
      real(kind=8), parameter :: C02 = 1.0D0/32.0D0
      real(kind=8), parameter :: C03 = 0.0625D0
      real(kind=8), parameter :: C04 = 1.0D0/6.0D0
      real(kind=8), parameter :: C05 = 0.25D0
      real(kind=8), parameter :: C06 = 0.5D0
      real(kind=8), parameter :: C07 = 0.3125D0
      real(kind=8), parameter :: C08 = 0.625D0
      real(kind=8), parameter :: C09 = 1.5D0
      real(kind=8), parameter :: C10 = 13.0D0/24.0D0
      real(kind=8), parameter :: LIMTR = 3.0D0
      real(kind=8), parameter :: spv = 0.0D0
      real(kind=8), dimension(-1:2,-1:2) :: dfx, dfy, ff
      DO j=Jstr,Jend
        DO i=Istr,Iend
           DO ii=LBx,(UBx-1)
             if ((Xinp(ii).le.Xout(i,j)).and.
     &           (Xinp(ii+1).gt.Xout(i,j))) then
               i1=ii
               i2=ii+1
               goto 10
             endif
           enddo
           print*, 'Did not find i1 and i2',
     &           Istr,Iend,Jstr,Jend,i,j,Xout(i,j),Xout(i-1,j)
           goto 100
 10        continue
           DO jj=LBy,UBy-1
             if ((Yinp(jj).le.Yout(i,j)).and.
     &           (Yinp(jj+1).gt.Yout(i,j))) then
               j1=jj
               j2=jj+1
               goto 20
             endif
           enddo
           print*, 'Did not find j1 and j2'
           goto 100
 20        continue
          IF (((LBx.le.i1).and.(i1.le.UBx)).and.
     &        ((LBy.le.j1).and.(j1.le.UBy))) THEN
            xy=Xinp(i2)-Xinp(i1)-Xinp(i2)+Xinp(i1)
            yx=Yinp(j2)-Yinp(j2)-Yinp(j1)+Yinp(j1)
            dx=Xout(i,j)-0.25D0*(Xinp(i2)+Xinp(i1)+
     &                         Xinp(i2)+Xinp(i1))
            dy=Yout(i,j)-0.25D0*(Yinp(j2)+Yinp(j2)+
     &                         Yinp(j1)+Yinp(j1))
            e11=0.5D0*(Xinp(i2)-Xinp(i1)+Xinp(i2)-Xinp(i1))
            e12=0.5D0*(Xinp(i2)+Xinp(i1)-Xinp(i2)-Xinp(i1))
            e21=0.5D0*(Yinp(j2)-Yinp(j2)+Yinp(j1)-Yinp(j1))
            e22=0.5D0*(Yinp(j2)+Yinp(j2)-Yinp(j1)-Yinp(j1))
            cff=1.0D0/(e11*e22-e12*e21)
            xi=cff*(e22*dx-e12*dy)
            eta=cff*(e11*dy-e21*dx)
            DO iter=1,4
              d1=dx-e11*xi-e12*eta-xy*xi*eta
              d2=dy-e21*xi-e22*eta-yx*xi*eta
              a11=e11+xy*eta
              a12=e12+xy*xi
              a21=e21+yx*eta
              a22=e22+yx*xi
              cff=1.0D0/(a11*a22-a12*a21)
              xi =xi +cff*(a22*d1-a12*d2)
              eta=eta+cff*(a11*d2-a21*d1)
            END DO
            DO jc=-1,2
              DO ic=-1,2
                ff(ic,jc)=Finp(MAX(1,MIN(UBx,i1+ic)),
     &                         MAX(1,MIN(UBy,j1+jc)))
              END DO
            END DO
            f0=C07*(ff(1,1)+ff(1,0)+ff(0,1)+ff(0,0))-
     &         C02*(ff(2,0)+ff(2,1)+ff(1,2)+ff(0,2)+
     &              ff(-1,1)+ff(-1,0)+ff(0,-1)+ff(1,-1))
            fx=C08*(ff(1,1)+ff(1,0)-ff(0,1)-ff(0,0))-
     &         C01*(ff(2,1)+ff(2,0)-ff(-1,1)-ff(-1,0))-
     &         C03*(ff(1,2)-ff(0,2)+ff(1,-1)-ff(0,-1))
            fy=C08*(ff(1,1)-ff(1,0)+ff(0,1)-ff(0,0))-
     &         C01*(ff(1,2)+ff(0,2)-ff(1,-1)-ff(0,-1))-
     &         C03*(ff(2,1)-ff(2,0)+ff(-1,1)-ff(-1,0))
            fxy=ff(1,1)-ff(1,0)-ff(0,1)+ff(0,0)
            fxx=C05*(ff(2,1)-ff(1,1)-ff(0,1)+ff(-1,1)+
     &               ff(2,0)-ff(1,0)-ff(0,0)+ff(-1,0))
            fyy=C05*(ff(1,2)-ff(1,1)-ff(1,0)+ff(1,-1)+
     &               ff(0,2)-ff(0,1)-ff(0,0)+ff(0,-1))
            fxxx=C06*(ff(2,1)+ff(2,0)-ff(-1,1)-ff(-1,0))-
     &           C09*(ff(1,1)+ff(1,0)-ff(0,1)-ff(0,0))
            fyyy=C06*(ff(1,2)+ff(0,2)-ff(1,-1)-ff(0,-1))-
     &           C09*(ff(1,1)-ff(1,0)+ff(0,1)-ff(0,0))
            fxxy=C06*(ff(2,1)-ff(1,1)-ff(0,1)+ff(-1,1)-
     &                ff(2,0)+ff(1,0)+ff(0,0)-ff(-1,0))
            fxyy=C06*(ff(1,2)-ff(1,1)-ff(1,0)+ff(1,-1)-
     &                ff(0,2)+ff(0,1)+ff(0,0)-ff(0,-1))
            Fout(i,j)=f0+
     &                fx*xi+
     &                fy*eta+
     &                C06*fxx*xi*xi+
     &                fxy*xi*eta+
     &                C06*fyy*eta*eta+
     &                C04*fxxx*xi*xi*xi+
     &                C06*fxxy*xi*xi*eta+
     &                C04*fyyy*eta*eta*eta+
     &                C06*fxyy*xi*eta*eta
          END IF
        END DO
      END DO
      RETURN
 100  continue
      print*, 'error in cinterp2d'
      END SUBROUTINE cinterp2d
