C NCLFORTSTART
      subroutine watermassvolume(ns,siglo,sighi,nx,ny,nz,potdens,missval
     +                   ,tvol,tarea,out1,out2)
      implicit none
      integer  ns,nx,ny,nz
      real    siglo(ns), sighi(ns)
      real*8  potdens(nx,ny,nz), tvol(nx,ny,nz), tarea(nx,ny), missval
      real*8  out1(ns)
      real*8  out2(ns)
C NCLEND
      integer  ix,iy,iz,isig,icnt1,icnt2
      real*8   areasum

c initialize out1 (total volume in this density class) and
c            out2 (mean thickness of this density class)
      do isig=1,ns
        out1(isig) = 0.0d0      
        out2(isig) = 0.0d0 
      end do

      do isig=1,ns
       icnt1 = 0
       areasum = 0.
       do iy=1,ny
        do ix=1,nx
         icnt2 = 0
         do iz=1,nz
           if (potdens(ix,iy,iz).ge.siglo(isig) .and.
     +         potdens(ix,iy,iz).lt.sighi(isig)) then
               icnt1=icnt1+1
               icnt2=icnt2+1
               out1(isig) = out1(isig) + tvol(ix,iy,iz)
           end if
         end do
         if (icnt2.gt.0) then
             areasum = areasum + tarea(ix,iy)
         end if
        end do
       end do
         if (icnt1.gt.0) then
             out2(isig) = out1(isig)/areasum
         end if
      end do

      return
      end

C NCLFORTSTART
      subroutine watermassvolume2(ns,siglo,sighi,nx,ny,nz,nt,potdens
     +                   ,missval,tvol,tarea,tlon,tlat,ht,x0,x1,y0,y1
     +                   ,hmin,out1,out2)
      implicit none
      integer  ns,nx,ny,nz,nt
      real    siglo(ns), sighi(ns),x0,x1,y0,y1,hmin
      real*8  potdens(nx,ny,nz,nt), tvol(nx,ny), tarea(nx,ny)
      real*8  tlon(nx,ny), tlat(nx,ny), ht(nx,ny)
      real*8  out1(ns,nt), missval
      real*8  out2(ns,nt)
C NCLEND
      integer  ix,iy,iz,it,isig,icnt1,icnt2
      real*8   areasum

c initialize out1 (total volume in this density class) and
c            out2 (mean thickness of this density class)
      do it=1,nt
      do isig=1,ns
        out1(isig,it) = 0.0d0      
        out2(isig,it) = 0.0d0 
      end do
      end do

      do it=1,nt
      do isig=1,ns
       icnt1 = 0
       areasum = 0.
       do iy=1,ny
        do ix=1,nx
         if (tlon(ix,iy).ge.x0.and.tlon(ix,iy).le.x1.and.
     +       tlat(ix,iy).ge.y0.and.tlat(ix,iy).le.y1.and.
     +       ht(ix,iy).gt.hmin) then
           icnt2 = 0
           do iz=1,nz
             if (potdens(ix,iy,iz,it).ge.siglo(isig) .and.
     +           potdens(ix,iy,iz,it).lt.sighi(isig)) then
                 icnt1=icnt1+1
                 icnt2=icnt2+1
                 out1(isig,it) = out1(isig,it) + tvol(ix,iy)
             end if
           end do
           if (icnt2.gt.0) then
               areasum = areasum + tarea(ix,iy)
           end if
         end if
        end do
       end do
         if (icnt1.gt.0) then
             out2(isig,it) = out1(isig,it)/areasum
         end if
      end do
      end do

      return
      end