c   <@(#)	 filter.F 4.9 4/22/94 PSS>
      subroutine FILTER_H( len, km, im, k1, mlk, extmode,
     &                     H, Area_h, Bh, D_min, D_max, hmin, hmax,
     &                     ieast,iwest,inorth,isouth,
     &                     order, strength_ratio )
      implicit none
      integer len
      integer km
      integer im
      integer k1
      integer mlk
      logical extmode
      real    H( len, km)
      real    Area_h( len )
      integer Bh( len )
      real    D_min(len,km)   
      real    D_max(len,km)
      real    hmin(km)
      real    hmax(km)
      integer ieast( len )
      integer iwest( len )
      integer inorth( len )
      integer isouth( len )
      integer order
      real    strength_ratio

c==> Local Variables
      integer i,k,m,nf
      real mass,tarea
            
      real Wk1(LEN)
      real Z(LEN,KM)
      integer Msk(LEN)
      real SSUM
      real factor
      

CFPP$ EXPAND( BYNCY     ) R
       
c---------------------------------------------------------------------

           

        if (extmode) then
        
         
         ! We rely on D_max(*,km) = D_min(*,km) = -Zbottom(*)
         ! Compute up from fixed bottom, Wk1 will be surface height
         
          do k=1,km
           do i=im+1,len
            Z(i,k) = 0.0
           enddo
          enddo
          do i=1,im
           Z(i,km) = D_max(i,km) 
          enddo
          do k=km,2,-1
           do i=1,im
            Z(i,k-1) = Z(i,k) - H(i,k)
           enddo
          enddo
          do i=1,im
           Wk1(i) = Z(i,1) - H(i,1)
          enddo
          
           
          call SHAPCONS(im,Wk1,Area_h,Bh,
     &       ieast,iwest,inorth,isouth,order,strength_ratio)
          do k=1,km-1
           call SHAPCONS(im,Z(1,k),Area_h,Bh,
     &       ieast,iwest,inorth,isouth,order,strength_ratio)
          enddo


          do i=1,im
            Z(i,1) = AMAX1( Z(i,1), Wk1(i) + hmin(1) )
            Z(i,1) = AMIN1( Z(i,1), Wk1(i) + hmax(1) )
            Z(i,1) = AMIN1( Z(i,1), D_max(i,1) )
            Z(i,1) = AMAX1( Z(i,1), D_min(i,1) )
          enddo
          do k=2,km-1
           do i=1,im
            Z(i,k) = AMAX1( Z(i,k), Z(i,k-1) + hmin(k) )
            Z(i,k) = AMIN1( Z(i,k), Z(i,k-1) + hmax(k) )
            Z(i,k) = AMIN1( Z(i,k), D_max(i,k) )
            Z(i,k) = AMAX1( Z(i,k), D_min(i,k) )
           enddo
          enddo ! k
          
          do i=1,im
           H(i,1) = Z(i,1) - Wk1(i)
          enddo
          do k=2,km
           do i=1,im
            H(i,k) = Z(i,k) - Z(i,k-1)
           enddo
          enddo
          
          
     
        else ! extmode
        
        
         do i=1,im
          Z(i,1) = H(i,1)
         enddo
         do k=2,km
          do i=1,im
            Z(i,k) = Z(i,k-1) + H(i,k)
          enddo
         enddo
        
        do k=k1,km
         
         call SHAPCONS(im,Z(1,k),Area_h,Bh,
     &       ieast,iwest,inorth,isouth,order,strength_ratio)

         do i=1,im
          Z(i,k) = AMIN1( Z(i,k), D_max(i,k) )
          Z(i,k) = AMAX1( Z(i,k), D_min(i,k) )
         enddo
         
         if ( k.EQ.1) then
         
       ! Top layer
         
          do i=1,im
           Z(i,1) = AMIN1( Z(i,1), hmax(1) )
           Z(i,1) = AMAX1( Z(i,1), hmin(1) )
          enddo
          
          
         elseif ( k.eq.k1) then   
           
      ! We are not starting at mixed layer
      ! We are starting at bottom of buffer layer
      ! range for buffer layer is mlk+1 to k1
          
          if ( mlk .eq. 1 ) then
            do i=1,im
              Z(i,k) = AMAX1(Z(i,k),Z(i,1)+SSUM(k1-mlk,hmin(2),1))
              Z(i,k) = AMIN1(Z(i,k),Z(i,1)+SSUM(k1-mlk,hmax(2),1))
            enddo
           else
            do i=1,im
              Z(i,k) = AMAX1( Z(i,k), SSUM( k1, hmin,1) )
              Z(i,k) = AMIN1( Z(i,k), SSUM( k1, hmax,1) )
            enddo
           endif

       ! General Isopycnal Layers
          
         else
          do i=1,im
           Z(i,k) = AMIN1( Z(i,k), Z(i,k-1) + hmax(k) ) 
           Z(i,k) = AMAX1( Z(i,k), Z(i,k-1) + hmin(k) )
          enddo
         endif
         
        enddo ! k
        
        !  Re-evaluate buffer layers
        
        if ( mlk .eq. 0 ) then
          if ( k1 .gt. 1 ) then
            do k=1,k1-1
             factor = FLOAT(k) / k1
             do i=1,im
              Z(i,k) = factor * Z(i,k1)
             enddo
            enddo
           endif ! ( k1.gt.1 )
         else ! (mlk.eq.0)
          if ( k1 .gt. mlk+1 ) then
           do k=mlk+1, k1-1
            factor = ( k - mlk ) / ( k1 - mlk )
            do i=1,im
              Z(i,k) = Z(i,mlk) + factor*( Z(i,k1) - Z(i,mlk) )
            enddo
           enddo
          endif ! (k1.gt.mlk+1)
         endif ! (mlk.eq.0)
         
         
        !  All layer bottoms are now known, re-compute thicknesses 
         
        do i=1,im
         H(i,1) = Z(i,1)
        enddo
        do k=2,km
         do i=1,im
           H(i,k) = Z(i,k) - Z(i,k-1)
         enddo
        enddo
        
       endif ! extmode

      end