model/src/integr_continuity.F

C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.25 2011/05/23 00:39:20 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 bi,bj
      INTEGER i,j,k
#ifdef EXACT_CONSERV
      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
#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
#ifdef ALLOW_DEBUG
        IF (debugMode) CALL DEBUG_CALL('UPDATE_ETAH',myThid)
#endif
        CALL UPDATE_ETAH( myTime, myIter, myThid )
      ENDIF
#endif /* EXACT_CONSERV */

      RETURN
      END
1	jmc	1.26	C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.25 2011/05/23 00:39:20 jmc Exp $
2	jmc	1.1	C $Name: $
3
4	edhill	1.5	#include "PACKAGES_CONFIG.h"
5	jmc	1.1	#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: INTEGR_CONTINUITY
9			C !INTERFACE:
10			SUBROUTINE INTEGR_CONTINUITY(
11	jmc	1.26	I uFld, vFld,
12	jmc	1.1	I myTime, myIter, myThid )
13			C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.17	C \| SUBROUTINE INTEGR_CONTINUITY
16			C \| o Integrate the continuity Eq : compute vertical velocity
17	jmc	1.26	C \| and free surface "r-anomaly" (etaN,etaH) to satisfy
18			C \| exactly the conservation of the Total Volume
19	jmc	1.1	C ==========================================================
20			C \ev
21
22			C !USES:
23			IMPLICIT NONE
24			C == Global variables
25			#include "SIZE.h"
26			#include "EEPARAMS.h"
27			#include "PARAMS.h"
28			#include "DYNVARS.h"
29			#include "GRID.h"
30			#include "SURFACE.h"
31			#include "FFIELDS.h"
32
33			C !INPUT/OUTPUT PARAMETERS:
34			C == Routine arguments ==
35	jmc	1.26	C uFld :: Zonal velocity ( m/s )
36			C vFld :: Meridional velocity ( m/s )
37			C myTime :: Current time in simulation
38			C myIter :: Current iteration number in simulation
39			C myThid :: my Thread Id. number
40	jmc	1.1	_RL myTime
41			INTEGER myIter
42			INTEGER myThid
43			_RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
44	jmc	1.17	_RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
45	jmc	1.1
46			C !LOCAL VARIABLES:
47			C Local variables in common block
48
49			C Local variables
50	jmc	1.26	C bi,bj :: tile index
51			C i,j,k :: Loop counters
52			C uTrans :: Volume transports ( uVel.xA )
53			C vTrans :: Volume transports ( vVel.yA )
54	jmc	1.6	C hDivFlow :: Div. Barotropic Flow [transport unit m3/s]
55	jmc	1.26	INTEGER bi,bj
56	mmazloff	1.24	INTEGER i,j,k
57	jmc	1.19	#ifdef EXACT_CONSERV
58	mmazloff	1.24	INTEGER ks
59	jmc	1.1	_RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60			_RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61	jmc	1.6	_RL hDivFlow(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62	jmc	1.21	_RL facEmP, facMass
63			#endif /* EXACT_CONSERV */
64			#ifndef ALLOW_ADDFLUID
65			_RL addMass(1)
66			#endif /* ndef ALLOW_ADDFLUID */
67	jmc	1.25	#if (defined NONLIN_FRSURF) && !(defined DISABLE_RSTAR_CODE)
68			_RL rStarDhDt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69			#else
70			_RL rStarDhDt(1)
71			#endif
72	jmc	1.1	CEOP
73
74	jmc	1.26	C-- Start bi,bj loops
75			DO bj=myByLo(myThid),myByHi(myThid)
76			DO bi=myBxLo(myThid),myBxHi(myThid)
77
78	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
79
80			#ifdef EXACT_CONSERV
81			IF (exactConserv) THEN
82
83	jmc	1.21	facEmP = 0.
84			IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit
85			facMass = 0.
86			IF ( selectAddFluid.GE.1 ) facMass = mass2rUnit
87
88	jmc	1.1	C-- Compute the Divergence of The Barotropic Flow :
89
90	jmc	1.17	C- Initialise
91	jmc	1.21	DO j=1-Oly,sNy+Oly
92			DO i=1-Olx,sNx+Olx
93	jmc	1.6	hDivFlow(i,j) = 0. _d 0
94			utrans(i,j) = 0. _d 0
95			vtrans(i,j) = 0. _d 0
96	jmc	1.21	ENDDO
97	jmc	1.1	ENDDO
98
99	jmc	1.21	DO k=1,Nr
100	jmc	1.17
101	jmc	1.1	C- Calculate velocity field "volume transports" through tracer cell faces
102	jmc	1.18	C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport).
103	jmc	1.1	DO j=1,sNy+1
104			DO i=1,sNx+1
105			uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
106	jmc	1.18	& deepFacC(k)rhoFacC(k)
107	jmc	1.1	& drF(k)_hFacW(i,j,k,bi,bj)
108			vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
109	jmc	1.18	& deepFacC(k)rhoFacC(k)
110	jmc	1.1	& drF(k)_hFacS(i,j,k,bi,bj)
111			ENDDO
112			ENDDO
113
114	jmc	1.17	C- Integrate vertically the Horizontal Divergence
115	jmc	1.1	DO j=1,sNy
116			DO i=1,sNx
117	jmc	1.6	hDivFlow(i,j) = hDivFlow(i,j)
118	jmc	1.1	& +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
119			& +vTrans(i,j+1)-vTrans(i,j) )
120	jmc	1.21	#ifdef ALLOW_ADDFLUID
121			& -facMass*addMass(i,j,k,bi,bj)
122			#endif /* ALLOW_ADDFLUID */
123	jmc	1.1	ENDDO
124			ENDDO
125
126			C- End DO k=1,Nr
127	jmc	1.21	ENDDO
128	jmc	1.1
129	jmc	1.6	C------------------------------------
130	jmc	1.18	C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ;
131			C anelastic: always assumes that rhoFacF(1) = 1
132	jmc	1.26	IF ( myIter.EQ.nIter0 .AND. myIter.NE.0
133			& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
134	jmc	1.7
135	jmc	1.6	C needs previous time-step value of E-P-R, that has not been loaded
136	jmc	1.7	C and was not in old pickup-file ; try to use etaN & etaH instead.
137			IF ( usePickupBeforeC54 ) THEN
138	jmc	1.6	DO j=1,sNy
139			DO i=1,sNx
140	jmc	1.7	dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
141			& / (implicDiv2Dflow*deltaTfreesurf)
142			ENDDO
143			ENDDO
144			ENDIF
145
146			DO j=1,sNy
147			DO i=1,sNx
148	jmc	1.6	PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
149			& + hDivFlow(i,j)*recip_rA(i,j,bi,bj)
150	jmc	1.18	& *recip_deepFac2F(1)
151	jmc	1.20	PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass
152	jmc	1.6	ENDDO
153	jmc	1.7	ENDDO
154	jmc	1.26	ELSEIF ( myIter.EQ.nIter0 ) THEN
155	jmc	1.6	DO j=1,sNy
156			DO i=1,sNx
157	jmc	1.23	ks = kSurfC(I,J,bi,bj)
158	jmc	1.6	PmEpR(i,j,bi,bj) = 0. _d 0
159			dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
160	jmc	1.18	& *recip_deepFac2F(ks)
161	jmc	1.6	ENDDO
162	jmc	1.17	ENDDO
163	jmc	1.6	ELSE
164	jmc	1.17	C-- Needs to get valid PmEpR (for T,S forcing) also in overlap regions
165			C (e.g., if using KPP) => set over full index range
166			DO j=1-OLy,sNy+OLy
167			DO i=1-OLx,sNx+OLx
168			PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
169			ENDDO
170			ENDDO
171	jmc	1.6	DO j=1,sNy
172			DO i=1,sNx
173	jmc	1.23	ks = kSurfC(i,j,bi,bj)
174	jmc	1.6	dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
175	jmc	1.18	& *recip_deepFac2F(ks)
176	jmc	1.6	& -facEmP*EmPmR(i,j,bi,bj)
177			ENDDO
178			ENDDO
179			ENDIF
180			C------------------------------------
181
182	jmc	1.23	#ifdef ALLOW_OBCS
183			C-- reset dEtaHdt to zero outside the OB interior region
184			IF ( useOBCS ) THEN
185			DO j=1,sNy
186			DO i=1,sNx
187			dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj)
188			ENDDO
189			ENDDO
190			ENDIF
191			#endif /* ALLOW_OBCS */
192
193			C-- end if exactConserv block
194	jmc	1.1	ENDIF
195
196			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
197
198	jmc	1.26	IF ( exactConserv .AND. myIter.NE.nIter0 ) THEN
199	jmc	1.17	C-- Update etaN at the end of the time step :
200	jmc	1.1	C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
201
202			IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
203			DO j=1-Oly,sNy+Oly
204			DO i=1-Olx,sNx+Olx
205	jmc	1.17	etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
206	jmc	1.1	ENDDO
207			ENDDO
208			ELSE
209			DO j=1,sNy
210			DO i=1,sNx
211			etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
212	jmc	1.6	& + implicDiv2DflowdEtaHdt(i,j,bi,bj)deltaTfreesurf
213	jmc	1.1	ENDDO
214			ENDDO
215			ENDIF
216
217	dimitri	1.11	#ifdef ALLOW_OBCS
218	jmc	1.23	C-- Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer
219			C needed with proper masking in solver (matrix+cg2d_b,x) and masking
220			C of dEtaHdt above. etaN next to OB does not enter present algorithm but
221			C leave it commented out in case we would implement an aternative scheme.
222			c IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
223	jmc	1.17	#endif /* ALLOW_OBCS */
224	adcroft	1.10
225	jmc	1.26	C-- end if exactConserv and not myIter=nIter0 block
226	jmc	1.1	ENDIF
227	jmc	1.3
228			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
229
230	heimbach	1.16	# ifdef NONLIN_FRSURF
231	jmc	1.3	IF (select_rStar .NE. 0) THEN
232	heimbach	1.16	# ifndef DISABLE_RSTAR_CODE
233	jmc	1.25	C-- note: rStarDhDt is similar to rStarDhCDt from S/R CALC_R_STAR
234			C except for deep-factor and rho factor.
235	jmc	1.6	DO j=1,sNy
236			DO i=1,sNx
237	jmc	1.23	ks = kSurfC(i,j,bi,bj)
238	jmc	1.25	rStarDhDt(i,j) = dEtaHdt(i,j,bi,bj)
239			& deepFac2F(ks)rhoFacF(ks)
240			& *recip_Rcol(i,j,bi,bj)
241	jmc	1.3	ENDDO
242	heimbach	1.16	ENDDO
243			# endif /* DISABLE_RSTAR_CODE */
244	jmc	1.3	ENDIF
245	heimbach	1.16	# endif /* NONLIN_FRSURF */
246
247	jmc	1.4	#endif /* EXACT_CONSERV */
248	jmc	1.1
249			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
250
251			DO k=Nr,1,-1
252			C-- Integrate continuity vertically for vertical velocity
253
254			CALL INTEGRATE_FOR_W(
255	jmc	1.25	I bi, bj, k, uFld, vFld,
256			I addMass, rStarDhDt,
257	jmc	1.1	O wVel,
258			I myThid )
259	jmc	1.17
260	jmc	1.8	#ifdef EXACT_CONSERV
261	jmc	1.26	IF ( k.EQ.Nr .AND. myIter.NE.0 .AND. usingPCoords
262			& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
263	jmc	1.2
264	jmc	1.23	DO j=1,sNy
265			DO i=1,sNx
266			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
267	jmc	1.20	& +mass2rUnitPmEpR(i,j,bi,bj)maskC(i,j,k,bi,bj)
268	jmc	1.2	ENDDO
269	jmc	1.23	ENDDO
270	jmc	1.2
271			ENDIF
272	jmc	1.8	#endif /* EXACT_CONSERV */
273	jmc	1.2
274	jmc	1.1	#ifdef ALLOW_OBCS
275	jmc	1.23	C-- reset W to zero outside the OB interior region
276			IF ( useOBCS ) THEN
277			DO j=1,sNy
278			DO i=1,sNx
279			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
280			ENDDO
281			ENDDO
282			ENDIF
283			C-- Apply OBC to W (non-hydrostatic):
284			IF ( useOBCS.AND.nonHydrostatic )
285			& CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
286	jmc	1.17	#endif /* ALLOW_OBCS */
287	jmc	1.1
288			C- End DO k=Nr,1,-1
289			ENDDO
290
291			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
292
293	jmc	1.26	C-- End bi,bj loops
294			ENDDO
295			ENDDO
296
297			IF ( exactConserv .AND. myIter.NE.nIter0
298			& .AND. implicDiv2Dflow .NE. 0. _d 0 )
299			& _EXCH_XY_RL( etaN , myThid )
300			IF ( implicitIntGravWave .OR. myIter.EQ.nIter0 )
301			& _EXCH_XYZ_RL( wVel , myThid )
302
303			#ifdef EXACT_CONSERV
304			C-- Update etaH (from etaN and dEtaHdt):
305			IF (exactConserv) THEN
306			#ifdef ALLOW_DEBUG
307			IF (debugMode) CALL DEBUG_CALL('UPDATE_ETAH',myThid)
308			#endif
309			CALL UPDATE_ETAH( myTime, myIter, myThid )
310			ENDIF
311			#endif /* EXACT_CONSERV */
312
313	jmc	1.1	RETURN
314			END