subroutine layer2w(u,v,p,depthu,depthv,scpx,scpy,
     &                   ip,iu,iv, kpu,kpv, ii,jj,kk, uflux,vflux,
     &                   ztop,zbot,flag, wz)
      implicit none
c
      integer ii,jj,kk,
     &        ip(ii,jj),iu(ii,jj),iv(ii,jj),kpu(ii,jj),kpv(ii,jj)
      real    u(ii,jj,kk),v(ii,jj,kk),p(ii,jj,kk+1),
     &        depthu(ii,jj),depthv(ii,jj),
     &        scpx(ii,jj),scpy(ii,jj),uflux(ii,jj),vflux(ii,jj),
     &        ztop,zbot,flag, wz(ii,jj)
c
c**********
c*
c  1) return vertical velocity difference between ztop and zbot.
c
c  2) input arguments:
c       u      - u-velocity in layer space
c       v      - v-velocity in layer space
c       p      - layer interface depths (non-negative m)
c                  p(:,:,   1) is the surface
c                  p(:,:,kk+1) is the bathymetry (i.e. depthp)
c       depthu - depth on u-grid
c       depthv - depth on v-grid
c       scpx   - x grid spacing
c       scpy   - y grid spacing
c       ip     - p-grid mask
c       iu     - u-grid mask
c       iv     - v-grid mask
c       kpu    - u-grid layer containing ztop
c       kpv    - v-grid layer containing ztop
c       ii     - 1st array dimension
c       jj     - 2nd array dimension
c       kk     - 3rd array dimension (number of layers)
c       ztop   - top    of z-level cell
c       zbot   - bottom of z-level cell
c       flag   - data void (land) marker
c
c  3) output arguments:
c       kpu    - u-grid layer containing zbot
c       kpv    - v-grid layer containing zbot
c       uflux  - workspace
c       vflux  - workspace
c       wz     - vertical velocity difference 
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c           0 <= ztop <= zbot
c     note that zbot > p(i,j,kk+1) implies that wz(i,j)=flag,
c      since the z-level is then below the bathymetry.
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, February 2002.
c*
c**********
c
      integer i,j,k
      real    puk,pukp,pvk,pvkp,uflx,vflx
c
*     write(6,'(a,2g15.5)')
*    &  'layer2w - ztop,zbot =',ztop,zbot
*     call flush(6)
c
c     u-transport between ztop and zbot.
c
      do j= 1,jj
        do i= 2,ii
          if     (iu(i,j).ne.1 .or. ztop.ge.depthu(i,j)) then
            uflux(i,j) = 0.0
            kpu  (i,j) = kk
            cycle
          endif
          uflx = 0.0
          k    = kpu(i,j)-1
          pukp = min(depthu(i,j),.5*(p(i,j,k+1)+p(i-1,j,k+1)))
          do k= kpu(i,j),kk
            puk  = pukp
            pukp = min(depthu(i,j),.5*(p(i,j,k+1)+p(i-1,j,k+1)))
            if     (k.eq.kpu(i,j)) then
              if     (ztop.lt.puk .or. ztop.gt.pukp) then
                write(6,'(a,i5,a,i5,a,i3, a,f8.2, a,2f8.2,a)')
     &            'error - puk(',i,',',j,') = ',k,
     &            '  but ztop =',ztop,
     &            ' not in that layer (',puk,pukp,')'
                stop
              endif
            endif
            uflx = uflx + u(i,j,k)*(min(zbot,pukp) - max(ztop,puk))
*           if     (i.eq.ii/2 .and. j.eq.jj/2) then
*             write(6,'(a,i3,3g15.5)')
*    &          'layer2w - k,puk,pukp,uflx =',k,puk,pukp,uflx
*             call flush(6)
*           endif
            if     (zbot.ge.puk .and. zbot.le.pukp) then
              kpu(i,j) = k
              exit
            endif
          enddo
          uflux(i,j) = uflx
        enddo
      enddo
c
c     v-transport between ztop and zbot.
c
      do j= 2,jj
        do i= 1,ii
          if     (iv(i,j).ne.1 .or. ztop.ge.depthv(i,j)) then
            vflux(i,j) = 0.0
            kpv  (i,j) = kk
            cycle
          endif
          vflx = 0.0
          k    = kpv(i,j)-1
          pvkp = min(depthv(i,j),.5*(p(i,j,k+1)+p(i,j-1,k+1)))
          do k= kpv(i,j),kk
            pvk  = pvkp
            pvkp = min(depthv(i,j),.5*(p(i,j,k+1)+p(i,j-1,k+1)))
            if     (k.eq.kpv(i,j)) then
              if     (ztop.lt.pvk .or. ztop.gt.pvkp) then
                write(6,'(a,i5,a,i5,a,i3, a,f8.2, a,2f8.2,a)')
     &            'error - pvk(',i,',',j,') = ',k,
     &            '  but ztop =',ztop,
     &            ' not in that layer (',pvk,pvkp,')'
                stop
              endif
            endif
            vflx = vflx + v(i,j,k)*(min(zbot,pvkp) - max(ztop,pvk))
*           if     (i.eq.ii/2 .and. j.eq.jj/2) then
*             write(6,'(a,i3,3g15.5)')
*    &          'layer2w - k,pvk,pvkp,vflx =',k,pvk,pvkp,vflx
*             call flush(6)
*           endif
            if     (zbot.ge.pvk .and. zbot.le.pvkp) then
              kpv(i,j) = k
              exit
            endif
          enddo
          vflux(i,j) = vflx
        enddo
      enddo
c
c     vertical velocity.
c
      do j= 2,jj-1
        wz(1, j) = flag
        do i= 2,ii-1
          if     (ip(i,j).eq.1) then
            if     (ztop.lt.p(i,j,kk+1)) then
              wz(i,j) = (uflux(i+1,j)-uflux(i,j))/scpx(i,j) +
     &                  (vflux(i,j+1)-vflux(i,j))/scpy(i,j)
            else
              wz(i,j) = flag
            endif
          else
            wz(i,j) = flag
          endif
        enddo
        wz(ii,j) = flag
      enddo
      do i= 1,ii
        wz(i, 1) = flag
        wz(i,jj) = flag
      enddo
      return
      end

      subroutine infac2z(a,p,az,z,flag,ii,jj,kk,kz,itype)
      implicit none
c
      integer ii,jj,kk,kz,itype
      real    a(ii,jj,kk+1),p(ii,jj,kk+1),az(ii,jj,kz),z(kz),flag
c
c**********
c*
c  1) interpolate a field at layer interfaces to fixed z depths
c
c  2) input arguments:
c       a     - scalar field in layer interface space
c       p     - layer interface depths (non-negative m)
c                 p(:,:,   1) is the surface
c                 p(:,:,kk+1) is the bathymetry
c       z     - target z-level  depths (non-negative m)
c       flag  - data void (land) marker
c       ii    - 1st dimension of a,p,az
c       jj    - 2nd dimension of a,p,az
c       kk    - 3rd dimension of a  (number of layers)
c       kz    - 3rd dimension of az (number of levels)
c       itype - interpolation type
c                 =1; linear interpolation between interfaces
c
c  3) output arguments:
c       az    - scalar field in z-space
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c           0 <= z(k) <= z(k+1)
c     note that z(k) > p(i,j,kk+1) implies that az(i,j,k)=flag,
c      since the z-level is then below the bathymetry.
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, February 2002.
c*
c**********
c
      integer i,j,k,l,lf
      real    s,zk,z0,zm,zp
c
      do j= 1,jj
        do i= 1,ii
          if     (a(i,j,1).eq.flag) then
            do k= 1,kz
              az(i,j,k) = flag  ! land
            enddo
          else
            lf=1
            do k= 1,kz
              zk=z(k)
              do l= lf,kk
                if     (p(i,j,l).le.zk .and. p(i,j,l+1).ge.zk) then
c
c                 z(k) is in layer l.
c
*                 if     (itype.eq.1) then
c
c                   linear interpolation between interfaces
c
                    z0 = p(i,j,l+1)
                    zm = p(i,j,l)
                    s  = (z0 - zk)/(z0 - zm)
                    az(i,j,k) = s*a(i,j,l) + (1.0-s)*a(i,j,l+1)
*                 endif
                  lf = l
                  exit
                elseif (l.eq.kk) then
                  az(i,j,k) = flag  ! below the bottom
                  lf = l
                  exit
                endif
              enddo !l
            enddo !k
          endif
        enddo !i
      enddo !j
      return
      end

      subroutine layer2z(a,p,az,z,flag,ii,jj,kk,kz,itype)
      implicit none
c
      integer ii,jj,kk,kz,itype
      real    a(ii,jj,kk),p(ii,jj,kk+1),az(ii,jj,kz),z(kz),flag
c
c**********
c*
c  1) interpolate a layered field to fixed z depths
c
c  2) input arguments:
c       a     - scalar field in layer space
c       p     - layer interface depths (non-negative m)
c                 p(:,:,   1) is the surface
c                 p(:,:,kk+1) is the bathymetry
c       z     - target z-level  depths (non-negative m)
c       flag  - data void (land) marker
c       ii    - 1st dimension of a,p,az
c       jj    - 2nd dimension of a,p,az
c       kk    - 3rd dimension of a  (number of layers)
c       kz    - 3rd dimension of az (number of levels)
c       itype - interpolation type
c                 =-2; piecewise quadratic across each layer
c                 =-1; piecewise linear    across each layer
c                 = 0; sample the layer spaning each depth
c                 = 1; linear interpolation between layer centers
c                 = 2; linear interpolation between layer interfaces
c
c  3) output arguments:
c       az    - scalar field in z-space
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c           0 <= z(k) <= z(k+1)
c     note that z(k) > p(i,j,kk+1) implies that az(i,j,k)=flag,
c      since the z-level is then below the bathymetry.
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, February 2002.
c*
c**********
c
      integer i,j,k,l,lf
      real    s,zk,z0,zm,zp
      real    si(kk,1),pi(kk+1),so(kz,1)
c
      do j= 1,jj
        do i= 1,ii
          if     (a(i,j,1).eq.flag) then
            do k= 1,kz
              az(i,j,k) = flag  ! land
            enddo
          elseif (itype.lt.0) then
            do k= 1,kk
              si(k,1) = a(i,j,k)
              pi(k)   = p(i,j,k)
            enddo
            pi(kk+1) = p(i,j,kk+1)
            if     (itype.eq.-1) then
              call layer2z_plm(si,pi,kk,1,so,z,kz, flag)
            else
              call layer2z_ppm(si,pi,kk,1,so,z,kz, flag)
            endif
            do k= 1,kz
              az(i,j,k) = so(k,1)
            enddo
          else !itype.ge.0
            lf=1
            do k= 1,kz
              zk=z(k)
              do l= lf,kk
                if     (p(i,j,l).le.zk .and. p(i,j,l+1).ge.zk) then
c
c                 z(k) is in layer l.
c
                  if     (itype.eq.0) then
c
c                   sample the layer
c
                    az(i,j,k) = a(i,j,l)
                  elseif (itype.eq.1) then
c
c                   linear interpolation between layer centers
c
                    z0 = 0.5*(p(i,j,l)+p(i,j,l+1))
                    if     (zk.le.z0) then
c
c                     z(k) is in the upper half of the layer
c
                      if     (l.eq.1) then
                        az(i,j,k) = a(i,j,1)
                      else
                        zm = 0.5*(p(i,j,l-1)+p(i,j,l))
                        s  = (z0 - zk)/(z0 - zm)
                        az(i,j,k) = s*a(i,j,l-1) + (1.0-s)*a(i,j,l)
                      endif
                    else
c
c                     z(k) is in the lower half of the layer
c
                      if     (p(i,j,l+1).eq.p(i,j,kk+1)) then
                        az(i,j,k) = a(i,j,kk)
                      else
                        zp = 0.5*(p(i,j,l+1)+p(i,j,l+2))
                        s  = (zk - z0)/(zp - z0)
                        az(i,j,k) = s*a(i,j,l+1) + (1.0-s)*a(i,j,l)
                      endif
                    endif
                  else !itype.eq.2
c
c                   linear interpolation between layer interfaces
c                   a is a vertical integral from the surface
c
                    if     (l.eq.1) then
                      s  = zk/p(i,j,2)
                      az(i,j,k) = s*a(i,j,1)
                    else
                      s  = (zk - p(i,j,l))/(p(i,j,l+1) - p(i,j,l))
                      az(i,j,k) = (1.0-s)*a(i,j,l-1) + s*a(i,j,l)
                    endif
                  endif
                  lf = l
                  exit
                elseif (l.eq.kk) then
                  az(i,j,k) = flag  ! below the bottom
                  lf = l
                  exit
                endif
              enddo !l
            enddo !k
          endif !land:itype<0:else
        enddo !i
      enddo !j
      return
      end

      subroutine layer2z_bot(a,p,az,zbot,flag,ii,jj,kk,itype)
      implicit none
c
      integer ii,jj,kk,itype
      real    a(ii,jj,kk),p(ii,jj,kk+1),az(ii,jj),zbot,flag
c
c**********
c*
c  1) interpolate a layered field to a fixed z depth above the bottom
c
c  2) input arguments:
c       a     - scalar field in layer space
c       p     - layer interface depths (non-negative m)
c                 p(:,:,   1) is the surface
c                 p(:,:,kk+1) is the bathymetry
c       zbot  - target height above the bottom (positive m)
c       flag  - data void (land) marker
c       ii    - 1st dimension of a,p,az
c       jj    - 2nd dimension of a,p,az
c       kk    - 3rd dimension of a  (number of layers)
c       itype - interpolation type
c                 =-2; piecewise quadratic across each layer
c                 =-1; piecewise linear    across each layer
c                 =0; sample the layer spaning each depth
c                 =1; linear interpolation between layer centers
c                 =2; linear interpolation between layer interfaces
c
c  3) output arguments:
c       az    - scalar field
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, February 2006.
c*
c**********
c
      integer i,j,k,l
      real    s,zk,z0,zm,zp
      real    si(kk,1),pi(kk+1),so(1,1),zo(1)
c
      do j= 1,jj
        do i= 1,ii
          if     (a(i,j,1).eq.flag) then
            az(i,j) = flag  ! land
          elseif (itype.lt.0) then
            do k= 1,kk
              si(k,1) = a(i,j,k)
              pi(k)   = p(i,j,k)
            enddo
            pi(kk+1) = p(i,j,kk+1)
            zo(1)=max(0.0,p(i,j,kk+1)-zbot)
            if     (itype.eq.-1) then
              call layer2z_plm(si,pi,kk,1,so,zo,1, flag)
            else
              call layer2z_ppm(si,pi,kk,1,so,zo,1, flag)
            endif
            az(i,j) = so(1,1)
          else
            zk=max(0.0,p(i,j,kk+1)-zbot)
            do l= 1,kk
              if     (p(i,j,l).le.zk .and. p(i,j,l+1).ge.zk) then
c
c               z(k) is in layer l.
c
                if     (itype.eq.0) then
c
c                 sample the layer
c
                  az(i,j) = a(i,j,l)
                elseif (itype.eq.1) then
c
c                 linear interpolation between layer centers
c
                  z0 = 0.5*(p(i,j,l)+p(i,j,l+1))
                  if     (zk.le.z0) then
c
c                   z(k) is in the upper half of the layer
c
                    if     (l.eq.1) then
                      az(i,j) = a(i,j,1)
                    else
                      zm = 0.5*(p(i,j,l-1)+p(i,j,l))
                      s  = (z0 - zk)/(z0 - zm)
                      az(i,j) = s*a(i,j,l-1) + (1.0-s)*a(i,j,l)
                    endif
                  else
c
c                   z(k) is in the lower half of the layer
c
                    if     (p(i,j,l+1).eq.p(i,j,kk+1)) then
                      az(i,j) = a(i,j,kk)
                    else
                      zp = 0.5*(p(i,j,l+1)+p(i,j,l+2))
                      s  = (zk - z0)/(zp - z0)
                      az(i,j) = s*a(i,j,l+1) + (1.0-s)*a(i,j,l)
                    endif
                  endif
                else !itype.eq.2
c
c                 linear interpolation between layer interfaces
c                 a is a vertical integral from the surface
c
                  if     (l.eq.1) then
                    s  = zk/p(i,j,2)
                    az(i,j) = s*a(i,j,1)
                  else
                    s  = (zk - p(i,j,l))/(p(i,j,l+1) - p(i,j,l))
                    az(i,j) = (1.0-s)*a(i,j,l-1) + s*a(i,j,l)
                  endif
                endif
                exit
              elseif (l.eq.kk) then
                az(i,j) = flag  ! below the bottom
                exit
              endif
            enddo !l
          endif
        enddo !i
      enddo !j
      return
      end

      subroutine layer2c(a,p,az,zi,flag,ii,jj,kk,kz,itype)
      implicit none
c
      integer ii,jj,kk,kz,itype
      real    a(ii,jj,kk),p(ii,jj,kk+1),az(ii,jj,kz),zi(kz+1),flag
c
c**********
c*
c  1) interpolate a layered field to fixed z-cells.
c
c  2) input arguments:
c       a     - scalar field in layer space
c       p     - layer interface depths (non-negative m)
c                 p(:,:,   1) is the surface
c                 p(:,:,kk+1) is the bathymetry
c       zi    - target z-cell interface depths (non-negative m)
c       flag  - data void (land) marker
c       ii    - 1st dimension of a,p,az
c       jj    - 2nd dimension of a,p,az
c       kk    - 3rd dimension of a  (number of layers)
c       kz    - 3rd dimension of az (number of levels)
c       itype - interpolation type
c                 =0; piecewise constant  across each input cell
c                 =1; piecewise linear    across each input cell
c                 =2; piecewise parabolic across each input cell
c               result is the averaged input profile across each output cell
c
c  3) output arguments:
c       az    - scalar field in z-space
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c           0 <= zi(k) <= zi(k+1)
c     note that zi(k) > p(i,j,kk+1) implies that az(i,j,k)=flag,
c      since the entire cell is then below the bathymetry.
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, August 2005.
c*
c**********
c
      integer i,j,k
      real    si(kk,1),pi(kk+1),so(kz,1),po(kz+1)
c
      do j= 1,jj
        do i= 1,ii
          if     (a(i,j,1).eq.flag) then
            do k= 1,kz
              az(i,j,k) = flag  ! land
            enddo
          else
            do k= 1,kk
              si(k,1) = a(i,j,k)
              pi(k)   = p(i,j,k)
            enddo
            pi(kk+1) = p(i,j,kk+1)
            do k= 1,kz+1
              po(k) = zi(k)
            enddo
            if     (itype.eq.0) then
              call layer2c_pcm(si,pi,kk,1,so,po,kz, flag)
            elseif (itype.eq.1) then
              call layer2c_plm(si,pi,kk,1,so,po,kz, flag)
            else
              call layer2c_ppm(si,pi,kk,1,so,po,kz, flag)
            endif
            do k= 1,kz
              az(i,j,k) = so(k,1)
            enddo
          endif
        enddo !i
      enddo !j
      return
      end

      subroutine layer2c_bot(a,p,az,zbot,flag,ii,jj,kk,itype)
      implicit none
c
      integer ii,jj,kk,itype
      real    a(ii,jj,kk),p(ii,jj,kk+1),az(ii,jj),zbot,flag
c
c**********
c*
c  1) interpolate a layered field to a fixed z-cell above the bottom
c
c  2) input arguments:
c       a     - scalar field in layer space
c       p     - layer interface depths (non-negative m)
c                 p(:,:,   1) is the surface
c                 p(:,:,kk+1) is the bathymetry
c       zbot  - target interface height above the bottom (positive m)
c       flag  - data void (land) marker
c       ii    - 1st dimension of a,p,az
c       jj    - 2nd dimension of a,p,az
c       kk    - 3rd dimension of a  (number of layers)
c       itype - interpolation type
c                 =0; piecewise constant  across each input cell
c                 =1; piecewise linear    across each input cell
c                 =2; piecewise parabolic across each input cell
c               result is the averaged input profile across each output cell
c
c  3) output arguments:
c       az    - scalar field
c
c  4) except at data voids, must have:
c           p(:,:,   1) == zero (surface)
c           p(:,:, l+1) >= p(:,:,l)
c           p(:,:,kk+1) == bathymetry
c
c  5) Alan J. Wallcraft, Naval Research Laboratory, August 2005.
c*
c**********
c
      integer i,j,k
      real    si(kk,1),pi(kk+1),so(1,1),po(2)
c
      do j= 1,jj
        do i= 1,ii
          if     (a(i,j,1).eq.flag) then
            az(i,j) = flag  ! land
          else
            do k= 1,kk
              si(k,1) = a(i,j,k)
              pi(k)   = p(i,j,k)
            enddo
            pi(kk+1) = p(i,j,kk+1)
            po(1) = max( 0.0, p(i,j,kk+1) - zbot )
            po(2) =           p(i,j,kk+1)
            if     (itype.eq.0) then
              call layer2c_pcm(si,pi,kk,1,so,po,1, flag)
            elseif (itype.eq.1) then
              call layer2c_plm(si,pi,kk,1,so,po,1, flag)
            else
              call layer2c_ppm(si,pi,kk,1,so,po,1, flag)
            endif
            az(i,j) = so(1,1)
          endif
        enddo !i
      enddo !j
      return
      end

      subroutine layer2c_pcm(si,pi,ki,ks,
     &                       so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise constant across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry are set to flag.
c
c  5) Alan J. Wallcraft,  Naval Research Laboratory,  Sep. 2002 (Aug. 2005).
c*
c**********
c
      real       thin 
      parameter (thin=1.e-6)  ! minimum layer thickness (no division by 0.0)
c
      integer i,k,l,lf
      real    q,zb,zt,sok(ks)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo !k
      else
        lf=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = min( po(k+1), pi(ki+1) )  !limit for correct cell average
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zt.ge.pi(ki+1)) then
c
c ---       cell below the bottom, set to flag
c
            do i= 1,ks
              so(k,i) = flag
            enddo !i
          elseif (zb-zt.lt.thin) then
c
c ---       thin layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lf).gt.zt) then
              WRITE(6,*) 'bad lf = ',lf
              stop
            endif
            do i= 1,ks
              sok(i) = 0.0
            enddo !i
c           recalibrate lf, usualy exit loop immediately
            do l= lf,ki
              if     (pi(l+1).ge.zt) then
                exit
              elseif (k.eq.ki) then
                exit
              endif
            enddo !l
            lf=l
            do l= lf,ki
              if     (pi(l).gt.zb) then
*               WRITE(6,*) 'l,lf= ',l,lf,l-1
c               the input layer is below the output layer
                lf = l-1
                exit
              elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then
c
c               the input layer is completely inside the output layer
c
                q   = max(pi(l+1)-pi(l),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'L,q = ',l,q
              else
c
c               the input layer is partially inside the output layer
c
                q   = max(min(pi(l+1),zb)-max(pi(l),zt),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'l,q = ',l,q
              endif
            enddo !l
            do i= 1,ks
              so(k,i) = sok(i)
            enddo !i
          endif
        enddo !k
      endif
      return
      end subroutine layer2c_pcm

      subroutine layer2z_plm(si,p,kk,ks,
     &                       sz,z,kz,   flag)
      implicit none
c
      integer kk,ks,kz
      real    si(kk,ks),p(kk+1),
     &        sz(kz,ks),z(kz),flag
c
c**********
c*
c  1) interpolate a set of layered field to fixed z depths.
c     method: piecewise linear across each cell
c
c  2) input arguments:
c       si    - scalar fields in layer space
c       p     - layer interface depths (non-negative m)
c                 p(   1) is the surface
c                 p(kk+1) is the bathymetry
c       kk    - dimension of si (number of layers)
c       ks    - dimension of si (number of fields)
c       z     - target z-level  depths (non-negative m)
c       flag  - data void (land) marker
c       kz    - dimension of sz (number of levels)
c
c  3) output arguments:
c       sz    - scalar fields in z-space
c
c  4) except at data voids, must have:
c           p(   1) == zero (surface)
c           p( l+1) >= p(:,:,l)
c           p(kk+1) == bathymetry
c           0 <= z(k) <= z(k+1)
c     note that z(k) > p(kk+1) implies that az(k)=flag,
c      since the z-level is then below the bathymetry.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
c     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real, parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lf
      real    q,zk
      real    sis(kk,ks),pt(kk+1)
c
      if     (si(1,ks).eq.flag) then
        do k= 1,kz
          do i= 1,ks
            sz(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PLM slopes for input layers
        do k=1,kk
          pt(k)=max(p(k+1)-p(k),thin)
        enddo
        call plm(pt,si,sis,kk,ks)
c
        lf=1
        do k= 1,kz
          zk=z(k)
          do l= lf,kk
            if     (p(l).le.zk .and. p(l+1).ge.zk) then
c
c             z(k) is in layer l, sample the linear profile at zk.
c
              q = (zk-p(l))/pt(l) - 0.5
              do i= 1,ks
                sz(k,i) = si(l,i) + q*sis(l,i)
              enddo !i
              lf = l  ! z monotonic increasing, so z(k+1) in layers l:kk
              exit
            elseif (l.eq.kk) then
              do i= 1,ks
                sz(k,i) = flag  ! below the bottom
              enddo !i
              lf = l
              exit
            endif
          enddo !l
        enddo !k
      endif
      return
      end

      subroutine layer2c_plm(si,pi,ki,ks,
     &                       so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise linear across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry are set to flag.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
c     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real,parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lf
      real    q,qc,zb,zc,zt,sok(ks)
      real    sis(ki,ks),pit(ki+1)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PLM slopes for input layers
        do k=1,ki
          pit(k)=max(pi(k+1)-pi(k),thin)
        enddo
        call plm(pit,si,sis,ki,ks)
c ---   compute output layer averages
        lf=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = min( po(k+1), pi(ki+1) )  !limit for correct cell average
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zt.ge.pi(ki+1)) then
c
c ---       cell below the bottom, set to flag
c
            do i= 1,ks
              so(k,i) = flag
            enddo !i
          elseif (zb-zt.lt.thin) then
c
c ---       thin layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lf).gt.zt) then
              WRITE(6,*) 'bad lf = ',lf
              stop
            endif
            do i= 1,ks
              sok(i) = 0.0
            enddo !i
c           recalibrate lf, usualy exit loop immediately
            do l= lf,ki
              if     (pi(l+1).ge.zt) then
                exit
              elseif (k.eq.ki) then
                exit
              endif
            enddo !l
            lf=l
            do l= lf,ki
              if     (pi(l).gt.zb) then
*               WRITE(6,*) 'l,lf= ',l,lf,l-1
c               the input layer is below the output layer
                lf = l-1
                exit
              elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then
c
c               the input layer is completely inside the output layer
c
                q   = max(pi(l+1)-pi(l),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'L,q = ',l,q
              else
c
c               the input layer is partially inside the output layer
c               average of linear profile is its center value
c
                q   = max( min(pi(l+1),zb)-max(pi(l),zt), thin )/(zb-zt)
                zc  = 0.5*(min(pi(l+1),zb)+max(pi(l),zt))
                qc  = (zc-pi(l))/pit(l) - 0.5
                do i= 1,ks
                  sok(i) = sok(i) + q*(si(l,i) + qc*sis(l,i))
                enddo !i
*               WRITE(6,*) 'l,q,qc = ',l,q,qc
              endif
            enddo !l
            do i= 1,ks
              so(k,i) = sok(i)
            enddo !i
          endif
        enddo !k
      endif
      return
      end subroutine layer2c_plm

      subroutine plm(pt, s,ss,ki,ks)
      implicit none
c
      integer ki,ks
      real    pt(ki+1),s(ki,ks),ss(ki,ks)
c
c**********
c*
c  1) generate a monotonic PLM interpolation of a layered field
c
c  2) input arguments:
c       pt    - layer interface thicknesses (non-zero)
c       s     - scalar fields in layer space
c       ki    - 1st dimension of s (number of layers)
c       ks    - 2nd dimension of s (number of fields)
c
c  3) output arguments:
c       ss    - scalar field slopes for PLM interpolation
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(:,:,l)
c           pi(ki+1) == bathymetry
c
c  5) Tim Campbell, Mississippi State University, September 2002.
c*
c**********
c
      integer l
      real    ql(ki),qc(ki),qr(ki)
c
      !compute grid spacing ratios for slope computations
      ql(1)=0.0
      qc(1)=0.0
      qr(1)=0.0
      do l=2,ki-1
        ql(l)=2.0*pt(l)/(pt(l-1)+pt(l))
        qc(l)=2.0*pt(l)/(pt(l-1)+2.0*pt(l)+pt(l+1))
        qr(l)=2.0*pt(l)/(pt(l)+pt(l+1))
      enddo
      ql(ki)=0.0
      qc(ki)=0.0
      qr(ki)=0.0
      !compute normalized layer slopes
      do l=1,ks
        call slope(ql,qc,qr,s(1,l),ss(1,l),ki)
      enddo
      return
      end subroutine plm

      subroutine slope(rl,rc,rr,a,s,n)
      implicit none
c
      integer,intent(in)  :: n
      real,   intent(in)  :: rl(n),rc(n),rr(n),a(n)
      real,   intent(out) :: s(n)
c
c**********
c*
c  1) generate slopes for monotonic piecewise linear distribution
c
c  2) input arguments:
c       rl   - left grid spacing ratio
c       rc   - center grid spacing ratio
c       rr   - right grid spacing ratio
c       a    - scalar field zone averages
c       n    - number of zones
c
c  3) output arguments:
c       s    - zone slopes
c
c  4) Tim Campbell, Mississippi State University, September 2002.
c*
c**********
c
      integer,parameter :: ic=2, im=1, imax=100
      real,parameter :: fracmin=1e-6, dfac=0.5
c
      integer i,j
      real    sl,sc,sr
      real    dnp,dnn,dl,dr,ds,frac
c
c Compute zone slopes
c Campbell Eq(15) -- nonuniform grid
c
      s(1)=0.0
      do j=2,n-1
        sl=rl(j)*(a(j)-a(j-1))
        sr=rr(j)*(a(j+1)-a(j))
        if (sl*sr.gt.0.) then
          s(j)=sign(min(abs(sl),abs(sr)),sl)
        else
          s(j)=0.0
        endif
      enddo
      s(n)=0.0
c
c Minimize discontinuities between zones
c Apply single pass discontinuity minimization: Campbell Eq(19)
c
      do j=2,n-1
        if(s(j).ne.0.0) then
          dl=-0.5*(s(j)+s(j-1))+a(j)-a(j-1)
          dr=-0.5*(s(j+1)+s(j))+a(j+1)-a(j)
          ds=sign(min(abs(dl),abs(dr)),dl)
          s(j)=s(j)+2.0*ds
        endif
      enddo
      return
      end subroutine slope

      subroutine layer2z_ppm(si,p,kk,ks,
     &                       sz,z,kz,   flag)
      implicit none
c
      integer kk,ks,kz
      real    si(kk,ks),p(kk+1),
     &        sz(kz,ks),z(kz),flag
c
c**********
c*
c  1) interpolate a set of layered field to fixed z depths.
c     method: piecewise parabolic method across each cell
c
c  2) input arguments:
c       si    - scalar fields in layer space
c       p     - layer interface depths (non-negative m)
c                 p(   1) is the surface
c                 p(kk+1) is the bathymetry
c       kk    - dimension of si (number of layers)
c       ks    - dimension of si (number of fields)
c       z     - target z-level  depths (non-negative m)
c       flag  - data void (land) marker
c       kz    - dimension of az (number of levels)
c
c  3) output arguments:
c       sz    - scalar fields in z-space
c
c  4) except at data voids, must have:
c           p(   1) == zero (surface)
c           p( l+1) >= p(:,:,l)
c           p(kk+1) == bathymetry
c           0 <= z(k) <= z(k+1)
c     note that z(k) > p(kk+1) implies that az(k)=flag,
c      since the z-level is then below the bathymetry.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
c     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real, parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lf
      real    q,zk
      real    sic(kk,ks,3),pt(kk+1)
c
      if     (si(1,ks).eq.flag) then
        do k= 1,kz
          do i= 1,ks
            sz(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PPM coefficients for input layers
        do k=1,kk
          pt(k)=max(p(k+1)-p(k),thin)
        enddo
        call ppm(pt,si,sic,kk,ks)
c
        lf=1
        do k= 1,kz
          zk=z(k)
          do l= lf,kk
            if     (p(l).le.zk .and. p(l+1).ge.zk) then
c
c             z(k) is in layer l, sample the quadratic profile at zk.
c
              q = (zk-p(l))/pt(l)
              do i= 1,ks
                sz(k,i) =    sic(l,i,1) +
     &                    q*(sic(l,i,2) +
     &                       sic(l,i,3)*(1.0-q))
              enddo !i
              lf = l  ! z monotonic increasing, so z(k+1) in layers l:kk
              exit
            elseif (l.eq.kk) then
              do i= 1,ks
                sz(k,i) = flag  ! below the bottom
              enddo !i
              lf = l
              exit
            endif
          enddo !l
        enddo !k
      endif
      return
      end

      subroutine layer2c_ppm(si,pi,ki,ks,
     &                       so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise parabolic method across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry are set to flag.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
C     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real,parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lb,lt
      real    sz,xb,xt,zb,zt
      real    sic(ki,ks,3),pit(ki+1)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PPM coefficients for input layers
        do k=1,ki
          pit(k)=max(pi(k+1)-pi(k),thin)
        enddo
        call ppm(pit,si,sic,ki,ks)
c ---   compute output layer averages
        lb=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = min( po(k+1), pi(ki+1) )  !limit for correct cell average
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zt.ge.pi(ki+1)) then
c
c ---       cell below the bottom, set to flag
c
            do i= 1,ks
              so(k,i) = flag
            enddo !i
          elseif (zb-zt.lt.thin) then
c
c ---       thin layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lb).gt.zt) then
              WRITE(6,*) 'bad lb = ',lb
              stop
            endif
            lt=lb !top will always correspond to bottom of previous
            lb=lt !find input layer containing bottom output interface
            do while (pi(lb+1).lt.zb.and.lb.lt.ki)
              lb=lb+1
            enddo
            xt=(zt-pi(lt))/pit(lt)
            xb=(zb-pi(lb))/pit(lb)
            do i= 1,ks
              if     (lt.ne.lb) then
                sz=   pit(lt)*(      sic(lt,i,1)  *(1.-xt)
     &                         +0.5*(sic(lt,i,2)+
     &                               sic(lt,i,3) )*(1.-xt**2)
     &                              -sic(lt,i,3)  *(1.-xt**3)/3.0)
                do l=lt+1,lb-1
                  sz=sz+pit(l)*si(l,i)
                enddo !l
                sz=sz+pit(lb)*(      sic(lb,i,1)  *    xb
     &                         +0.5*(sic(lb,i,2)+
     &                               sic(lb,i,3) )*    xb**2
     &                              -sic(lb,i,3)  *    xb**3 /3.0)
              else
                sz=pit(lt)*(      sic(lt,i,1)  *(xb-xt)
     &                      +0.5*(sic(lt,i,2)+
     &                            sic(lt,i,3) )*(xb**2-xt**2)
     &                           -sic(lt,i,3)  *(xb**3-xt**3)/3.0)
              endif
              so(k,i) = sz/(zb-zt)
            enddo !i
          endif !thin:std layer
        enddo !k
      endif
      return
      end subroutine layer2c_ppm

      subroutine ppm(pt, s,sc,ki,ks)
      implicit none
c
      integer ki,ks
      real    pt(ki+1),s(ki,ks),sc(ki,ks,3)
c
c**********
c*
c  1) generate a monotonic PPM interpolation of a layered field:
c     Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201.
c
c  2) input arguments:
c       pt    - layer interface thicknesses (non-zero)
c       s     - scalar fields in layer space
c       ki    - 1st dimension of s (number of layers)
c       ks    - 2nd dimension of s (number of fields)
c
c  3) output arguments:
c       sc    - scalar field coefficients for PPM interpolation
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(:,:,l)
c           pi(ki+1) == bathymetry
c
c  5) Tim Campbell, Mississippi State University, September 2002;
C     Alan J. Wallcraft,  Naval Research Laboratory,  August 2007.
c*
c**********
c
      integer j,k,l
      real    da,a6,slj,scj,srj
      real    as(ki),al(ki),ar(ki)
      real     ptjp(ki), pt2jp(ki), ptj2p(ki),
     &        qptjp(ki),qpt2jp(ki),qptj2p(ki),ptq3(ki),qpt4(ki)
c
      !compute grid metrics
      do j=1,ki
         ptjp( j) = pt(j)   + pt(j+1)
         pt2jp(j) = pt(j)   + ptjp(j)
         ptj2p(j) = ptjp(j) + pt(j+1)
        qptjp( j) = 1.0/ptjp( j)
        qpt2jp(j) = 1.0/pt2jp(j)
        qptj2p(j) = 1.0/ptj2p(j)
      enddo !j
         ptq3(2) = pt(2)/(pt(1)+ptjp(2))
      do j=3,ki-1
         ptq3(j) = pt(j)/(pt(j-1)+ptjp(j))  !pt(j)/      (pt(j-1)+pt(j)+pt(j+1))
         qpt4(j) = 1.0/(ptjp(j-2)+ptjp(j))  !1.0/(pt(j-2)+pt(j-1)+pt(j)+pt(j+1))
      enddo !j
c
      do l= 1,ks
        !Compute average slopes: Colella, Eq. (1.8)
        as(1)=0.
        do j=2,ki-1
          slj=s(j,  l)-s(j-1,l)
          srj=s(j+1,l)-s(j,  l)
          if (slj*srj.gt.0.) then
            scj=ptq3(j)*( pt2jp(j-1)*srj*qptjp(j)
     &                   +ptj2p(j)  *slj*qptjp(j-1) )
            as(j)=sign(min(abs(2.0*slj),abs(scj),abs(2.0*srj)),scj)
          else
            as(j)=0.
          endif
        enddo !j
        as(ki)=0.
        !Compute "first guess" edge values: Colella, Eq. (1.6)
        al(1)=s(1,l)  !1st layer PCM
        ar(1)=s(1,l)  !1st layer PCM
        al(2)=s(1,l)  !1st layer PCM
        do j=3,ki-1
          al(j)=s(j-1,l)+pt(j-1)*(s(j,l)-s(j-1,l))*qptjp(j-1)
     &         +qpt4(j)*(
     &            2.*pt(j)*pt(j-1)*qptjp(j-1)*(s(j,l)-s(j-1,l))*
     &            ( ptjp(j-2)*qpt2jp(j-1)
     &             -ptjp(j)  *qptj2p(j-1) )
     &            -pt(j-1)*as(j)  *ptjp(j-2)*qpt2jp(j-1)
     &            +pt(j)  *as(j-1)*ptjp(j)  *qptj2p(j-1)
     &              )
          ar(j-1)=al(j)
        enddo !j
        ar(ki-1)=s(ki,l)  !last layer PCM
        al(ki)  =s(ki,l)  !last layer PCM
        ar(ki)  =s(ki,l)  !last layer PCM
        !Impose monotonicity: Colella, Eq. (1.10)
        do j=2,ki-1
          if ((s(j+1,l)-s(j,l))*(s(j,l)-s(j-1,l)).le.0.) then !local extremum
            al(j)=s(j,l)
            ar(j)=s(j,l)
          else
            da=ar(j)-al(j)
            a6=6.0*s(j,l)-3.0*(al(j)+ar(j))
            if     (da*a6 .gt.  da*da) then !peak in right half of zone
              al(j)=3.0*s(j,l)-2.0*ar(j)
            elseif (da*a6 .lt. -da*da) then !peak in left half of zone
              ar(j)=3.0*s(j,l)-2.0*al(j)
            endif
          endif
        enddo !j
        !Set coefficients
        do j=1,ki
          sc(j,l,1)=al(j)
          sc(j,l,2)=ar(j)-al(j)
          sc(j,l,3)=6.0*s(j,l)-3.0*(al(j)+ar(j))
        enddo !j
      enddo !l
      return
      end subroutine ppm