C    <@(#)	 hydro.F	4.8 3/28/94>
      subroutine HYDRO( len, km, nx, nwater, nland, extmode,
     &                  Hp,Up,Vp,Tp,Sp,Bp,
     &                  H, U, V, T, S, B,
     &                  Hn,Un,Vn,Tn,Sn,Bn,
     &                  Dx_h, Dy_h, Area_h , Area_u, Zbottom,
     &                  ieast,iwest,inorth,isouth, Bh, Bu,
     &                  F, U_metric, V_metric,taux, tauy, 
     &                  Woc, Qoc, Fso, Qrad_net, 
     &                  D_min, D_max,hmin,hmax,binit, materl, mlk,
     &                  pendepth,penfrac, T_deep,S_deep,
     &                  Entrainment, M_O_Depth, dtime, W_ent, 
     *                  Area_h_r, Area_u_r,Ut,Vt,Wt)
c
c  Computes ocean forcing and performs advection.  Modified to allow
c  use on a Sun system.

      implicit none
      
C==> INPUT PARAMETERS

c==>  Dimension info:

      integer len
      integer km
      integer nx
      integer nwater
      integer nland
      logical extmode

c==> Past State
      real Hp(LEN,KM)
      real Up(LEN,KM)
      real Vp(LEN,KM)
      real Tp(LEN,KM)
      real Sp(LEN,KM)
      real Bp(LEN,KM)
      
c==> Present State
      real H(LEN,KM)
      real U(LEN,KM)
      real V(LEN,KM)
      real T(LEN,KM)
      real S(LEN,KM)
      real B(LEN,KM)
      
c==> Future State
      real Hn(LEN,KM)
      real Un(LEN,KM)
      real Vn(LEN,KM)
      real Tn(LEN,KM)
      real Sn(LEN,KM)
      real Bn(LEN,KM)
      
      
c==> Horizontal Grid
     
      real Dx_h(LEN)
      real Dy_h(LEN)
      real Area_h(LEN)
      real Area_u(LEN)
      real Zbottom(len)
      integer ieast(LEN)
      integer iwest(LEN)
      integer inorth(LEN)
      integer isouth(LEN)
      integer Bh(LEN)
      integer Bu(LEN)
      
c==> Coriolis and Metric terms

      real F(LEN)
      real U_metric(LEN)
      real V_metric(LEN)

c==> Forcing

      real taux(LEN)
      real tauy(LEN)
      real Woc(LEN)
      real Qoc(LEN)
      real Fso(LEN)
      real Qrad_net(LEN)
      
c==> Layer control
      
      real D_min(LEN,KM)
      real D_max(LEN,KM)
      real hmin(KM)
      real hmax(KM)
      real binit(KM)
      integer materl
      integer mlk

c==> Mixed layer parameters and extras
      
      real    pendepth
      real    penfrac
      real    T_deep
      real    S_deep
      real    Entrainment(LEN)
      real    M_O_Depth(LEN)
      real dtime
      
      real W_ent(LEN,KM)

#include "param.h"
c==>  Local Storage
      
      real H_tend(LEN,KM)
      real U_tend(LEN,KM)
      real V_tend(LEN,KM)
      real T_tend(LEN,KM)
      real S_tend(LEN,KM)
      real Area_h_r(LEN)
      real Area_u_r(LEN)
      real Hnew(LEN,KM)
      real Tnew
      real Snew
      real Bnew(LEN,KM)
      real Pt(LEN)
      real Zt(LEN)
      real Ze(LEN)
      integer i,k,m
      
      real BYNCY
      real Ut(len,km),Vt(len,km),Wt(len,km)

c---------------------------------------------------------------------


      DO k=1,KM
       DO i=1,LEN
        H_tend(i,k) = 0.
        U_tend(i,k) = 0.
        V_tend(i,k) = 0.
        T_tend(i,k) = 0.
        S_tend(i,k) = 0.
       ENDDO
      ENDDO


 
c  This part computed in setophys.F and save in a common in ogrid.F.h
c  Passed in here from advec.F
c      do i=1,nwater+nland
c       Area_h_r(i) = Bh(i) / ( Area_h(i) + 1.E-36)
c       Area_u_r(i) = Bu(i) / ( Area_u(i) + 1.E-36)
c      enddo

      call SURF_HEIGHT( len, km, nwater, extmode, H, B, Zbottom, Zt )
 
      if ( extmode ) then
       do i=1,nwater
        Pt(i) = Zt(i) * GRAVITY
        Ze(i) = Zbottom(i)
       enddo
      else ! extmode
       do i=1,len
        Pt(i) = Zt(i) * GRAVITY
        Ze(i) = Zt(i)
        Zt(i) = 0.
       enddo
      endif ! extmode

      
    
       call ZONADV( LEN, KM, nx, nwater, nland,
     &               H, U, V, T, S, 
     &               ieast, iwest, inorth, isouth,
     &               Dy_h, Bh, Bu,
     &               H_tend,U_tend,V_tend,T_tend,S_tend,
     * Ut)

       
       call MERADV(LEN, KM, nx, nwater, nland,
     &               H, U, V, T, S, 
     &               ieast, iwest, inorth, isouth,
     &               Dx_h, Bh, Bu,
     &               H_tend,U_tend,V_tend,T_tend,S_tend,
     * Vt)
 
         

       call CORIOL( LEN, KM, nwater, nland,
     &              Hp,U, V,
     &              Area_h, F, U_metric, V_metric,
     &              ieast, inorth,
     &              U_tend,V_tend)


       call SFLUX(  LEN, KM, nwater, nland,
     &              taux, tauy, Woc,Qoc,Fso,
     &              Area_h, Area_u, 
     &              H_tend , U_tend , V_tend , T_tend, S_tend )

       call PEN_RADIATION(  LEN, KM, nwater,
     &              Hp, Qrad_net, 
     &              Area_h, PENDEPTH, PENFRAC,
     &              T_tend )


      call PRESS( LEN, KM, nwater, nland,
     &               H, B, Pt, Zt,
     &               ieast, iwest, inorth, isouth,
     &               Dx_h, Dy_h,
     &               U_tend,V_tend)


       do k=1,km
       do i=1,nwater
        Hnew(i,k) = Hp(i,k) + H_tend(i,k)* dtime * Area_h_r(i)
        Tnew      = ( Tp(i,k)*Hp(i,k) + T_tend(i,k)*dtime*Area_h_r(i) )
     &               / Hnew(i,k)
        Snew      = ( Sp(i,k)*Hp(i,k) + S_tend(i,k)*dtime*Area_h_r(i) )
     &               / Hnew(i,k)
        Bnew(i,k) = BYNCY( Tnew, Snew, 0. )
       enddo
       do i=nwater+1,LEN
        Hnew(i,k) = 0.
        Bnew(i,k) = 0.
       enddo
       enddo

       call VVELOC( LEN, nwater, KM , extmode, Hnew ,  
     &             Bnew , Binit , W_ent , 
     &             MATERL , MLK , Entrainment, M_O_Depth , Area_h ,
     &             D_min, D_max, hmin , hmax , dtime )
      
       
       call VERADV( LEN, KM, nx, nwater, nland, extmode,
     &             Hp,Up,Vp,Tp,Sp,
     &             H,U,V,T,S,
     &             ieast, inorth, 
     &             W_ent, Area_h,T_deep,S_deep,
     &             H_tend,U_tend,V_tend,T_tend,S_tend,
     * Wt)

      call UPDATE(  LEN, KM, nx, nwater, nland, dtime,
     &              Hp, Up, Vp, Tp, Sp, 
     &              H,  U,  V,  T,  S,  
     &              Hn, Un, Vn, Tn, Sn, 
     &              B, Bn,
     &              ieast, inorth, Bu,
     &              Area_h, Area_h_r, Area_u_r,
     &              H_tend,U_tend,V_tend,T_tend,S_tend)

       
      END
      
     
