pkg/mom_fluxform/mom_fluxform.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/mom_fluxform/mom_fluxform.F,v 1.2 2001/08/17 18:40:30 adcroft Exp $
C $Name:  $

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

#include "CPP_OPTIONS.h"

CBOP
C !ROUTINE: MOM_FLUXFORM

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

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

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

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  iMin,iMax,jMin,jMAx  :: loop ranges
C  k                    :: vertical level
C  kUp                  :: =1 or 2 for consecutive k
C  kDown                :: =2 or 1 for consecutive k
C  phi_hyd              :: hydrostatic pressure (perturbation)
C  KappaRU              :: vertical viscosity
C  KappaRV              :: vertical viscosity
C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining
C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining
C  myCurrentTime        :: current time
C  myIter               :: current time-step number
C  myThid               :: thread number
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k,kUp,kDown
      _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RL     myCurrentTime
      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)
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.
       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---- 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--   Free surface correction term (flux at k=1)
      IF (momAdvection.AND.k.EQ.1) THEN
       CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,af,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         fVerU(i,j,kUp) = af(i,j)
        ENDDO
       ENDDO
      ENDIF
C     Mean flow component of vertical flux (at k+1) -> aF
      IF (momAdvection)
     & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,af,myThid)

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

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

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

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

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     myCurrentTime,myThid)

C--   Metric terms for curvilinear grid systems
      IF (usingSphericalPolarMTerms) THEN
C      o Spherical polar grid 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
       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--   Free surface correction term (flux at k=1)
      IF (momAdvection.AND.k.EQ.1) THEN
       CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,af,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         fVerV(i,j,kUp) = af(i,j)
        ENDDO
       ENDDO
      ENDIF
C     o Mean flow component of vertical flux
      IF (momAdvection)
     & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,af,myThid)

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

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

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

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

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     myCurrentTime,myThid)

C--   Metric terms for curvilinear grid systems
      IF (usingSphericalPolarMTerms) 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
       CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
       DO j=jMin,jMax
        DO i=iMin,iMax
         gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
        ENDDO
       ENDDO
      ENDIF

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

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

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