pkg/mom_fluxform/mom_fluxform.F

C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.14 2003/10/10 23:00:01 heimbach 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 "MOM_FLUXFORM_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.
#ifdef ALLOW_AUTODIFF_TAMC
C- jmc: this is wrong, but at least with #ifdef/endif TAMC, it does not break 
C       the forward code ; (same thing in mom_vectinv)
        fVerU(i,j,1) = 0. _d 0
        fVerU(i,j,2) = 0. _d 0
        fVerV(i,j,1) = 0. _d 0
        fVerV(i,j,2) = 0. _d 0
#endif
       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	jmc	1.15	C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.14 2003/10/10 23:00:01 heimbach Exp $
2	adcroft	1.2	C $Name: $
3	adcroft	1.1
4	adcroft	1.3	CBOI
5			C !TITLE: pkg/mom\_advdiff
6			C !AUTHORS: adcroft@mit.edu
7	adcroft	1.4	C !INTRODUCTION: Flux-form Momentum Equations Package
8	adcroft	1.3	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	edhill	1.13	#include "MOM_FLUXFORM_OPTIONS.h"
29	adcroft	1.1
30	adcroft	1.3	CBOP
31			C !ROUTINE: MOM_FLUXFORM
32
33			C !INTERFACE: ==========================================================
34	adcroft	1.1	SUBROUTINE MOM_FLUXFORM(
35			I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
36	jmc	1.11	I dPhihydX,dPhiHydY,KappaRU,KappaRV,
37	adcroft	1.1	U fVerU, fVerV,
38	jmc	1.8	I myTime,myIter,myThid)
39	adcroft	1.3
40			C !DESCRIPTION:
41			C Calculates all the horizontal accelerations except for the implicit surface
42			C pressure gradient and implciit vertical viscosity.
43	adcroft	1.1
44	adcroft	1.3	C !USES: ===============================================================
45	adcroft	1.1	C == Global variables ==
46	adcroft	1.3	IMPLICIT NONE
47	adcroft	1.1	#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	adcroft	1.3	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	jmc	1.9	C dPhiHydX,Y :: Gradient (X & Y dir.) of Hydrostatic Potential
62	adcroft	1.3	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	jmc	1.8	C myTime :: current time
67	adcroft	1.3	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	jmc	1.9	_RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72			_RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73	adcroft	1.1	_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	jmc	1.8	_RL myTime
78	adcroft	1.2	INTEGER myIter
79	adcroft	1.1	INTEGER myThid
80
81	adcroft	1.3	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	adcroft	1.1	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	jmc	1.8	_RL rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124			_RL rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125	adcroft	1.1	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	adcroft	1.3	CEOP
160	adcroft	1.1
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	jmc	1.8	rTransU(i,j) = 0.
181			rTransV(i,j) = 0.
182	jmc	1.15	#ifdef ALLOW_AUTODIFF_TAMC
183			C- jmc: this is wrong, but at least with #ifdef/endif TAMC, it does not break
184			C the forward code ; (same thing in mom_vectinv)
185	heimbach	1.14	fVerU(i,j,1) = 0. _d 0
186			fVerU(i,j,2) = 0. _d 0
187			fVerV(i,j,1) = 0. _d 0
188			fVerV(i,j,2) = 0. _d 0
189	jmc	1.15	#endif
190	adcroft	1.1	ENDDO
191			ENDDO
192
193			C-- Term by term tracer parmeters
194			C o U momentum equation
195			uDudxFac = afFacMom*1.
196			AhDudxFac = vfFacMom*1.
197			A4DuxxdxFac = vfFacMom*1.
198			vDudyFac = afFacMom*1.
199			AhDudyFac = vfFacMom*1.
200			A4DuyydyFac = vfFacMom*1.
201			rVelDudrFac = afFacMom*1.
202			ArDudrFac = vfFacMom*1.
203			mTFacU = mtFacMom*1.
204			fuFac = cfFacMom*1.
205			phxFac = pfFacMom*1.
206			C o V momentum equation
207			uDvdxFac = afFacMom*1.
208			AhDvdxFac = vfFacMom*1.
209			A4DvxxdxFac = vfFacMom*1.
210			vDvdyFac = afFacMom*1.
211			AhDvdyFac = vfFacMom*1.
212			A4DvyydyFac = vfFacMom*1.
213			rVelDvdrFac = afFacMom*1.
214			ArDvdrFac = vfFacMom*1.
215			mTFacV = mtFacMom*1.
216			fvFac = cfFacMom*1.
217			phyFac = pfFacMom*1.
218			vForcFac = foFacMom*1.
219
220			IF ( no_slip_bottom
221			& .OR. bottomDragQuadratic.NE.0.
222			& .OR. bottomDragLinear.NE.0.) THEN
223			bottomDragTerms=.TRUE.
224			ELSE
225			bottomDragTerms=.FALSE.
226			ENDIF
227
228			C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP
229			IF (staggerTimeStep) THEN
230			phxFac = 0.
231			phyFac = 0.
232			ENDIF
233
234			C-- Calculate open water fraction at vorticity points
235			CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
236
237			C---- Calculate common quantities used in both U and V equations
238			C Calculate tracer cell face open areas
239			DO j=1-OLy,sNy+OLy
240			DO i=1-OLx,sNx+OLx
241			xA(i,j) = _dyG(i,j,bi,bj)
242			& drF(k)_hFacW(i,j,k,bi,bj)
243			yA(i,j) = _dxG(i,j,bi,bj)
244			& drF(k)_hFacS(i,j,k,bi,bj)
245			ENDDO
246			ENDDO
247
248			C Make local copies of horizontal flow field
249			DO j=1-OLy,sNy+OLy
250			DO i=1-OLx,sNx+OLx
251			uFld(i,j) = uVel(i,j,k,bi,bj)
252			vFld(i,j) = vVel(i,j,k,bi,bj)
253			ENDDO
254			ENDDO
255
256			C Calculate velocity field "volume transports" through tracer cell faces.
257			DO j=1-OLy,sNy+OLy
258			DO i=1-OLx,sNx+OLx
259			uTrans(i,j) = uFld(i,j)*xA(i,j)
260			vTrans(i,j) = vFld(i,j)*yA(i,j)
261			ENDDO
262			ENDDO
263
264			CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
265
266	jmc	1.8	C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
267			IF (momAdvection.AND.k.EQ.1) THEN
268
269			C- Calculate vertical transports above U & V points (West & South face):
270			CALL MOM_CALC_RTRANS( k, bi, bj,
271			O rTransU, rTransV,
272			I myTime, myIter, myThid)
273
274			C- Free surface correction term (flux at k=1)
275			CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
276			DO j=jMin,jMax
277			DO i=iMin,iMax
278			fVerU(i,j,kUp) = af(i,j)
279			ENDDO
280			ENDDO
281
282			CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
283			DO j=jMin,jMax
284			DO i=iMin,iMax
285			fVerV(i,j,kUp) = af(i,j)
286			ENDDO
287			ENDDO
288
289			C--- endif momAdvection & k=1
290			ENDIF
291
292
293			C--- Calculate vertical transports (at k+1) below U & V points :
294			IF (momAdvection) THEN
295			CALL MOM_CALC_RTRANS( k+1, bi, bj,
296			O rTransU, rTransV,
297			I myTime, myIter, myThid)
298			ENDIF
299
300
301	adcroft	1.1	C---- Zonal momentum equation starts here
302
303			C Bi-harmonic term del^2 U -> v4F
304	jmc	1.10	IF (momViscosity .AND. viscA4.NE.0. )
305	adcroft	1.1	& CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
306
307			C--- Calculate mean and eddy fluxes between cells for zonal flow.
308
309			C-- Zonal flux (fZon is at east face of "u" cell)
310
311			C Mean flow component of zonal flux -> aF
312			IF (momAdvection)
313			& CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,aF,myThid)
314
315			C Laplacian and bi-harmonic terms -> vF
316			IF (momViscosity)
317			& CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,vF,myThid)
318
319			C Combine fluxes -> fZon
320			DO j=jMin,jMax
321			DO i=iMin,iMax
322			fZon(i,j) = uDudxFacaF(i,j) + AhDudxFacvF(i,j)
323			ENDDO
324			ENDDO
325
326			C-- Meridional flux (fMer is at south face of "u" cell)
327
328			C Mean flow component of meridional flux
329			IF (momAdvection)
330			& CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,aF,myThid)
331
332			C Laplacian and bi-harmonic term
333			IF (momViscosity)
334			& CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
335
336			C Combine fluxes -> fMer
337	jmc	1.10	DO j=jMin,jMax+1
338	adcroft	1.1	DO i=iMin,iMax
339			fMer(i,j) = vDudyFacaF(i,j) + AhDudyFacvF(i,j)
340			ENDDO
341			ENDDO
342
343			C-- Vertical flux (fVer is at upper face of "u" cell)
344
345			C Mean flow component of vertical flux (at k+1) -> aF
346			IF (momAdvection)
347	jmc	1.8	& CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,rTransU,af,myThid)
348	adcroft	1.1
349			C Eddy component of vertical flux (interior component only) -> vrF
350			IF (momViscosity.AND..NOT.implicitViscosity)
351			& CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)
352
353			C Combine fluxes
354			DO j=jMin,jMax
355			DO i=iMin,iMax
356			fVerU(i,j,kDown) = rVelDudrFacaF(i,j) + ArDudrFacvrF(i,j)
357			ENDDO
358			ENDDO
359
360			C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
361			DO j=jMin,jMax
362			DO i=iMin,iMax
363			gU(i,j,k,bi,bj) =
364			#ifdef OLD_UV_GEOM
365			& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
366			& ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
367			#else
368			& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
369			& *recip_rAw(i,j,bi,bj)
370			#endif
371			& *(fZon(i,j ) - fZon(i-1,j)
372			& +fMer(i,j+1) - fMer(i ,j)
373			& +fVerU(i,j,kUp)rkFac - fVerU(i,j,kDown)rkFac
374			& )
375	jmc	1.9	& - phxFac*dPhiHydX(i,j)
376	adcroft	1.1	ENDDO
377			ENDDO
378
379	jmc	1.8	#ifdef NONLIN_FRSURF
380			C-- account for 3.D divergence of the flow in rStar coordinate:
381			IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
382			DO j=jMin,jMax
383			DO i=iMin,iMax
384			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
385			& - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
386			& *uVel(i,j,k,bi,bj)
387			ENDDO
388			ENDDO
389			ENDIF
390			IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
391			DO j=jMin,jMax
392			DO i=iMin,iMax
393			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
394			& - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
395			ENDDO
396			ENDDO
397			ENDIF
398			#endif /* NONLIN_FRSURF */
399
400	adcroft	1.1	C-- No-slip and drag BCs appear as body forces in cell abutting topography
401			IF (momViscosity.AND.no_slip_sides) THEN
402			C- No-slip BCs impose a drag at walls...
403			CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
404			DO j=jMin,jMax
405			DO i=iMin,iMax
406			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
407			ENDDO
408			ENDDO
409			ENDIF
410			C- No-slip BCs impose a drag at bottom
411			IF (momViscosity.AND.bottomDragTerms) THEN
412			CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
413			DO j=jMin,jMax
414			DO i=iMin,iMax
415			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
416			ENDDO
417			ENDDO
418			ENDIF
419
420	jmc	1.12	C-- Forcing term (moved to timestep.F)
421			c IF (momForcing)
422			c & CALL EXTERNAL_FORCING_U(
423			c I iMin,iMax,jMin,jMax,bi,bj,k,
424			c I myTime,myThid)
425	adcroft	1.1
426			C-- Metric terms for curvilinear grid systems
427	adcroft	1.5	IF (useNHMTerms) THEN
428			C o Non-hydrosatic metric terms
429	adcroft	1.1	CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,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	adcroft	1.5	ENDIF
436			IF (usingSphericalPolarMTerms) THEN
437	adcroft	1.1	CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
438			DO j=jMin,jMax
439			DO i=iMin,iMax
440			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
441			ENDDO
442			ENDDO
443			ENDIF
444
445			C-- Set du/dt on boundaries to zero
446			DO j=jMin,jMax
447			DO i=iMin,iMax
448			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
449			ENDDO
450			ENDDO
451
452
453			C---- Meridional momentum equation starts here
454
455			C Bi-harmonic term del^2 V -> v4F
456	jmc	1.10	IF (momViscosity .AND. viscA4.NE.0. )
457	adcroft	1.1	& CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
458
459			C--- Calculate mean and eddy fluxes between cells for meridional flow.
460
461			C-- Zonal flux (fZon is at west face of "v" cell)
462
463			C Mean flow component of zonal flux -> aF
464			IF (momAdvection)
465			& CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,af,myThid)
466
467			C Laplacian and bi-harmonic terms -> vF
468			IF (momViscosity)
469			& CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,vf,myThid)
470
471			C Combine fluxes -> fZon
472			DO j=jMin,jMax
473	jmc	1.10	DO i=iMin,iMax+1
474	adcroft	1.1	fZon(i,j) = uDvdxFacaF(i,j) + AhDvdxFacvF(i,j)
475			ENDDO
476			ENDDO
477
478			C-- Meridional flux (fMer is at north face of "v" cell)
479
480			C Mean flow component of meridional flux
481			IF (momAdvection)
482			& CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,af,myThid)
483
484			C Laplacian and bi-harmonic term
485			IF (momViscosity)
486			& CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,vf,myThid)
487
488			C Combine fluxes -> fMer
489			DO j=jMin,jMax
490			DO i=iMin,iMax
491			fMer(i,j) = vDvdyFacaF(i,j) + AhDvdyFacvF(i,j)
492			ENDDO
493			ENDDO
494
495			C-- Vertical flux (fVer is at upper face of "v" cell)
496
497			C o Mean flow component of vertical flux
498			IF (momAdvection)
499	jmc	1.8	& CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,rTransV,af,myThid)
500	adcroft	1.1
501			C Eddy component of vertical flux (interior component only) -> vrF
502			IF (momViscosity.AND..NOT.implicitViscosity)
503			& CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)
504
505			C Combine fluxes -> fVerV
506			DO j=jMin,jMax
507			DO i=iMin,iMax
508			fVerV(i,j,kDown) = rVelDvdrFacaF(i,j) + ArDvdrFacvrF(i,j)
509			ENDDO
510			ENDDO
511
512			C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
513			DO j=jMin,jMax
514			DO i=iMin,iMax
515			gV(i,j,k,bi,bj) =
516			#ifdef OLD_UV_GEOM
517			& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
518			& ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
519			#else
520			& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
521			& *recip_rAs(i,j,bi,bj)
522			#endif
523			& *(fZon(i+1,j) - fZon(i,j )
524			& +fMer(i,j ) - fMer(i,j-1)
525			& +fVerV(i,j,kUp)rkFac - fVerV(i,j,kDown)rkFac
526			& )
527	jmc	1.9	& - phyFac*dPhiHydY(i,j)
528	adcroft	1.1	ENDDO
529			ENDDO
530
531	jmc	1.8	#ifdef NONLIN_FRSURF
532			C-- account for 3.D divergence of the flow in rStar coordinate:
533			IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
534			DO j=jMin,jMax
535			DO i=iMin,iMax
536			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
537			& - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
538			& *vVel(i,j,k,bi,bj)
539			ENDDO
540			ENDDO
541			ENDIF
542			IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
543			DO j=jMin,jMax
544			DO i=iMin,iMax
545			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
546			& - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
547			ENDDO
548			ENDDO
549			ENDIF
550			#endif /* NONLIN_FRSURF */
551
552	adcroft	1.1	C-- No-slip and drag BCs appear as body forces in cell abutting topography
553			IF (momViscosity.AND.no_slip_sides) THEN
554			C- No-slip BCs impose a drag at walls...
555			CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid)
556			DO j=jMin,jMax
557			DO i=iMin,iMax
558			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
559			ENDDO
560			ENDDO
561			ENDIF
562			C- No-slip BCs impose a drag at bottom
563			IF (momViscosity.AND.bottomDragTerms) THEN
564			CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
565			DO j=jMin,jMax
566			DO i=iMin,iMax
567			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
568			ENDDO
569			ENDDO
570			ENDIF
571
572	jmc	1.12	C-- Forcing term (moved to timestep.F)
573			c IF (momForcing)
574			c & CALL EXTERNAL_FORCING_V(
575			c I iMin,iMax,jMin,jMax,bi,bj,k,
576			c I myTime,myThid)
577	adcroft	1.1
578			C-- Metric terms for curvilinear grid systems
579	adcroft	1.5	IF (useNHMTerms) THEN
580	adcroft	1.1	C o Spherical polar grid metric terms
581			CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,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	adcroft	1.5	ENDIF
588			IF (usingSphericalPolarMTerms) THEN
589	adcroft	1.1	CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
590			DO j=jMin,jMax
591			DO i=iMin,iMax
592			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
593			ENDDO
594			ENDDO
595			ENDIF
596
597			C-- Set dv/dt on boundaries to zero
598			DO j=jMin,jMax
599			DO i=iMin,iMax
600			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
601			ENDDO
602			ENDDO
603
604			C-- Coriolis term
605			C Note. As coded here, coriolis will not work with "thin walls"
606	jmc	1.12	c IF (useCDscheme) THEN
607			c CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
608			c ELSE
609			IF (.NOT.useCDscheme) THEN
610			CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
611			DO j=jMin,jMax
612			DO i=iMin,iMax
613			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
614			ENDDO
615			ENDDO
616			CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
617			DO j=jMin,jMax
618			DO i=iMin,iMax
619			gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
620			ENDDO
621			ENDDO
622			ENDIF
623
624	adcroft	1.7	IF (nonHydrostatic.OR.quasiHydrostatic) THEN
625	adcroft	1.6	CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
626			DO j=jMin,jMax
627			DO i=iMin,iMax
628			gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
629			ENDDO
630			ENDDO
631			ENDIF
632	adcroft	1.1
633			RETURN
634			END