model/src/integr_continuity.F

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------------------------------------

#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. 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--    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 */

      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--    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-|--+----|

      RETURN
      END
1	jmc	1.23	C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.22 2009/11/28 20:59:18 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			I bi, bj, uFld, vFld,
12			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			C \| and free surface "r-anomaly" (etaN) to satisfy
18			C \| exactly the convervation 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			C uFld :: Zonal velocity ( m/s )
36			C vFld :: Meridional velocity ( m/s )
37			C bi,bj :: tile index
38			C myTime :: Current time in simulation
39			C myIter :: Current iteration number in simulation
40			C myThid :: Thread number for this instance of the routine.
41			_RL myTime
42			INTEGER myIter
43			INTEGER myThid
44			INTEGER bi,bj
45			_RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46	jmc	1.17	_RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
47	jmc	1.1
48			C !LOCAL VARIABLES:
49			C Local variables in common block
50
51			C Local variables
52			C i,j,k :: Loop counters
53			C uTrans :: Volume transports ( uVel.xA )
54			C vTrans :: Volume transports ( vVel.yA )
55	jmc	1.6	C hDivFlow :: Div. Barotropic Flow [transport unit m3/s]
56	jmc	1.19	INTEGER k
57			#ifdef EXACT_CONSERV
58			INTEGER i,j, 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.1	CEOP
68
69			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
70
71			#ifdef EXACT_CONSERV
72			IF (exactConserv) THEN
73
74	jmc	1.21	facEmP = 0.
75			IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit
76			facMass = 0.
77			IF ( selectAddFluid.GE.1 ) facMass = mass2rUnit
78
79	jmc	1.1	C-- Compute the Divergence of The Barotropic Flow :
80
81	jmc	1.17	C- Initialise
82	jmc	1.21	DO j=1-Oly,sNy+Oly
83			DO i=1-Olx,sNx+Olx
84	jmc	1.6	hDivFlow(i,j) = 0. _d 0
85			utrans(i,j) = 0. _d 0
86			vtrans(i,j) = 0. _d 0
87	jmc	1.21	ENDDO
88	jmc	1.1	ENDDO
89
90	jmc	1.21	DO k=1,Nr
91	jmc	1.17
92	jmc	1.1	C- Calculate velocity field "volume transports" through tracer cell faces
93	jmc	1.18	C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport).
94	jmc	1.1	DO j=1,sNy+1
95			DO i=1,sNx+1
96			uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
97	jmc	1.18	& deepFacC(k)rhoFacC(k)
98	jmc	1.1	& drF(k)_hFacW(i,j,k,bi,bj)
99			vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
100	jmc	1.18	& deepFacC(k)rhoFacC(k)
101	jmc	1.1	& drF(k)_hFacS(i,j,k,bi,bj)
102			ENDDO
103			ENDDO
104
105	jmc	1.17	C- Integrate vertically the Horizontal Divergence
106	jmc	1.1	DO j=1,sNy
107			DO i=1,sNx
108	jmc	1.6	hDivFlow(i,j) = hDivFlow(i,j)
109	jmc	1.1	& +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
110			& +vTrans(i,j+1)-vTrans(i,j) )
111	jmc	1.21	#ifdef ALLOW_ADDFLUID
112			& -facMass*addMass(i,j,k,bi,bj)
113			#endif /* ALLOW_ADDFLUID */
114	jmc	1.1	ENDDO
115			ENDDO
116
117			C- End DO k=1,Nr
118	jmc	1.21	ENDDO
119	jmc	1.1
120	jmc	1.6	C------------------------------------
121	jmc	1.18	C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ;
122			C anelastic: always assumes that rhoFacF(1) = 1
123	jmc	1.12	IF ( myTime.EQ.startTime .AND. myTime.NE.baseTime
124	jmc	1.21	& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
125	jmc	1.7
126	jmc	1.6	C needs previous time-step value of E-P-R, that has not been loaded
127	jmc	1.7	C and was not in old pickup-file ; try to use etaN & etaH instead.
128			IF ( usePickupBeforeC54 ) THEN
129	jmc	1.6	DO j=1,sNy
130			DO i=1,sNx
131	jmc	1.7	dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
132			& / (implicDiv2Dflow*deltaTfreesurf)
133			ENDDO
134			ENDDO
135			ENDIF
136
137			DO j=1,sNy
138			DO i=1,sNx
139	jmc	1.6	PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
140			& + hDivFlow(i,j)*recip_rA(i,j,bi,bj)
141	jmc	1.18	& *recip_deepFac2F(1)
142	jmc	1.20	PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass
143	jmc	1.6	ENDDO
144	jmc	1.7	ENDDO
145			ELSEIF ( myTime.EQ.startTime ) THEN
146	jmc	1.6	DO j=1,sNy
147			DO i=1,sNx
148	jmc	1.23	ks = kSurfC(I,J,bi,bj)
149	jmc	1.6	PmEpR(i,j,bi,bj) = 0. _d 0
150			dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
151	jmc	1.18	& *recip_deepFac2F(ks)
152	jmc	1.6	ENDDO
153	jmc	1.17	ENDDO
154	jmc	1.6	ELSE
155	jmc	1.17	C-- Needs to get valid PmEpR (for T,S forcing) also in overlap regions
156			C (e.g., if using KPP) => set over full index range
157			DO j=1-OLy,sNy+OLy
158			DO i=1-OLx,sNx+OLx
159			PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
160			ENDDO
161			ENDDO
162	jmc	1.6	DO j=1,sNy
163			DO i=1,sNx
164	jmc	1.23	ks = kSurfC(i,j,bi,bj)
165	jmc	1.6	dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj)
166	jmc	1.18	& *recip_deepFac2F(ks)
167	jmc	1.6	& -facEmP*EmPmR(i,j,bi,bj)
168			ENDDO
169			ENDDO
170			ENDIF
171			C------------------------------------
172
173	jmc	1.23	#ifdef ALLOW_OBCS
174			C-- reset dEtaHdt to zero outside the OB interior region
175			IF ( useOBCS ) THEN
176			DO j=1,sNy
177			DO i=1,sNx
178			dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj)
179			ENDDO
180			ENDDO
181			ENDIF
182			#endif /* ALLOW_OBCS */
183
184			C-- end if exactConserv block
185	jmc	1.1	ENDIF
186
187			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
188
189			IF (exactConserv .AND. myTime.NE.startTime) THEN
190	jmc	1.17	C-- Update etaN at the end of the time step :
191	jmc	1.1	C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
192
193			IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
194			DO j=1-Oly,sNy+Oly
195			DO i=1-Olx,sNx+Olx
196	jmc	1.17	etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
197	jmc	1.1	ENDDO
198			ENDDO
199			ELSE
200			DO j=1,sNy
201			DO i=1,sNx
202			etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
203	jmc	1.6	& + implicDiv2DflowdEtaHdt(i,j,bi,bj)deltaTfreesurf
204	jmc	1.1	ENDDO
205			ENDDO
206			ENDIF
207
208	dimitri	1.11	#ifdef ALLOW_OBCS
209	jmc	1.23	C-- Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer
210			C needed with proper masking in solver (matrix+cg2d_b,x) and masking
211			C of dEtaHdt above. etaN next to OB does not enter present algorithm but
212			C leave it commented out in case we would implement an aternative scheme.
213			c IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid )
214	jmc	1.17	#endif /* ALLOW_OBCS */
215	adcroft	1.10
216	jmc	1.1	ENDIF
217	jmc	1.3
218			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
219
220	heimbach	1.16	# ifdef NONLIN_FRSURF
221	jmc	1.3	IF (select_rStar .NE. 0) THEN
222	heimbach	1.16	# ifndef DISABLE_RSTAR_CODE
223	jmc	1.22	C-- note: final value of rStarDhCDt from S/R CALC_R_STAR
224			C will be different (no deep-factor, no rho factor)
225	jmc	1.6	DO j=1,sNy
226			DO i=1,sNx
227	jmc	1.23	ks = kSurfC(i,j,bi,bj)
228	jmc	1.22	rStarDhCDt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)
229			& deepFac2F(ks)rhoFacF(ks)
230			& *recip_Rcol(i,j,bi,bj)
231	jmc	1.3	ENDDO
232	heimbach	1.16	ENDDO
233			# endif /* DISABLE_RSTAR_CODE */
234	jmc	1.3	ENDIF
235	heimbach	1.16	# endif /* NONLIN_FRSURF */
236
237	jmc	1.4	#endif /* EXACT_CONSERV */
238	jmc	1.1
239			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
240
241			DO k=Nr,1,-1
242			C-- Integrate continuity vertically for vertical velocity
243
244			CALL INTEGRATE_FOR_W(
245	jmc	1.21	I bi, bj, k, uFld, vFld, addMass,
246	jmc	1.1	O wVel,
247			I myThid )
248	jmc	1.17
249	jmc	1.8	#ifdef EXACT_CONSERV
250	jmc	1.21	IF ( k.EQ.Nr .AND. myTime.NE.baseTime .AND. usingPCoords
251			& .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN
252	jmc	1.2
253	jmc	1.23	DO j=1,sNy
254			DO i=1,sNx
255			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
256	jmc	1.20	& +mass2rUnitPmEpR(i,j,bi,bj)maskC(i,j,k,bi,bj)
257	jmc	1.2	ENDDO
258	jmc	1.23	ENDDO
259	jmc	1.2
260			ENDIF
261	jmc	1.8	#endif /* EXACT_CONSERV */
262	jmc	1.2
263	jmc	1.1	#ifdef ALLOW_OBCS
264	jmc	1.23	C-- reset W to zero outside the OB interior region
265			IF ( useOBCS ) THEN
266			DO j=1,sNy
267			DO i=1,sNx
268			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj)
269			ENDDO
270			ENDDO
271			ENDIF
272			C-- Apply OBC to W (non-hydrostatic):
273			IF ( useOBCS.AND.nonHydrostatic )
274			& CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
275	jmc	1.17	#endif /* ALLOW_OBCS */
276	jmc	1.1
277			C- End DO k=Nr,1,-1
278			ENDDO
279
280			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
281
282			RETURN
283			END