C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.22 2009/11/28 20:59:18 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: INTEGR_CONTINUITY
C !INTERFACE:
SUBROUTINE INTEGR_CONTINUITY(
I bi, bj, uFld, vFld,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INTEGR_CONTINUITY
C | o Integrate the continuity Eq : compute vertical velocity
C | and free surface "r-anomaly" (etaN) to satisfy
C | exactly the convervation of the Total Volume
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C uFld :: Zonal velocity ( m/s )
C vFld :: Meridional velocity ( m/s )
C bi,bj :: tile index
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: Thread number for this instance of the routine.
_RL myTime
INTEGER myIter
INTEGER myThid
INTEGER bi,bj
_RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
C !LOCAL VARIABLES:
C Local variables in common block
C Local variables
C i,j,k :: Loop counters
C uTrans :: Volume transports ( uVel.xA )
C vTrans :: Volume transports ( vVel.yA )
C hDivFlow :: Div. Barotropic Flow [transport unit m3/s]
INTEGER k
#ifdef EXACT_CONSERV
INTEGER i,j, ks
_RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL hDivFlow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
_RL facEmP, facMass
#endif /* EXACT_CONSERV */
#ifndef ALLOW_ADDFLUID
_RL addMass(1)
#endif /* ndef ALLOW_ADDFLUID */
CEOP
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef EXACT_CONSERV
IF (exactConserv) THEN
facEmP = 0.
IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit
facMass = 0.
IF ( selectAddFluid.GE.1 ) facMass = mass2rUnit
C-- Compute the Divergence of The Barotropic Flow :
C- Initialise
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
hDivFlow(i,j) = 0. _d 0
utrans(i,j) = 0. _d 0
vtrans(i,j) = 0. _d 0
ENDDO
ENDDO
DO k=1,Nr
C- Calculate velocity field "volume transports" through tracer cell faces
C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport).
DO j=1,sNy+1
DO i=1,sNx+1
uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
& *deepFacC(k)*rhoFacC(k)
& *drF(k)*_hFacW(i,j,k,bi,bj)
vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
& *deepFacC(k)*rhoFacC(k)
& *drF(k)*_hFacS(i,j,k,bi,bj)
ENDDO
ENDDO
C- Integrate vertically the Horizontal Divergence
DO j=1,sNy
DO i=1,sNx
hDivFlow(i,j) = hDivFlow(i,j)
& +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
& +vTrans(i,j+1)-vTrans(i,j) )
#ifdef ALLOW_ADDFLUID
& -facMass*addMass(i,j,k,bi,bj)
#endif /* ALLOW_ADDFLUID */
ENDDO
ENDDO
C- End DO k=1,Nr
ENDDO
C------------------------------------
C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ;
C anelastic: always assumes that rhoFacF(1) = 1
IF ( myTime.EQ.startTime .AND. myTime.NE.baseTime
& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
C needs previous time-step value of E-P-R, that has not been loaded
C and was not in old pickup-file ; try to use etaN & etaH instead.
IF ( usePickupBeforeC54 ) THEN
DO j=1,sNy
DO i=1,sNx
dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
& / (implicDiv2Dflow*deltaTfreesurf)
ENDDO
ENDDO
ENDIF
DO j=1,sNy
DO i=1,sNx
PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
& + hDivFlow(i,j)*recip_rA(i,j,bi,bj)
& *recip_deepFac2F(1)
PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass
ENDDO
ENDDO
ELSEIF ( myTime.EQ.startTime ) THEN
DO j=1,sNy
DO i=1,sNx
ks = ksurfC(I,J,bi,bj)
PmEpR(i,j,bi,bj) = 0. _d 0
dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
& *recip_deepFac2F(ks)
ENDDO
ENDDO
ELSE
C-- Needs to get valid PmEpR (for T,S forcing) also in overlap regions
C (e.g., if using KPP) => set over full index range
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
ENDDO
ENDDO
DO j=1,sNy
DO i=1,sNx
ks = ksurfC(i,j,bi,bj)
dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
& *recip_deepFac2F(ks)
& -facEmP*EmPmR(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
C------------------------------------
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF (exactConserv .AND. myTime.NE.startTime) THEN
C-- Update etaN at the end of the time step :
C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
ENDDO
ENDDO
ELSE
DO j=1,sNy
DO i=1,sNx
etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
& + implicDiv2Dflow*dEtaHdt(i,j,bi,bj)*deltaTfreesurf
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_OBCS
C-- Apply OBC to etaN if NonLin-FreeSurf, reset to zero otherwise:
IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
#endif /* ALLOW_OBCS */
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# ifdef NONLIN_FRSURF
IF (select_rStar .NE. 0) THEN
# ifndef DISABLE_RSTAR_CODE
C-- note: final value of rStarDhCDt from S/R CALC_R_STAR
C will be different (no deep-factor, no rho factor)
DO j=1,sNy
DO i=1,sNx
ks = ksurfC(i,j,bi,bj)
rStarDhCDt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
& *deepFac2F(ks)*rhoFacF(ks)
& *recip_Rcol(i,j,bi,bj)
ENDDO
ENDDO
# endif /* DISABLE_RSTAR_CODE */
ENDIF
# endif /* NONLIN_FRSURF */
#endif /* EXACT_CONSERV */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
DO k=Nr,1,-1
C-- Integrate continuity vertically for vertical velocity
CALL INTEGRATE_FOR_W(
I bi, bj, k, uFld, vFld, addMass,
O wVel,
I myThid )
#ifdef EXACT_CONSERV
IF ( k.EQ.Nr .AND. myTime.NE.baseTime .AND. usingPCoords
& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
DO j=1,sNy
DO i=1,sNx
wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
& +mass2rUnit*PmEpR(i,j,bi,bj)*maskC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* EXACT_CONSERV */
#ifdef ALLOW_OBCS
C-- Apply OBC to W if in N-H mode, reset to zero otherwise:
IF ( useOBCS ) CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
#endif /* ALLOW_OBCS */
C- End DO k=Nr,1,-1
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
RETURN
END