pkg/mom_fluxform/mom_fluxform.F

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