C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.28 2011/12/22 19:00:58 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INTEGR_CONTINUITY
C     !INTERFACE:
      SUBROUTINE INTEGR_CONTINUITY(
     I                             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,etaH) to satisfy
C     |   exactly the conservation 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     myTime   :: Current time in simulation
C     myIter   :: Current iteration number in simulation
C     myThid   :: my Thread Id. number
      _RL myTime
      INTEGER myIter
      INTEGER myThid
      _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     bi,bj    :: tile index
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,bi,bj
#ifdef EXACT_CONSERV
      INTEGER i, j
      INTEGER 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
#else /* EXACT_CONSERV */
# ifdef ALLOW_OBCS
      INTEGER i, j
# endif
#endif /* EXACT_CONSERV */
#ifndef ALLOW_ADDFLUID
      _RL addMass(1)
#endif /* ndef ALLOW_ADDFLUID */
#if (defined NONLIN_FRSURF) && !(defined DISABLE_RSTAR_CODE)
      _RL rStarDhDt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#else
      _RL rStarDhDt(1)
#endif
CEOP

C--   Start bi,bj loops
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

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 ( myIter.EQ.nIter0 .AND. myIter.NE.0
     &    .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 ( myIter.EQ.nIter0 ) 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------------------------------------

#ifdef ALLOW_OBCS
C--    reset dEtaHdt to zero outside the OB interior region
       IF ( useOBCS ) THEN
         DO j=1,sNy
          DO i=1,sNx
            dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj)
          ENDDO
         ENDDO
       ENDIF
#endif /* ALLOW_OBCS */

C-- end if exactConserv block
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( exactConserv .AND. myIter.NE.nIter0 ) 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--    Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer
C      needed with proper masking in solver (matrix+cg2d_b,x) and masking
C      of dEtaHdt above. etaN next to OB does not enter present algorithm but
C      leave it commented out in case we would implement an aternative scheme.
c      IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
#endif /* ALLOW_OBCS */

C-- end if exactConserv and not myIter=nIter0 block
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

# ifdef NONLIN_FRSURF
      IF (select_rStar .NE. 0) THEN
#  ifndef DISABLE_RSTAR_CODE
C-- note: rStarDhDt is similar to rStarDhCDt from S/R CALC_R_STAR
C         except for deep-factor and rho factor.
        DO j=1,sNy
         DO i=1,sNx
           ks = kSurfC(i,j,bi,bj)
           rStarDhDt(i,j) = 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,
     I                       addMass, rStarDhDt,
     O                       wVel,
     I                       myThid )

#ifdef EXACT_CONSERV
       IF ( k.EQ.Nr .AND. myIter.NE.0 .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--    reset W to zero outside the OB interior region
       IF ( useOBCS ) THEN
         DO j=1,sNy
          DO i=1,sNx
             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
          ENDDO
         ENDDO
       ENDIF
C--    Apply OBC to W (non-hydrostatic):
       IF ( useOBCS.AND.nonHydrostatic )
     &                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-|--+----|

C--   End bi,bj loops
       ENDDO
      ENDDO

      IF ( exactConserv .AND. myIter.NE.nIter0
     &                  .AND. implicDiv2Dflow .NE. 0. _d 0 )
     &    _EXCH_XY_RL( etaN , myThid )
      IF ( implicitIntGravWave .OR. myIter.EQ.nIter0 )
     &    _EXCH_XYZ_RL( wVel , myThid )

#ifdef EXACT_CONSERV
C--   Update etaH (from etaN and dEtaHdt):
      IF (exactConserv) THEN
c       IF ( useRealFreshWaterFlux .AND. myIter.EQ.nIter0 )
c    &    _EXCH_XY_RL( PmEpR, myThid )
#ifdef ALLOW_DEBUG
        IF (debugMode) CALL DEBUG_CALL('UPDATE_ETAH',myThid)
#endif
        CALL UPDATE_ETAH( myTime, myIter, myThid )
      ENDIF

#ifdef NONLIN_FRSURF
# ifndef DISABLE_SIGMA_CODE
      IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN
        CALL EXCH_XY_RL( dEtaHdt, myThid )
        CALL UPDATE_ETAWS( myTime, myIter, myThid )
      ENDIF
# endif /* DISABLE_SIGMA_CODE */
#endif /* NONLIN_FRSURF */

#endif /* EXACT_CONSERV */

      RETURN
      END