pkg/mom_fluxform/mom_fluxform.F

C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.11 2003/02/18 15:36:45 jmc 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        dPhihydX,dPhiHydY,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  dPhiHydX,Y           :: Gradient (X & Y dir.) of Hydrostatic Potential 
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 dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _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 .AND. viscA4.NE.0. ) 
     & 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+1
       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--   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
     &   )
     &  - phxFac*dPhiHydX(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 (moved to timestep.F)
c     IF (momForcing)
c    &  CALL EXTERNAL_FORCING_U(
c    I     iMin,iMax,jMin,jMax,bi,bj,k,
c    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 .AND. viscA4.NE.0. ) 
     & 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+1
        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--   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
     &   )
     &  - phyFac*dPhiHydY(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 (moved to timestep.F)
c     IF (momForcing)
c    & CALL EXTERNAL_FORCING_V(
c    I     iMin,iMax,jMin,jMax,bi,bj,k,
c    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"
c     IF (useCDscheme) THEN
c       CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
c     ELSE
      IF (.NOT.useCDscheme) THEN
        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

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