pkg/mom_fluxform/mom_fluxform.F

C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.7 2002/11/07 21:51:15 adcroft Exp $
C $Name:  $

CBOI
C !TITLE: pkg/mom\_advdiff
C !AUTHORS: adcroft@mit.edu
C !INTRODUCTION: Flux-form Momentum Equations Package
C
C Package "mom\_fluxform" provides methods for calculating explicit terms
C in the momentum equation cast in flux-form:
C \begin{eqnarray*}
C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
C           -\nabla \cdot {\bf v} u
C           -fv
C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
C           + \mbox{metrics}
C \\
C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
C           -\nabla \cdot {\bf v} v
C           +fu
C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
C           + \mbox{metrics}
C \end{eqnarray*}
C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
C stresses as well as internal viscous stresses.
CEOI

#include "CPP_OPTIONS.h"

CBOP
C !ROUTINE: MOM_FLUXFORM

C !INTERFACE: ==========================================================
      SUBROUTINE MOM_FLUXFORM( 
     I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
     I        phi_hyd,KappaRU,KappaRV,
     U        fVerU, fVerV,
     I        myTime,myIter,myThid)

C !DESCRIPTION:
C Calculates all the horizontal accelerations except for the implicit surface
C pressure gradient and implciit vertical viscosity.

C !USES: ===============================================================
C     == Global variables ==
      IMPLICIT NONE
#include "SIZE.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  iMin,iMax,jMin,jMAx  :: loop ranges
C  k                    :: vertical level
C  kUp                  :: =1 or 2 for consecutive k
C  kDown                :: =2 or 1 for consecutive k
C  phi_hyd              :: hydrostatic pressure (perturbation)
C  KappaRU              :: vertical viscosity
C  KappaRV              :: vertical viscosity
C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining
C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining
C  myTime               :: current time
C  myIter               :: current time-step number
C  myThid               :: thread number
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k,kUp,kDown
      _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C None - updates gU() and gV() in common blocks

C !LOCAL VARIABLES: ====================================================
C  i,j                  :: loop indices
C  aF                   :: advective flux
C  vF                   :: viscous flux
C  v4F                  :: bi-harmonic viscous flux
C  vrF                  :: vertical viscous flux
C  cF                   :: Coriolis acceleration
C  mT                   :: Metric terms
C  pF                   :: Pressure gradient
C  fZon                 :: zonal fluxes
C  fMer                 :: meridional fluxes
      INTEGER i,j
      _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C     wMaskOverride - Land sea flag override for top layer.
C     afFacMom      - Tracer parameters for turning terms
C     vfFacMom        on and off.
C     pfFacMom        afFacMom - Advective terms 
C     cfFacMom        vfFacMom - Eddy viscosity terms
C     mTFacMom        pfFacMom - Pressure terms
C                     cfFacMom - Coriolis terms
C                     foFacMom - Forcing
C                     mTFacMom - Metric term
C     uDudxFac, AhDudxFac, etc ... individual term tracer parameters
      _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS      xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS      yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C     I,J,K - Loop counters
C     rVelMaskOverride - Factor for imposing special surface boundary conditions
C                        ( set according to free-surface condition ).
C     hFacROpen        - Lopped cell factos used tohold fraction of open
C     hFacRClosed        and closed cell wall.
      _RL  rVelMaskOverride
C     xxxFac - On-off tracer parameters used for switching terms off.
      _RL  uDudxFac
      _RL  AhDudxFac
      _RL  A4DuxxdxFac
      _RL  vDudyFac
      _RL  AhDudyFac
      _RL  A4DuyydyFac
      _RL  rVelDudrFac
      _RL  ArDudrFac
      _RL  fuFac
      _RL  phxFac
      _RL  mtFacU
      _RL  uDvdxFac
      _RL  AhDvdxFac
      _RL  A4DvxxdxFac
      _RL  vDvdyFac
      _RL  AhDvdyFac
      _RL  A4DvyydyFac
      _RL  rVelDvdrFac
      _RL  ArDvdrFac
      _RL  fvFac
      _RL  phyFac
      _RL  vForcFac
      _RL  mtFacV
      INTEGER km1,kp1
      _RL wVelBottomOverride
      LOGICAL bottomDragTerms
      _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

      km1=MAX(1,k-1)
      kp1=MIN(Nr,k+1)
      rVelMaskOverride=1.
      IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac
      wVelBottomOverride=1.
      IF (k.EQ.Nr) wVelBottomOverride=0.

C     Initialise intermediate terms
      DO J=1-OLy,sNy+OLy
       DO I=1-OLx,sNx+OLx
        aF(i,j)   = 0.
        vF(i,j)   = 0.
        v4F(i,j)  = 0.
        vrF(i,j)  = 0.
        cF(i,j)   = 0.
        mT(i,j)   = 0.
        pF(i,j)   = 0.
        fZon(i,j) = 0.
        fMer(i,j) = 0.
        rTransU(i,j) = 0.
        rTransV(i,j) = 0.
       ENDDO
      ENDDO

C--   Term by term tracer parmeters
C     o U momentum equation
      uDudxFac     = afFacMom*1.
      AhDudxFac    = vfFacMom*1.
      A4DuxxdxFac  = vfFacMom*1.
      vDudyFac     = afFacMom*1.
      AhDudyFac    = vfFacMom*1.
      A4DuyydyFac  = vfFacMom*1.
      rVelDudrFac  = afFacMom*1.
      ArDudrFac    = vfFacMom*1.
      mTFacU       = mtFacMom*1.
      fuFac        = cfFacMom*1.
      phxFac       = pfFacMom*1.
C     o V momentum equation
      uDvdxFac     = afFacMom*1.
      AhDvdxFac    = vfFacMom*1.
      A4DvxxdxFac  = vfFacMom*1.
      vDvdyFac     = afFacMom*1.
      AhDvdyFac    = vfFacMom*1.
      A4DvyydyFac  = vfFacMom*1.
      rVelDvdrFac  = afFacMom*1.
      ArDvdrFac    = vfFacMom*1.
      mTFacV       = mtFacMom*1.
      fvFac        = cfFacMom*1.
      phyFac       = pfFacMom*1.
      vForcFac     = foFacMom*1.

      IF (     no_slip_bottom
     &    .OR. bottomDragQuadratic.NE.0.
     &    .OR. bottomDragLinear.NE.0.) THEN
       bottomDragTerms=.TRUE.
      ELSE
       bottomDragTerms=.FALSE.
      ENDIF

C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP
      IF (staggerTimeStep) THEN
        phxFac = 0.
        phyFac = 0.
      ENDIF

C--   Calculate open water fraction at vorticity points
      CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)

C---- Calculate common quantities used in both U and V equations
C     Calculate tracer cell face open areas
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        xA(i,j) = _dyG(i,j,bi,bj)
     &   *drF(k)*_hFacW(i,j,k,bi,bj)
        yA(i,j) = _dxG(i,j,bi,bj)
     &   *drF(k)*_hFacS(i,j,k,bi,bj)
       ENDDO
      ENDDO

C     Make local copies of horizontal flow field
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        uFld(i,j) = uVel(i,j,k,bi,bj)
        vFld(i,j) = vVel(i,j,k,bi,bj)
       ENDDO
      ENDDO

C     Calculate velocity field "volume transports" through tracer cell faces.
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        uTrans(i,j) = uFld(i,j)*xA(i,j)
        vTrans(i,j) = vFld(i,j)*yA(i,j)
       ENDDO
      ENDDO

      CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)

C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
      IF (momAdvection.AND.k.EQ.1) THEN

C-    Calculate vertical transports above U & V points (West & South face):
       CALL MOM_CALC_RTRANS( k, bi, bj,
     O                       rTransU, rTransV,
     I                       myTime, myIter, myThid)

C-    Free surface correction term (flux at k=1)
       CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         fVerU(i,j,kUp) = af(i,j)
        ENDDO
       ENDDO

       CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         fVerV(i,j,kUp) = af(i,j)
        ENDDO
       ENDDO

C---  endif momAdvection & k=1
      ENDIF


C---  Calculate vertical transports (at k+1) below U & V points :
      IF (momAdvection) THEN
       CALL MOM_CALC_RTRANS( k+1, bi, bj,
     O                       rTransU, rTransV,
     I                       myTime, myIter, myThid)
      ENDIF


C---- Zonal momentum equation starts here

C     Bi-harmonic term del^2 U -> v4F
      IF (momViscosity)
     & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)

C---  Calculate mean and eddy fluxes between cells for zonal flow.

C--   Zonal flux (fZon is at east face of "u" cell)

C     Mean flow component of zonal flux -> aF
      IF (momAdvection)
     & CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,aF,myThid)

C     Laplacian and bi-harmonic terms -> vF
      IF (momViscosity)
     & CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,vF,myThid)

C     Combine fluxes -> fZon
      DO j=jMin,jMax
       DO i=iMin,iMax
        fZon(i,j) = uDudxFac*aF(i,j) + AhDudxFac*vF(i,j)
       ENDDO
      ENDDO

C--   Meridional flux (fMer is at south face of "u" cell)

C     Mean flow component of meridional flux
      IF (momAdvection)
     & CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,aF,myThid)

C     Laplacian and bi-harmonic term
      IF (momViscosity)
     & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)

C     Combine fluxes -> fMer
      DO j=jMin,jMax
       DO i=iMin,iMax
        fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)
       ENDDO
      ENDDO

C--   Vertical flux (fVer is at upper face of "u" cell)

C     Mean flow component of vertical flux (at k+1) -> aF
      IF (momAdvection)
     & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,rTransU,af,myThid)

C     Eddy component of vertical flux (interior component only) -> vrF
      IF (momViscosity.AND..NOT.implicitViscosity)
     & CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)

C     Combine fluxes
      DO j=jMin,jMax
       DO i=iMin,iMax
        fVerU(i,j,kDown) = rVelDudrFac*aF(i,j) + ArDudrFac*vrF(i,j)
       ENDDO
      ENDDO

C---  Hydrostatic term ( -1/rhoConst . dphi/dx )
      IF (momPressureForcing) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         pf(i,j) = - _recip_dxC(i,j,bi,bj)
     &    *(phi_hyd(i,j,k)-phi_hyd(i-1,j,k))
        ENDDO
       ENDDO
      ENDIF

C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
      DO j=jMin,jMax
       DO i=iMin,iMax
        gU(i,j,k,bi,bj) =
#ifdef OLD_UV_GEOM
     &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
     &    ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
#else
     &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
     &   *recip_rAw(i,j,bi,bj)
#endif
     &  *(fZon(i,j  )          - fZon(i-1,j)
     &   +fMer(i,j+1)          - fMer(i  ,j)
     &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
     &   )
     & _PHM( +phxFac * pf(i,j) )
       ENDDO
      ENDDO

#ifdef NONLIN_FRSURF
C-- account for 3.D divergence of the flow in rStar coordinate:
      IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
     &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
     &       *uVel(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
      IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
     &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
#endif /* NONLIN_FRSURF */

C-- No-slip and drag BCs appear as body forces in cell abutting topography 
      IF (momViscosity.AND.no_slip_sides) THEN
C-     No-slip BCs impose a drag at walls...
       CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
        ENDDO
       ENDDO
      ENDIF
C-    No-slip BCs impose a drag at bottom
      IF (momViscosity.AND.bottomDragTerms) THEN
       CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   Forcing term
      IF (momForcing)
     &  CALL EXTERNAL_FORCING_U(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)

C--   Metric terms for curvilinear grid systems
      IF (useNHMTerms) THEN
C      o Non-hydrosatic metric terms
       CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
        ENDDO
       ENDDO
      ENDIF
      IF (usingSphericalPolarMTerms) THEN
       CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   Set du/dt on boundaries to zero
      DO j=jMin,jMax
       DO i=iMin,iMax
        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
       ENDDO
      ENDDO


C---- Meridional momentum equation starts here

C     Bi-harmonic term del^2 V -> v4F
      IF (momViscosity)
     & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)

C---  Calculate mean and eddy fluxes between cells for meridional flow.

C--   Zonal flux (fZon is at west face of "v" cell)

C     Mean flow component of zonal flux -> aF
      IF (momAdvection)
     & CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,af,myThid)

C     Laplacian and bi-harmonic terms -> vF
      IF (momViscosity)
     & CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,vf,myThid)

C     Combine fluxes -> fZon
      DO j=jMin,jMax
       DO i=iMin,iMax
        fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)
       ENDDO
      ENDDO

C--   Meridional flux (fMer is at north face of "v" cell)

C     Mean flow component of meridional flux
      IF (momAdvection)
     & CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,af,myThid)

C     Laplacian and bi-harmonic term
      IF (momViscosity)
     & CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,vf,myThid)

C     Combine fluxes -> fMer
      DO j=jMin,jMax
       DO i=iMin,iMax
        fMer(i,j) = vDvdyFac*aF(i,j) + AhDvdyFac*vF(i,j)
       ENDDO
      ENDDO

C--   Vertical flux (fVer is at upper face of "v" cell)

C     o Mean flow component of vertical flux
      IF (momAdvection)
     & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,rTransV,af,myThid)

C     Eddy component of vertical flux (interior component only) -> vrF
      IF (momViscosity.AND..NOT.implicitViscosity)
     & CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)

C     Combine fluxes -> fVerV
      DO j=jMin,jMax
       DO i=iMin,iMax
        fVerV(i,j,kDown) = rVelDvdrFac*aF(i,j) + ArDvdrFac*vrF(i,j)
       ENDDO
      ENDDO

C---  Hydorstatic term (-1/rhoConst . dphi/dy )
      IF (momPressureForcing) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         pF(i,j) = -_recip_dyC(i,j,bi,bj)
     &    *(phi_hyd(i,j,k)-phi_hyd(i,j-1,k))
        ENDDO
       ENDDO
      ENDIF

C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
      DO j=jMin,jMax
       DO i=iMin,iMax
        gV(i,j,k,bi,bj) =
#ifdef OLD_UV_GEOM
     &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
     &    ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
#else
     &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rAs(i,j,bi,bj)
#endif
     &  *(fZon(i+1,j)          - fZon(i,j  )
     &   +fMer(i,j  )          - fMer(i,j-1)
     &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
     &   )
     & _PHM( +phyFac*pf(i,j) )
       ENDDO
      ENDDO

#ifdef NONLIN_FRSURF
C-- account for 3.D divergence of the flow in rStar coordinate:
      IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
     &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
     &       *vVel(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
      IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
     &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
#endif /* NONLIN_FRSURF */

C-- No-slip and drag BCs appear as body forces in cell abutting topography 
      IF (momViscosity.AND.no_slip_sides) THEN
C-     No-slip BCs impose a drag at walls...
       CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
        ENDDO
       ENDDO
      ENDIF
C-    No-slip BCs impose a drag at bottom
      IF (momViscosity.AND.bottomDragTerms) THEN
       CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   Forcing term
      IF (momForcing)
     & CALL EXTERNAL_FORCING_V(
     I     iMin,iMax,jMin,jMax,bi,bj,k,
     I     myTime,myThid)

C--   Metric terms for curvilinear grid systems
      IF (useNHMTerms) THEN
C      o Spherical polar grid metric terms
       CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
        ENDDO
       ENDDO
      ENDIF
      IF (usingSphericalPolarMTerms) THEN
       CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   Set dv/dt on boundaries to zero
      DO j=jMin,jMax
       DO i=iMin,iMax
        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
       ENDDO
      ENDDO

C--   Coriolis term
C     Note. As coded here, coriolis will not work with "thin walls"
#ifdef INCLUDE_CD_CODE
      CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,myThid)
#else
      CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
      DO j=jMin,jMax
       DO i=iMin,iMax
        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
       ENDDO
      ENDDO
      CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
      DO j=jMin,jMax
       DO i=iMin,iMax
        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
       ENDDO
      ENDDO
#endif /* INCLUDE_CD_CODE */
      IF (nonHydrostatic.OR.quasiHydrostatic) THEN
       CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
        ENDDO
       ENDDO
      ENDIF

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.7 2002/11/07 21:51:15 adcroft Exp $
2	C $Name: $
3
4	CBOI
5	C !TITLE: pkg/mom\_advdiff
6	C !AUTHORS: adcroft@mit.edu
7	C !INTRODUCTION: Flux-form Momentum Equations Package
8	C
9	C Package "mom\_fluxform" provides methods for calculating explicit terms
10	C in the momentum equation cast in flux-form:
11	C \begin{eqnarray*}
12	C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
13	C -\nabla \cdot {\bf v} u
14	C -fv
15	C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
16	C + \mbox{metrics}
17	C \\
18	C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
19	C -\nabla \cdot {\bf v} v
20	C +fu
21	C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
22	C + \mbox{metrics}
23	C \end{eqnarray*}
24	C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
25	C stresses as well as internal viscous stresses.
26	CEOI
27
28	#include "CPP_OPTIONS.h"
29
30	CBOP
31	C !ROUTINE: MOM_FLUXFORM
32
33	C !INTERFACE: ==========================================================
34	SUBROUTINE MOM_FLUXFORM(
35	I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
36	I phi_hyd,KappaRU,KappaRV,
37	U fVerU, fVerV,
38	I myTime,myIter,myThid)
39
40	C !DESCRIPTION:
41	C Calculates all the horizontal accelerations except for the implicit surface
42	C pressure gradient and implciit vertical viscosity.
43
44	C !USES: ===============================================================
45	C == Global variables ==
46	IMPLICIT NONE
47	#include "SIZE.h"
48	#include "DYNVARS.h"
49	#include "FFIELDS.h"
50	#include "EEPARAMS.h"
51	#include "PARAMS.h"
52	#include "GRID.h"
53	#include "SURFACE.h"
54
55	C !INPUT PARAMETERS: ===================================================
56	C bi,bj :: tile indices
57	C iMin,iMax,jMin,jMAx :: loop ranges
58	C k :: vertical level
59	C kUp :: =1 or 2 for consecutive k
60	C kDown :: =2 or 1 for consecutive k
61	C phi_hyd :: hydrostatic pressure (perturbation)
62	C KappaRU :: vertical viscosity
63	C KappaRV :: vertical viscosity
64	C fVerU :: vertical flux of U, 2 1/2 dim for pipe-lining
65	C fVerV :: vertical flux of V, 2 1/2 dim for pipe-lining
66	C myTime :: current time
67	C myIter :: current time-step number
68	C myThid :: thread number
69	INTEGER bi,bj,iMin,iMax,jMin,jMax
70	INTEGER k,kUp,kDown
71	_RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72	_RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73	_RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74	_RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
75	_RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
76	_RL myTime
77	INTEGER myIter
78	INTEGER myThid
79
80	C !OUTPUT PARAMETERS: ==================================================
81	C None - updates gU() and gV() in common blocks
82
83	C !LOCAL VARIABLES: ====================================================
84	C i,j :: loop indices
85	C aF :: advective flux
86	C vF :: viscous flux
87	C v4F :: bi-harmonic viscous flux
88	C vrF :: vertical viscous flux
89	C cF :: Coriolis acceleration
90	C mT :: Metric terms
91	C pF :: Pressure gradient
92	C fZon :: zonal fluxes
93	C fMer :: meridional fluxes
94	INTEGER i,j
95	_RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96	_RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97	_RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98	_RL vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99	_RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100	_RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101	_RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102	_RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103	_RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104	C wMaskOverride - Land sea flag override for top layer.
105	C afFacMom - Tracer parameters for turning terms
106	C vfFacMom on and off.
107	C pfFacMom afFacMom - Advective terms
108	C cfFacMom vfFacMom - Eddy viscosity terms
109	C mTFacMom pfFacMom - Pressure terms
110	C cfFacMom - Coriolis terms
111	C foFacMom - Forcing
112	C mTFacMom - Metric term
113	C uDudxFac, AhDudxFac, etc ... individual term tracer parameters
114	_RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115	_RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116	_RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117	_RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118	_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119	_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120	_RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121	_RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122	_RL rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123	_RL rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124	C I,J,K - Loop counters
125	C rVelMaskOverride - Factor for imposing special surface boundary conditions
126	C ( set according to free-surface condition ).
127	C hFacROpen - Lopped cell factos used tohold fraction of open
128	C hFacRClosed and closed cell wall.
129	_RL rVelMaskOverride
130	C xxxFac - On-off tracer parameters used for switching terms off.
131	_RL uDudxFac
132	_RL AhDudxFac
133	_RL A4DuxxdxFac
134	_RL vDudyFac
135	_RL AhDudyFac
136	_RL A4DuyydyFac
137	_RL rVelDudrFac
138	_RL ArDudrFac
139	_RL fuFac
140	_RL phxFac
141	_RL mtFacU
142	_RL uDvdxFac
143	_RL AhDvdxFac
144	_RL A4DvxxdxFac
145	_RL vDvdyFac
146	_RL AhDvdyFac
147	_RL A4DvyydyFac
148	_RL rVelDvdrFac
149	_RL ArDvdrFac
150	_RL fvFac
151	_RL phyFac
152	_RL vForcFac
153	_RL mtFacV
154	INTEGER km1,kp1
155	_RL wVelBottomOverride
156	LOGICAL bottomDragTerms
157	_RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158	CEOP
159
160	km1=MAX(1,k-1)
161	kp1=MIN(Nr,k+1)
162	rVelMaskOverride=1.
163	IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac
164	wVelBottomOverride=1.
165	IF (k.EQ.Nr) wVelBottomOverride=0.
166
167	C Initialise intermediate terms
168	DO J=1-OLy,sNy+OLy
169	DO I=1-OLx,sNx+OLx
170	aF(i,j) = 0.
171	vF(i,j) = 0.
172	v4F(i,j) = 0.
173	vrF(i,j) = 0.
174	cF(i,j) = 0.
175	mT(i,j) = 0.
176	pF(i,j) = 0.
177	fZon(i,j) = 0.
178	fMer(i,j) = 0.
179	rTransU(i,j) = 0.
180	rTransV(i,j) = 0.
181	ENDDO
182	ENDDO
183
184	C-- Term by term tracer parmeters
185	C o U momentum equation
186	uDudxFac = afFacMom*1.
187	AhDudxFac = vfFacMom*1.
188	A4DuxxdxFac = vfFacMom*1.
189	vDudyFac = afFacMom*1.
190	AhDudyFac = vfFacMom*1.
191	A4DuyydyFac = vfFacMom*1.
192	rVelDudrFac = afFacMom*1.
193	ArDudrFac = vfFacMom*1.
194	mTFacU = mtFacMom*1.
195	fuFac = cfFacMom*1.
196	phxFac = pfFacMom*1.
197	C o V momentum equation
198	uDvdxFac = afFacMom*1.
199	AhDvdxFac = vfFacMom*1.
200	A4DvxxdxFac = vfFacMom*1.
201	vDvdyFac = afFacMom*1.
202	AhDvdyFac = vfFacMom*1.
203	A4DvyydyFac = vfFacMom*1.
204	rVelDvdrFac = afFacMom*1.
205	ArDvdrFac = vfFacMom*1.
206	mTFacV = mtFacMom*1.
207	fvFac = cfFacMom*1.
208	phyFac = pfFacMom*1.
209	vForcFac = foFacMom*1.
210
211	IF ( no_slip_bottom
212	& .OR. bottomDragQuadratic.NE.0.
213	& .OR. bottomDragLinear.NE.0.) THEN
214	bottomDragTerms=.TRUE.
215	ELSE
216	bottomDragTerms=.FALSE.
217	ENDIF
218
219	C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP
220	IF (staggerTimeStep) THEN
221	phxFac = 0.
222	phyFac = 0.
223	ENDIF
224
225	C-- Calculate open water fraction at vorticity points
226	CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
227
228	C---- Calculate common quantities used in both U and V equations
229	C Calculate tracer cell face open areas
230	DO j=1-OLy,sNy+OLy
231	DO i=1-OLx,sNx+OLx
232	xA(i,j) = _dyG(i,j,bi,bj)
233	& drF(k)_hFacW(i,j,k,bi,bj)
234	yA(i,j) = _dxG(i,j,bi,bj)
235	& drF(k)_hFacS(i,j,k,bi,bj)
236	ENDDO
237	ENDDO
238
239	C Make local copies of horizontal flow field
240	DO j=1-OLy,sNy+OLy
241	DO i=1-OLx,sNx+OLx
242	uFld(i,j) = uVel(i,j,k,bi,bj)
243	vFld(i,j) = vVel(i,j,k,bi,bj)
244	ENDDO
245	ENDDO
246
247	C Calculate velocity field "volume transports" through tracer cell faces.
248	DO j=1-OLy,sNy+OLy
249	DO i=1-OLx,sNx+OLx
250	uTrans(i,j) = uFld(i,j)*xA(i,j)
251	vTrans(i,j) = vFld(i,j)*yA(i,j)
252	ENDDO
253	ENDDO
254
255	CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
256
257	C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
258	IF (momAdvection.AND.k.EQ.1) THEN
259
260	C- Calculate vertical transports above U & V points (West & South face):
261	CALL MOM_CALC_RTRANS( k, bi, bj,
262	O rTransU, rTransV,
263	I myTime, myIter, myThid)
264
265	C- Free surface correction term (flux at k=1)
266	CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
267	DO j=jMin,jMax
268	DO i=iMin,iMax
269	fVerU(i,j,kUp) = af(i,j)
270	ENDDO
271	ENDDO
272
273	CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
274	DO j=jMin,jMax
275	DO i=iMin,iMax
276	fVerV(i,j,kUp) = af(i,j)
277	ENDDO
278	ENDDO
279
280	C--- endif momAdvection & k=1
281	ENDIF
282
283
284	C--- Calculate vertical transports (at k+1) below U & V points :
285	IF (momAdvection) THEN
286	CALL MOM_CALC_RTRANS( k+1, bi, bj,
287	O rTransU, rTransV,
288	I myTime, myIter, myThid)
289	ENDIF
290
291
292	C---- Zonal momentum equation starts here
293
294	C Bi-harmonic term del^2 U -> v4F
295	IF (momViscosity)
296	& CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
297
298	C--- Calculate mean and eddy fluxes between cells for zonal flow.
299
300	C-- Zonal flux (fZon is at east face of "u" cell)
301
302	C Mean flow component of zonal flux -> aF
303	IF (momAdvection)
304	& CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,aF,myThid)
305
306	C Laplacian and bi-harmonic terms -> vF
307	IF (momViscosity)
308	& CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,vF,myThid)
309
310	C Combine fluxes -> fZon
311	DO j=jMin,jMax
312	DO i=iMin,iMax
313	fZon(i,j) = uDudxFacaF(i,j) + AhDudxFacvF(i,j)
314	ENDDO
315	ENDDO
316
317	C-- Meridional flux (fMer is at south face of "u" cell)
318
319	C Mean flow component of meridional flux
320	IF (momAdvection)
321	& CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,aF,myThid)
322
323	C Laplacian and bi-harmonic term
324	IF (momViscosity)
325	& CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
326
327	C Combine fluxes -> fMer
328	DO j=jMin,jMax
329	DO i=iMin,iMax
330	fMer(i,j) = vDudyFacaF(i,j) + AhDudyFacvF(i,j)
331	ENDDO
332	ENDDO
333
334	C-- Vertical flux (fVer is at upper face of "u" cell)
335
336	C Mean flow component of vertical flux (at k+1) -> aF
337	IF (momAdvection)
338	& CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,rTransU,af,myThid)
339
340	C Eddy component of vertical flux (interior component only) -> vrF
341	IF (momViscosity.AND..NOT.implicitViscosity)
342	& CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)
343
344	C Combine fluxes
345	DO j=jMin,jMax
346	DO i=iMin,iMax
347	fVerU(i,j,kDown) = rVelDudrFacaF(i,j) + ArDudrFacvrF(i,j)
348	ENDDO
349	ENDDO
350
351	C--- Hydrostatic term ( -1/rhoConst . dphi/dx )
352	IF (momPressureForcing) THEN
353	DO j=jMin,jMax
354	DO i=iMin,iMax
355	pf(i,j) = - _recip_dxC(i,j,bi,bj)
356	& *(phi_hyd(i,j,k)-phi_hyd(i-1,j,k))
357	ENDDO
358	ENDDO
359	ENDIF
360
361	C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
362	DO j=jMin,jMax
363	DO i=iMin,iMax
364	gU(i,j,k,bi,bj) =
365	#ifdef OLD_UV_GEOM
366	& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
367	& ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
368	#else
369	& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
370	& *recip_rAw(i,j,bi,bj)
371	#endif
372	& *(fZon(i,j ) - fZon(i-1,j)
373	& +fMer(i,j+1) - fMer(i ,j)
374	& +fVerU(i,j,kUp)rkFac - fVerU(i,j,kDown)rkFac
375	& )
376	& _PHM( +phxFac * pf(i,j) )
377	ENDDO
378	ENDDO
379
380	#ifdef NONLIN_FRSURF
381	C-- account for 3.D divergence of the flow in rStar coordinate:
382	IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
383	DO j=jMin,jMax
384	DO i=iMin,iMax
385	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
386	& - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
387	& *uVel(i,j,k,bi,bj)
388	ENDDO
389	ENDDO
390	ENDIF
391	IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
392	DO j=jMin,jMax
393	DO i=iMin,iMax
394	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
395	& - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
396	ENDDO
397	ENDDO
398	ENDIF
399	#endif /* NONLIN_FRSURF */
400
401	C-- No-slip and drag BCs appear as body forces in cell abutting topography
402	IF (momViscosity.AND.no_slip_sides) THEN
403	C- No-slip BCs impose a drag at walls...
404	CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
405	DO j=jMin,jMax
406	DO i=iMin,iMax
407	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
408	ENDDO
409	ENDDO
410	ENDIF
411	C- No-slip BCs impose a drag at bottom
412	IF (momViscosity.AND.bottomDragTerms) THEN
413	CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
414	DO j=jMin,jMax
415	DO i=iMin,iMax
416	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
417	ENDDO
418	ENDDO
419	ENDIF
420
421	C-- Forcing term
422	IF (momForcing)
423	& CALL EXTERNAL_FORCING_U(
424	I iMin,iMax,jMin,jMax,bi,bj,k,
425	I myTime,myThid)
426
427	C-- Metric terms for curvilinear grid systems
428	IF (useNHMTerms) THEN
429	C o Non-hydrosatic metric terms
430	CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
431	DO j=jMin,jMax
432	DO i=iMin,iMax
433	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
434	ENDDO
435	ENDDO
436	ENDIF
437	IF (usingSphericalPolarMTerms) THEN
438	CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
439	DO j=jMin,jMax
440	DO i=iMin,iMax
441	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
442	ENDDO
443	ENDDO
444	ENDIF
445
446	C-- Set du/dt on boundaries to zero
447	DO j=jMin,jMax
448	DO i=iMin,iMax
449	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
450	ENDDO
451	ENDDO
452
453
454	C---- Meridional momentum equation starts here
455
456	C Bi-harmonic term del^2 V -> v4F
457	IF (momViscosity)
458	& CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
459
460	C--- Calculate mean and eddy fluxes between cells for meridional flow.
461
462	C-- Zonal flux (fZon is at west face of "v" cell)
463
464	C Mean flow component of zonal flux -> aF
465	IF (momAdvection)
466	& CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,af,myThid)
467
468	C Laplacian and bi-harmonic terms -> vF
469	IF (momViscosity)
470	& CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,vf,myThid)
471
472	C Combine fluxes -> fZon
473	DO j=jMin,jMax
474	DO i=iMin,iMax
475	fZon(i,j) = uDvdxFacaF(i,j) + AhDvdxFacvF(i,j)
476	ENDDO
477	ENDDO
478
479	C-- Meridional flux (fMer is at north face of "v" cell)
480
481	C Mean flow component of meridional flux
482	IF (momAdvection)
483	& CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,af,myThid)
484
485	C Laplacian and bi-harmonic term
486	IF (momViscosity)
487	& CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,vf,myThid)
488
489	C Combine fluxes -> fMer
490	DO j=jMin,jMax
491	DO i=iMin,iMax
492	fMer(i,j) = vDvdyFacaF(i,j) + AhDvdyFacvF(i,j)
493	ENDDO
494	ENDDO
495
496	C-- Vertical flux (fVer is at upper face of "v" cell)
497
498	C o Mean flow component of vertical flux
499	IF (momAdvection)
500	& CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,rTransV,af,myThid)
501
502	C Eddy component of vertical flux (interior component only) -> vrF
503	IF (momViscosity.AND..NOT.implicitViscosity)
504	& CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)
505
506	C Combine fluxes -> fVerV
507	DO j=jMin,jMax
508	DO i=iMin,iMax
509	fVerV(i,j,kDown) = rVelDvdrFacaF(i,j) + ArDvdrFacvrF(i,j)
510	ENDDO
511	ENDDO
512
513	C--- Hydorstatic term (-1/rhoConst . dphi/dy )
514	IF (momPressureForcing) THEN
515	DO j=jMin,jMax
516	DO i=iMin,iMax
517	pF(i,j) = -_recip_dyC(i,j,bi,bj)
518	& *(phi_hyd(i,j,k)-phi_hyd(i,j-1,k))
519	ENDDO
520	ENDDO
521	ENDIF
522
523	C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
524	DO j=jMin,jMax
525	DO i=iMin,iMax
526	gV(i,j,k,bi,bj) =
527	#ifdef OLD_UV_GEOM
528	& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
529	& ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
530	#else
531	& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
532	& *recip_rAs(i,j,bi,bj)
533	#endif
534	& *(fZon(i+1,j) - fZon(i,j )
535	& +fMer(i,j ) - fMer(i,j-1)
536	& +fVerV(i,j,kUp)rkFac - fVerV(i,j,kDown)rkFac
537	& )
538	& _PHM( +phyFac*pf(i,j) )
539	ENDDO
540	ENDDO
541
542	#ifdef NONLIN_FRSURF
543	C-- account for 3.D divergence of the flow in rStar coordinate:
544	IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
545	DO j=jMin,jMax
546	DO i=iMin,iMax
547	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
548	& - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
549	& *vVel(i,j,k,bi,bj)
550	ENDDO
551	ENDDO
552	ENDIF
553	IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
554	DO j=jMin,jMax
555	DO i=iMin,iMax
556	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
557	& - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
558	ENDDO
559	ENDDO
560	ENDIF
561	#endif /* NONLIN_FRSURF */
562
563	C-- No-slip and drag BCs appear as body forces in cell abutting topography
564	IF (momViscosity.AND.no_slip_sides) THEN
565	C- No-slip BCs impose a drag at walls...
566	CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid)
567	DO j=jMin,jMax
568	DO i=iMin,iMax
569	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
570	ENDDO
571	ENDDO
572	ENDIF
573	C- No-slip BCs impose a drag at bottom
574	IF (momViscosity.AND.bottomDragTerms) THEN
575	CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
576	DO j=jMin,jMax
577	DO i=iMin,iMax
578	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
579	ENDDO
580	ENDDO
581	ENDIF
582
583	C-- Forcing term
584	IF (momForcing)
585	& CALL EXTERNAL_FORCING_V(
586	I iMin,iMax,jMin,jMax,bi,bj,k,
587	I myTime,myThid)
588
589	C-- Metric terms for curvilinear grid systems
590	IF (useNHMTerms) THEN
591	C o Spherical polar grid metric terms
592	CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
593	DO j=jMin,jMax
594	DO i=iMin,iMax
595	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
596	ENDDO
597	ENDDO
598	ENDIF
599	IF (usingSphericalPolarMTerms) THEN
600	CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
601	DO j=jMin,jMax
602	DO i=iMin,iMax
603	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
604	ENDDO
605	ENDDO
606	ENDIF
607
608	C-- Set dv/dt on boundaries to zero
609	DO j=jMin,jMax
610	DO i=iMin,iMax
611	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
612	ENDDO
613	ENDDO
614
615	C-- Coriolis term
616	C Note. As coded here, coriolis will not work with "thin walls"
617	#ifdef INCLUDE_CD_CODE
618	CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,myThid)
619	#else
620	CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
621	DO j=jMin,jMax
622	DO i=iMin,iMax
623	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
624	ENDDO
625	ENDDO
626	CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
627	DO j=jMin,jMax
628	DO i=iMin,iMax
629	gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
630	ENDDO
631	ENDDO
632	#endif /* INCLUDE_CD_CODE */
633	IF (nonHydrostatic.OR.quasiHydrostatic) THEN
634	CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
635	DO j=jMin,jMax
636	DO i=iMin,iMax
637	gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
638	ENDDO
639	ENDDO
640	ENDIF
641
642	RETURN
643	END