Subversion Repositories lagranto.ecmwf

      real function int2d(ar,n1,n2,rid,rjd)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 2d-array to an arbitrary
c        location within the grid.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2)
c        n1,n2     int   input   dimensions of ar
c        ri,rj     real  input   grid location to be interpolated to
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2
      real      ar(n1,n2), rid,rjd

c     local declarations
      integer   i,j,ip1,jp1,ih,jh
      real      frac0i,frac0j,frac1i,frac1j,ri,rj

c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      ih=nint(ri)
      jh=nint(rj)

c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif

c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif

      if ((i.eq.ip1).and.(j.eq.jp1)) then
c        no interpolation at all
         int2d=ar(i,j)
      else
         frac0i=ri-float(i)
         frac0j=rj-float(j)
         frac1i=1.-frac0i
         frac1j=1.-frac0j
         int2d = ar(i  ,j  ) * frac1i * frac1j
     &         + ar(i  ,jp1) * frac1i * frac0j
     &         + ar(ip1,j  ) * frac0i * frac1j
     &         + ar(ip1,jp1) * frac0i * frac0j
      endif
      end
      real function int2dm(ar,n1,n2,rid,rjd,misdat)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 2d-array to an arbitrary
c        location within the grid. The interpolation includes the
c        testing of the missing data flag 'misdat'.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2)
c        n1,n2     int   input   dimensions of ar
c        ri,rj     real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested
c     History:
c-----------------------------------------------------------------------
 
c     argument declarations
      integer   n1,n2
      real      ar(n1,n2), rid,rjd, misdat
 
c     local declarations
      integer   i,j,ip1,jp1,ih,jh
      real      frac0i,frac0j,frac1i,frac1j,ri,rj,int2d
 
c     check if routine without missing data checking can be called instead
      if (misdat.eq.0.) then
        int2dm=int2d(ar,n1,n2,rid,rjd)
        return
      endif
 
c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      ih=nint(ri)
      jh=nint(rj)
 
c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif
 
c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif
 
      if ((i.eq.ip1).and.(j.eq.jp1)) then
c        no interpolation at all
         int2dm=ar(i,j)
      else
         if ((misdat.eq.ar(i  ,j  )).or.
     &       (misdat.eq.ar(i  ,jp1)).or.
     &       (misdat.eq.ar(ip1,j  )).or.
     &       (misdat.eq.ar(ip1,jp1))) then
           int2dm=misdat
         else
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int2dm = ar(i  ,j  ) * frac1i * frac1j
     &            + ar(i  ,jp1) * frac1i * frac0j
     &            + ar(ip1,j  ) * frac0i * frac1j
     &            + ar(ip1,jp1) * frac0i * frac0j
         endif
      endif
      end
      real function int2dp(ar,n1,n2,rid,rjd)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 2d-array to an arbitrary
c        location within the grid. The 2d-array is periodic in the
c        i-direction: ar(0,.)=ar(n1,.) and ar(1,.)=ar(n1+1,.).
c        Therefore rid can take values in the range 0,...,n1+1
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2)
c        n1,n2     int   input   dimensions of ar
c        ri,rj     real  input   grid location to be interpolated to
c     History:
c-----------------------------------------------------------------------
 
c     argument declarations
      integer   n1,n2
      real      ar(0:n1+1,n2), rid,rjd
 
c     local declarations
      integer   i,j,ip1,jp1,ih,jh
      real      frac0i,frac0j,frac1i,frac1j,ri,rj
 
c     do linear interpolation
      ri=amax1(0.,amin1(float(n1+1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      ih=nint(ri)
      jh=nint(rj)
 
c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-5) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif
 
c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-5) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif
 
      if ((i.eq.ip1).and.(j.eq.jp1)) then
c        no interpolation at all
         int2dp=ar(i,j)
      else
         frac0i=ri-float(i)
         frac0j=rj-float(j)
         frac1i=1.-frac0i
         frac1j=1.-frac0j
         int2dp = ar(i  ,j  ) * frac1i * frac1j
     &         + ar(i  ,jp1) * frac1i * frac0j
     &         + ar(ip1,j  ) * frac0i * frac1j
     &         + ar(ip1,jp1) * frac0i * frac0j
      endif
      end
      real function cint2d(ar,n1,n2,n3,rid,rjd,ikd)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location in the horizontal. ikd specifies the level (must
c        be an integer). A bicubic method is applied (following
c        the Numerical Recipes).
c     Arguments:
c        ar           real  input   field, define as ar(n1,n2,n3)
c        n1,n2        int   input   dimensions of ar
c        rid,rjd      real  input   grid location to be interpolated to
c        ikd          int   input   level for interpolation
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3,ikd
      real      ar(n1,n2,n3), rid,rjd

c     local declarations
      integer   i,j,k
      real      y(4),y1(4),y2(4),y12(4)

c     indices of lower left corner of interpolation grid

      i=int(rid)
      j=int(rjd)
      k=ikd

c     do not interpolate if i or j are at the boundaries

      if ((i.eq.1).or.(j.eq.1).or.(i.eq.n1).or.(j.eq.n2)) then
        cint2d=ar(i,j,k)
        return
      endif

c     define the arrays y,y1,y2,y12 necessary for the bicubic
c     interpolation (following the Numerical Recipes).

      y(1)=ar(i,j,k)
      y(2)=ar(i+1,j,k)
      y(3)=ar(i+1,j+1,k)
      y(4)=ar(i,j+1,k)

      y1(1)=-(ar(i-1,j,k)-ar(i+1,j,k))/2.
      y1(2)=-(ar(i,j,k)-ar(i+2,j,k))/2.
      y1(3)=-(ar(i,j+1,k)-ar(i+2,j+1,k))/2.
      y1(4)=-(ar(i-1,j+1,k)-ar(i+1,j+1,k))/2.

      y2(1)=-(ar(i,j-1,k)-ar(i,j+1,k))/2.
      y2(2)=-(ar(i+1,j-1,k)-ar(i+1,j+1,k))/2.
      y2(3)=-(ar(i+1,j,k)-ar(i+1,j+2,k))/2.
      y2(4)=-(ar(i,j,k)-ar(i,j+2,k))/2.

      y12(1)=(ar(i+1,j+1,k)-ar(i+1,j-1,k)-ar(i-1,j+1,k)+ar(i-1,j-1,k))/4.
      y12(2)=(ar(i+2,j+1,k)-ar(i+2,j-1,k)-ar(i,j+1,k)+ar(i,j-1,k))/4.
      y12(3)=(ar(i+2,j+2,k)-ar(i+2,j,k)-ar(i,j+2,k)+ar(i,j,k))/4.
      y12(4)=(ar(i+1,j+2,k)-ar(i+1,j,k)-ar(i-1,j+2,k)+ar(i-1,j,k))/4.

      call bcuint(y,y1,y2,y12,i,j,rid,rjd,cint2d)
      return
      end
      real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd

c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk

c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)

c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif

c     Check for interpolation in j
      if (abs(float(jh)-rj).lt.1.e-3) then
        j  =jh
        jp1=jh
      else
        j =min0(int(rj),n2-1)
        jp1=j+1
      endif

c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif

      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           int3d=ar(i,j,k)
c          print *,'int3d 00: ',rid,rjd,rkd,int3d
        else
c          horizontal interpolation only
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j
     &           + ar(i  ,jp1,k  ) * frac1i * frac0j
     &           + ar(ip1,j  ,k  ) * frac0i * frac1j
     &           + ar(ip1,jp1,k  ) * frac0i * frac0j
c          print *,'int3d 10: ',rid,rjd,rkd,int3d
        endif
      else 
        frac0k=rk-float(k)
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           int3d = ar(i  ,j  ,k  ) * frac1k
     &           + ar(i  ,j  ,kp1) * frac0k
c          print *,'int3d 01: ',rid,rjd,rkd,int3d
        else
c          full 3d interpolation
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
     &           + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k 
     &           + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
     &           + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
     &           + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
     &           + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k 
     &           + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
     &           + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c          print *,'int3d 11: ',rid,rjd,rkd,int3d
        endif
      endif
      end
      real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The interpolation includes the 
c        testing of the missing data flag 'misdat'. 
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat

c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d

c     check if routine without missing data checking can be called instead
      if (misdat.eq.0.) then
        int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)
        return
      endif

c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)

c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif

c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif

c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif

      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           if (misdat.eq.ar(i,j,k)) then
             int3dm=misdat
           else
             int3dm=ar(i,j,k)
           endif
c          print *,'int3dm 00: ',rid,rjd,rkd,int3dm
        else
c          horizontal interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  ))) then
             int3dm=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
c            print *,'int3dm 10: ',rid,rjd,rkd,int3dm
           endif
        endif
      else 
        frac0k=rk-float(k)
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1))) then
             int3dm=misdat
           else
             int3dm = ar(i  ,j  ,k  ) * frac1k
     &              + ar(i  ,j  ,kp1) * frac0k
c            print *,'int3dm 01: ',rid,rjd,rkd,int3dm
           endif
        else
c          full 3d interpolation
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
     &         (misdat.eq.ar(ip1,jp1,kp1))) then
             int3dm=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
     &              + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
     &              + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
     &              + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
     &              + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c            print *,'int3dm 11: ',rid,rjd,rkd,int3dm
           endif
        endif
      endif
      end
      real function int3dmlog(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The interpolation includes the
c        testing of the missing data flag 'misdat'.
c        Prior to vertical interpolations the log is taken from the array.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat

c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d

c     print*,'hallo in SR int3dmlog'
c     check if routine without missing data checking can be called instead
      if (misdat.eq.0.) then
        int3dmlog=int3d(ar,n1,n2,n3,rid,rjd,rkd)
        return
      endif

c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)

c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif

c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif

c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif

      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           if (misdat.eq.ar(i,j,k)) then
             int3dmlog=misdat
           else
             int3dmlog=ar(i,j,k)
           endif
c          print *,'int3dmlog 00: ',rid,rjd,rkd,int3dmlog
        else
c          horizontal interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  ))) then
             int3dmlog=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dmlog = ar(i  ,j  ,k  ) * frac1i * frac1j
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
c            print *,'int3dmlog 10: ',rid,rjd,rkd,int3dmlog
           endif
        endif
      else
        frac0k=rk-float(k)
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1))) then
             int3dmlog=misdat
           else
             int3dmlog = log(ar(i  ,j  ,k  )) * frac1k
     &                 + log(ar(i  ,j  ,kp1)) * frac0k
             int3dmlog = exp(int3dmlog)
c            print *,'int3dmlog 01: ',rid,rjd,rkd,int3dmlog
           endif
        else
c          full 3d interpolation
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
     &         (misdat.eq.ar(ip1,jp1,kp1))) then
             int3dmlog=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dmlog = log(ar(i  ,j  ,k  )) * frac1i * frac1j * frac1k
     &                 + log(ar(i  ,jp1,k  )) * frac1i * frac0j * frac1k
     &                 + log(ar(ip1,j  ,k  )) * frac0i * frac1j * frac1k
     &                 + log(ar(ip1,jp1,k  )) * frac0i * frac0j * frac1k
     &                 + log(ar(i  ,j  ,kp1)) * frac1i * frac1j * frac0k
     &                 + log(ar(i  ,jp1,kp1)) * frac1i * frac0j * frac0k
     &                 + log(ar(ip1,j  ,kp1)) * frac0i * frac1j * frac0k
     &                 + log(ar(ip1,jp1,kp1)) * frac0i * frac0j * frac0k
             int3dmlog = exp(int3dmlog)
c            print *,'int3dmlog 11: ',rid,rjd,rkd,int3dmlog
           endif
        endif
      endif
      end

      real function int3dmvc(ar,n1,n2,n3,rid,rjd,rkd,misdat)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The interpolation includes the
c        testing of the missing data flag 'misdat'.
c        In the vertical a Lagrangian cubic interpolation is
c        performed.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
 
c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh,klow,n
      real      frac0i,frac0j,frac1i,frac1j,ri,rj,rk
      real      int2d(4)

      real      int3dm

c     if n3 < 4 then do call linear interpolation in the vertical
      if (n3.lt.4) then
        int3dmvc=int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
        return
      endif
 
c     do linear interpolation in the horizontal
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)
 
c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif
 
c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif
 
c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif
 
      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           int3dmvc=ar(i,j,k)
c          print *,'int3dmvc 00: ',rid,rjd,rkd,int3dmvc
        else
c          horizontal interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  ))) then
             int3dmvc=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dmvc = ar(i  ,j  ,k  ) * frac1i * frac1j
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
c            print *,'int3dmvc 10: ',rid,rjd,rkd,int3dmvc
           endif
        endif
      else
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1))) then
             int3dmvc=misdat
           else
             if (k-1.lt.1) then
               klow=1
             else if (k+2.gt.n3) then
               klow=n3-3
             else
               klow=k-1
             endif
             call cubint(ar(i,j,klow),ar(i,j,klow+1),ar(i,j,klow+2),
     &                   ar(i,j,klow+3),klow,rk,int3dmvc)
           endif
        else
c          full 3d interpolation
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
     &         (misdat.eq.ar(ip1,jp1,kp1))) then
             int3dmvc=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             if (k-1.lt.1) then
               klow=1
             else if (k+2.gt.n3) then
               klow=n3-3
             else
               klow=k-1
             endif
             do n=1,4
               int2d(n) = ar(i  ,j  ,klow-1+n  ) * frac1i * frac1j
     &                  + ar(i  ,jp1,klow-1+n  ) * frac1i * frac0j
     &                  + ar(ip1,j  ,klow-1+n  ) * frac0i * frac1j
     &                  + ar(ip1,jp1,klow-1+n  ) * frac0i * frac0j
             enddo
             call cubint(int2d(1),int2d(2),int2d(3),int2d(4),
     &                   klow,rk,int3dmvc)
           endif
        endif
      endif
      end
      real function int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 4d-array to an arbitrary
c        location within the grid.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,..,n4  int   input   dimensions of ar
c        ri,..,rl  real  input   grid location to be interpolated to
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3,n4
      real      ar(n1,n2,n3,n4), rid,rjd,rkd,rld

c     local declarations
      integer   l,lp1,lh
      real      frac0l,frac1l,rl,int3d

c     do linear interpolation in l-direction
      rl=amax1(1.,amin1(float(n4),rld))
      lh=nint(rl)

c     Check for interpolation in l
*     if (abs(float(lh)-rl).lt.1.e-3) then
*       l  =lh
*       lp1=lh
*     else
        l =min0(int(rl),n4-1)
        lp1=l+1
*     endif

      if (l.eq.lp1) then
c       no interpolation in l
        int4d=int3d(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd)
      else
c       interpolation in l
        frac0l=rl-float(l)
        frac1l=1.-frac0l
        int4d = int3d(ar(1,1,1,l  ),n1,n2,n3,rid,rjd,rkd) * frac1l
     &        + int3d(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd) * frac0l
      endif
      end
      real function int4dm(ar,n1,n2,n3,n4,rid,rjd,rkd,rld,misdat)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 4d-array to an arbitrary
c        location within the grid. The interpolation includes the 
c        testing of the missing data flag 'misdat'. 
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,..,n4  int   input   dimensions of ar
c        ri,..,rl  real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested.
c     History:
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3,n4
      real      ar(n1,n2,n3,n4), rid,rjd,rkd,rld, misdat

c     local declarations
      integer   l,lp1,lh
      real      frac0l,frac1l,rl,rint0,rint1,int4d,int3dm

c     check whether missing data checking is required
      if (misdat.eq.0.) then
        int4dm=int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
        return
      endif

c     do linear interpolation in l-direction
      rl=amax1(1.,amin1(float(n4),rld))
      lh=nint(rl)

c     Check for interpolation in l
*     if (abs(float(lh)-rl).lt.1.e-3) then
*       l  =lh
*       lp1=lh
*     else
        l =min0(int(rl),n4-1)
        lp1=l+1
*     endif

      if (l.eq.lp1) then
c       no interpolation in l
        int4dm = int3dm(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd,misdat)
      else
c       interpolation in l
        frac0l=rl-float(l)
        frac1l=1.-frac0l
        rint0 = int3dm(ar(1,1,1,l  ),n1,n2,n3,rid,rjd,rkd,misdat)
        rint1 = int3dm(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd,misdat)
        if ((rint0.eq.misdat).or.(rint1.eq.misdat)) then
          int4dm = misdat
        else
          int4dm = rint0*frac1l + rint1*frac0l
        endif
      endif
      end
      real function int3dl(ar,n1,n2,n3,levels,rid,rjd,rkd)
c-----------------------------------------------------------------------
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The vertical interpolation is linear
c        in log(pressure).
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        levels    real  input   array contains pressure levels for ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c     History:
c        Based on int3d                 July 93
c-----------------------------------------------------------------------

c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd
      real      levels(n3)

c     local declarations
      real      pval
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk

c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)

c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif

c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif

c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif

      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           int3dl=ar(i,j,k)
c          print *,'int3dl 00: ',rid,rjd,rkd,int3dl
        else
c          horizontal interpolation only
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3dl = ar(i  ,j  ,k  ) * frac1i * frac1j
     &            + ar(i  ,jp1,k  ) * frac1i * frac0j
     &            + ar(ip1,j  ,k  ) * frac0i * frac1j
     &            + ar(ip1,jp1,k  ) * frac0i * frac0j
c          print *,'int3dl 10: ',rid,rjd,rkd,int3dl
        endif
      else
*       frac0k=rk-float(k)
c       calculate the pressure value to be interpolated to
        pval=levels(int(rk))
     >            -(rk-aint(rk))*(levels(int(rk))-levels(int(rk)+1))
        frac0k=log(levels(int(rk))/pval)
     &        /log(levels(int(rk))/levels(int(rk)+1))
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           int3dl = ar(i  ,j  ,k  ) * frac1k
     &            + ar(i  ,j  ,kp1) * frac0k
c          print *,'int3dl 01: ',rid,rjd,rkd,int3dl
        else
c          full 3d interpolation
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3dl = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
     &            + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
     &            + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
     &            + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
     &            + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
     &            + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
     &            + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
     &            + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c          print *,'int3dl 11: ',rid,rjd,rkd,int3dl
        endif
      endif
      end
      subroutine bcucof(y,y1,y2,y12,c)
c-----------------------------------------------------------------------
c       Given arrays y,y1,y2 and y12, each of length 4, containing the
c       function, gradients, and cross derivative at the four grid points
c       of a rectangular grid cell (numbered counterclockwise from the 
c       lower left), and given d1 and d2, the length of the grid cell in
c       the 1- and 2-directions, this routine returns the table c that is
c       used by routine bcuint for biqubic interpolation.
c     Source: Numerical Recipes, Fortran Version, p.99
c-----------------------------------------------------------------------
      real      c(4,4),y(4),y1(4),y2(4),y12(4),cl(16),x(16)
      integer   wt(16,16)

      data      wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,
     >          8*0,3,0,-9,6,-2,0,6,-4,
     >          10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,
     >          4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,
     >          10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,
     >          0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,
     >          10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,
     >          5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,
     >          10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/

      real      xx
      integer   i,j,k,l

      do i=1,4                  ! pack a temporary vector x
        x(i)=y(i)
        x(i+4)=y1(i)
        x(i+8)=y2(i)
        x(i+12)=y12(i)
      enddo     

      do i=1,16                 ! matrix multiply by the stored table
        xx=0.
        do k=1,16
          xx=xx+wt(i,k)*x(k)
        enddo
        cl(i)=xx
      enddo

      l=0
      do i=1,4                  ! unpack the result into the output table
      do j=1,4
        l=l+1
        c(i,j)=cl(l)
      enddo
      enddo

      return
      end
      subroutine bcuint(y,y1,y2,y12,x1l,x2l,x1,x2,ansy)
c-----------------------------------------------------------------------
c       Bicubic interpolation within a grid square. Input quantities are
c       y,y1,y2,y12 (as described in bcucof); x1l and x1u, the lower and
c       upper coordinates of the grid square in the 1-direction; x2l and
c       x2u likewise for the 2-direction; and x1,x2, the coordinates of
c       the desired point for the interpolation. The interplated function
c       value is returned as ansy. This routine calls bcucof.
c     Source: Numerical Recipes, Fortran Version, p.99/100
c     !!! changed the proposed code !!!
c-----------------------------------------------------------------------

      real      y(4),y1(4),y2(4),y12(4),c(4,4)
      real      ansy,x1,x2,t,u
      integer   i,x1l,x2l

      call bcucof(y,y1,y2,y12,c)

      t=x1-real(x1l)
      u=x2-real(x2l)

      ansy=0.

      do i=4,1,-1
        ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
      enddo

      return
      end
      subroutine cubint(ya1,ya2,ya3,ya4,k,x,y)
c-----------------------------------------------------------------------
c       Interface routine for SR polint for special case of cubic
c       interpolation in index space, with xa=k,k+1,k+2,k+3
c-----------------------------------------------------------------------
 
      integer   k
      real      ya1,ya2,ya3,ya4,x,y

      integer   n
      real      xa(4),ya(4),dy

      do n=1,4
        xa(1)=real(k)
        xa(2)=real(k+1)
        xa(3)=real(k+2)
        xa(4)=real(k+3)
  
        ya(1)=ya1
        ya(2)=ya2
        ya(3)=ya3
        ya(4)=ya4
      enddo

      call polint(xa,ya,4,x,y,dy)

      return
      end
      subroutine polint(xa,ya,n,x,y,dy)
c-----------------------------------------------------------------------
c       Given arrays xa and ya, each of length n, and given a value x, this
c       routine returns a value y, and an error estimate dy. If P(x) is the
c       polynomial of degree n-1 such that p(xa(i))=ya(i),i=1,...,n, then
c       the returned value y=p(x)
c     Source: Numerical Recipes, Fortran Version, p.82
c-----------------------------------------------------------------------

      integer   nmax,n
      parameter(nmax=10)
      real      xa(n),ya(n),x,y,dy
      integer   i,m,ns
      real      c(nmax),d(nmax),den,dif,dift,ho,hp,w

      ns=1
      dif=abs(x-xa(1))

      do i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
      enddo

      y=ya(ns)
      ns=ns-1
      do m=1,n-1
        do i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
        enddo
        if (2*ns.lt.n-m) then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
      enddo

      return
      end
      subroutine filt2d (a,af,f1,f2,nx,ny,fil,misdat,
     &                                     iperx,ipery,ispol,inpol)
c     =============================================================
c     Apply a conservative diffusion operator onto the 2d field a,
c     with full missing data checking.
c
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
c     af    real   out  filtered array, dimensioned (nx,ny), can be
c                       equivalenced with array a in the calling routine
c     f1    real        workarray, dimensioned (nx+1,ny)
c     f2    real        workarray, dimensioned (nx,ny+1)
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
c                       corresponds to one application of the linear filter.
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
c                       the corresponding value is not available. The
c                       misdat-checking can be switched off with with misdat=0.
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
c
c     Christoph Schaer, 1993
 
c     argument declaration
      integer     nx,ny
      real        a(nx,ny),af(nx,ny),f1(nx+1,ny),f2(nx,ny+1),fil,misdat
      integer     iperx,ipery,inpol,ispol
 
c     local variable declaration
      integer     i,j,is
      real        fh
 
c     compute constant fh
      fh=0.125*fil
 
c     compute fluxes in x-direction
      if (misdat.eq.0.) then
        do j=1,ny
        do i=2,nx
          f1(i,j)=a(i-1,j)-a(i,j)
        enddo
        enddo
      else
        do j=1,ny
        do i=2,nx
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
            f1(i,j)=0.
          else
            f1(i,j)=a(i-1,j)-a(i,j)
          endif
        enddo
        enddo
      endif
      if (iperx.eq.1) then
c       do periodic boundaries in the x-direction
        do j=1,ny
          f1(1,j)=f1(nx,j)
          f1(nx+1,j)=f1(2,j)
        enddo
      else
c       set boundary-fluxes to zero
        do j=1,ny
          f1(1,j)=0.
          f1(nx+1,j)=0.
        enddo
      endif
 
c     compute fluxes in y-direction
      if (misdat.eq.0.) then
        do j=2,ny
        do i=1,nx
          f2(i,j)=a(i,j-1)-a(i,j)
        enddo
        enddo
      else
        do j=2,ny
        do i=1,nx
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
            f2(i,j)=0.
          else
            f2(i,j)=a(i,j-1)-a(i,j)
          endif
        enddo
        enddo
      endif
c     set boundary-fluxes to zero
      do i=1,nx
        f2(i,1)=0.
        f2(i,ny+1)=0.
      enddo
      if (ipery.eq.1) then
c       do periodic boundaries in the x-direction
        do i=1,nx
          f2(i,1)=f2(i,ny)
          f2(i,ny+1)=f2(i,2)
        enddo
      endif
      if (iperx.eq.1) then
        if (ispol.eq.1) then
c         do south-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
          enddo
        endif
        if (inpol.eq.1) then
c         do north-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
          enddo
        endif
      endif
 
c     compute flux-convergence -> filter
      if (misdat.eq.0.) then
        do j=1,ny
        do i=1,nx
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
        enddo
        enddo
      else
        do j=1,ny
        do i=1,nx
          if (a(i,j).eq.misdat) then
            af(i,j)=misdat
          else
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
          endif
        enddo
        enddo
      endif
      end

      subroutine pipo(var3d,p3d,lev,var,nx,ny,nz,mdv,mode)
C     ====================================================

C     Interpolates the 3d variable var3d on the pressure surface
C     defined by lev. The interpolated field is returned as var.
C     p3d denotes the 3d pressure array.
C     mode determines the way of vertical interpolation:
C       mode=0 is for linear interpolation
C       mode=1 is for logarithmic interpolation

      integer   nx,ny,nz,mode
      real      lev,mdv
      real      var3d(nx,ny,nz),p3d(nx,ny,nz),var(nx,ny)

      integer   i,j,k
      real      kind
      real      int3dm

      do i=1,nx
      do j=1,ny

        kind=0.
        do k=1,nz-1
          if ((p3d(i,j,k).ge.lev).and.(p3d(i,j,k+1).le.lev)) then
            kind=float(k)+(p3d(i,j,k)-lev)/
     >                   (p3d(i,j,k)-p3d(i,j,k+1))
            goto 100
          endif
        enddo
 100    continue

        if (kind.eq.0.) then
          var(i,j)=mdv
        else
          var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
        endif

      enddo
      enddo

      return
      end

      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
C     ======================================================

C     Interpolates the 3d variable var3d on the isentropic surface
C     defined by lev. The interpolated field is returned as var.
C     th3d denotes the 3d theta array.
C     mode determines the way of vertical interpolation:
C       mode=0 is for linear interpolation
C       mode=1 is for logarithmic interpolation

      integer   nx,ny,nz,mode
      real      lev,mdv
      real      var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)

      integer   i,j,k
      real      kind
      real      int3dm

      do i=1,nx
      do j=1,ny

        kind=0
        do k=1,nz-1
          if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then
            kind=float(k)+(th3d(i,j,k)-lev)/
     >                   (th3d(i,j,k)-th3d(i,j,k+1))
            goto 100
          endif
        enddo
 100    continue

        if (kind.eq.0) then
          var(i,j)=mdv
        else
          var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
        endif

      enddo
      enddo

      return
      end