model/src/pre_cg3d.F

C $Header: /u/gcmpack/MITgcm/model/src/pre_cg3d.F,v 1.5 2012/04/11 15:49:38 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: PRE_CG3D
C     !INTERFACE:
      SUBROUTINE PRE_CG3D(
     I                     oldFreeSurfTerm,
     I                     cg2d_x,
     U                     cg3d_b,
     I                     myTime, myIter, myThid )

C     !DESCRIPTION:
C     Called from SOLVE_FOR_PRESSURE, before 3-D solver (cg3d):
C     Finish calculation of 3-D RHS after 2-D inversionis done.

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

C     === Functions ====
c     LOGICAL  DIFFERENT_MULTIPLE
c     EXTERNAL DIFFERENT_MULTIPLE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     oldFreeSurfTerm :: Treat free-surface term in the old way (no exactConserv)
C     cg2d_x          :: Solution vector of the 2-D solver equation a.x=b
C     cg3d_b          :: Right Hand side vector of the 3-D solver equation A.X=B
C     myTime          :: Current time in simulation
C     myIter          :: Current iteration number in simulation
C     myThid          :: My Thread Id number
      LOGICAL oldFreeSurfTerm
      _RL     cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
C     wSurfP2d :: surface vertical velocity after 2-D solver
      INTEGER i,j,k,bi,bj
      INTEGER ks, kp1
c     CHARACTER*10 sufx
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL     locGamma, surfFac, tmpFac
      _RL     wFacKm, wFacKp
      _RL     uf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     vf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     wSurfP2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

c     IF ( use3Dsolver ) THEN

C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
C     see CG3D.h for the interface to this routine.
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

C--   Calculate updated (after 2-D solver) vertical velocity at the surface
         IF ( oldFreeSurfTerm .OR. implicDiv2DFlow.EQ.zeroRL ) THEN
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              wSurfP2d(i,j) = 0. _d 0
            ENDDO
           ENDDO
         ELSE
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
               wSurfP2d(i,j) = ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) )
     &                       / ( implicDiv2DFlow*deltaTFreeSurf )
            ENDDO
           ENDDO
         ENDIF

C--   Add EmPmR contribution to top level cg3d_b or to wSurfP2d:
C      (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
         IF ( useRealFreshWaterFlux.AND.fluidIsWater ) THEN
          IF ( oldFreeSurfTerm .OR. usingPCoords ) THEN
           tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow/deltaTMom
           ks = 1
           IF ( usingPCoords ) ks = Nr
           DO j=1,sNy
            DO i=1,sNx
              cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &          + tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
     &                                 *maskInC(i,j,bi,bj)
            ENDDO
           ENDDO
          ELSE
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              wSurfP2d(i,j) = wSurfP2d(i,j)
     &                      + EmPmR(i,j,bi,bj)*mass2rUnit
     &                                 *maskInC(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDIF
         ENDIF

C--   Update or Add free-surface contribution to cg3d_b:
         surfFac = 0.
         IF ( selectNHfreeSurf.GE.1 ) THEN
           tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
           DO j=1,sNy
            DO i=1,sNx
              locGamma = drC(1)*recip_Bo(i,j,bi,bj)
     &                 /( deltaTMom*deltaTFreeSurf
     &                   *implicitNHPress*implicDiv2DFlow )
              ks = 1
c             ks = kSurfC(i,j,bi,bj)
c             IF ( ks.LE.Nr ) THEN
               cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &                      + tmpFac*( wSurfP2d(i,j)
     &                               + locGamma*wVel(i,j,ks,bi,bj) )
     &                              /( 1. _d 0 + locGamma )
     &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
c             ENDIF
C-    Save wSurfP2d (used in POST_CG3D) into dPhiNH :
              dPhiNH(i,j,bi,bj) = wSurfP2d(i,j)
            ENDDO
           ENDDO
         ELSEIF ( .NOT.oldFreeSurfTerm ) THEN
           tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
           DO j=1,sNy
            DO i=1,sNx
              ks = kSurfC(i,j,bi,bj)
              IF ( ks.LE.Nr ) THEN
               cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &                      + tmpFac*wSurfP2d(i,j)
     &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
              ENDIF
            ENDDO
           ENDDO
         ELSEIF ( uniformFreeSurfLev ) THEN
C-       Z coordinate: assume surface @ level k=1
           surfFac = freeSurfFac*deepFac2F(1)
         ELSE
C-       Other than Z coordinate: no assumption on surface level index
           DO j=1,sNy
            DO i=1,sNx
              ks = kSurfC(i,j,bi,bj)
              IF ( ks.LE.Nr ) THEN
               cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
     &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom
              ENDIF
            ENDDO
           ENDDO
         ENDIF

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

C--   Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
C                             2) add vertical velocity contribution.
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j) = -_recip_dxC(i,j,bi,bj)
     &             * implicSurfPress*implicDiv2DFlow
     &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
#ifdef ALLOW_OBCS
     &             *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
#endif
           vf(i,j) = -_recip_dyC(i,j,bi,bj)
     &             * implicSurfPress*implicDiv2DFlow
     &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
#ifdef ALLOW_OBCS
     &             *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
#endif
          ENDDO
         ENDDO

C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
C       uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
C       always recomputed from continuity eq (like eta when exactConserv=T)
         k=1
         kp1 = MIN(k+1,Nr)
         wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
         IF (k.GE.Nr) wFacKp = 0.
         DO j=1,sNy
          DO i=1,sNx
            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
     &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
     &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
     &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
     &       +( surfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
     &         -wVel(i,j,kp1,bi,bj)*wFacKp
     &        )*_rA(i,j,bi,bj)/deltaTMom
          ENDDO
         ENDDO
         DO k=2,Nr
          kp1 = MIN(k+1,Nr)
C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
C        both appear in wVel term, but at 2 different levels
          wFacKm = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k )
          wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
          IF (k.GE.Nr) wFacKp = 0.
          DO j=1,sNy
           DO i=1,sNx
            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
     &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
     &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
     &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
     &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
     &         -wVel(i,j,kp1,bi,bj)*wFacKp
     &        )*_rA(i,j,bi,bj)/deltaTMom
           ENDDO
          ENDDO
         ENDDO

#ifdef ALLOW_OBCS
C- Note: solver matrix is trivial outside OB region (main diagonal only)
C     => no real need to reset RHS (=cg3d_b) & cg3d_x, except that:
C    a) normalisation is fct of Max(RHS), which can be large ouside OB region
C      (would be different if we were solving for increment of phi_nh
C       instead of directly for phi_nh).
C       => need to reset RHS to ensure that interior solution does not depend
C       on ouside OB region.
C    b) provide directly the trivial solution cg3d_x == 0 for outside OB region
C      (=> no residual => no effect on solver convergence and interior solution)
         IF (useOBCS) THEN
          DO k=1,Nr
           DO j=1,sNy
            DO i=1,sNx
              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &                            *maskInC(i,j,bi,bj)
              phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
     &                            *maskInC(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDIF
#endif /* ALLOW_OBCS */

C-    end bi,bj loops
        ENDDO
       ENDDO

c     ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/model/src/pre_cg3d.F,v 1.5 2012/04/11 15:49:38 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: PRE_CG3D
9			C !INTERFACE:
10			SUBROUTINE PRE_CG3D(
11			I oldFreeSurfTerm,
12			I cg2d_x,
13			U cg3d_b,
14			I myTime, myIter, myThid )
15
16			C !DESCRIPTION:
17			C Called from SOLVE_FOR_PRESSURE, before 3-D solver (cg3d):
18			C Finish calculation of 3-D RHS after 2-D inversionis done.
19
20			C !USES:
21			IMPLICIT NONE
22			C == Global variables
23			#include "SIZE.h"
24			#include "EEPARAMS.h"
25			#include "PARAMS.h"
26			#include "GRID.h"
27			#include "SURFACE.h"
28			#include "FFIELDS.h"
29			#include "DYNVARS.h"
30			#ifdef ALLOW_NONHYDROSTATIC
31			#include "NH_VARS.h"
32			#endif
33
34			C === Functions ====
35			c LOGICAL DIFFERENT_MULTIPLE
36			c EXTERNAL DIFFERENT_MULTIPLE
37
38			C !INPUT/OUTPUT PARAMETERS:
39			C == Routine arguments ==
40			C oldFreeSurfTerm :: Treat free-surface term in the old way (no exactConserv)
41			C cg2d_x :: Solution vector of the 2-D solver equation a.x=b
42			C cg3d_b :: Right Hand side vector of the 3-D solver equation A.X=B
43			C myTime :: Current time in simulation
44			C myIter :: Current iteration number in simulation
45			C myThid :: My Thread Id number
46			LOGICAL oldFreeSurfTerm
47			_RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48			_RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49			_RL myTime
50			INTEGER myIter
51			INTEGER myThid
52
53			#ifdef ALLOW_NONHYDROSTATIC
54			C !LOCAL VARIABLES:
55			C == Local variables ==
56	jmc	1.6	C wSurfP2d :: surface vertical velocity after 2-D solver
57	jmc	1.1	INTEGER i,j,k,bi,bj
58			INTEGER ks, kp1
59			c CHARACTER*10 sufx
60			c CHARACTER*(MAX_LEN_MBUF) msgBuf
61	jmc	1.6	_RL locGamma, surfFac, tmpFac
62	jmc	1.1	_RL wFacKm, wFacKp
63	jmc	1.4	_RL uf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64			_RL vf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65	jmc	1.6	_RL wSurfP2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66	jmc	1.1	CEOP
67
68			c IF ( use3Dsolver ) THEN
69
70			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
71			C see CG3D.h for the interface to this routine.
72			DO bj=myByLo(myThid),myByHi(myThid)
73			DO bi=myBxLo(myThid),myBxHi(myThid)
74
75	jmc	1.6	C-- Calculate updated (after 2-D solver) vertical velocity at the surface
76			IF ( oldFreeSurfTerm .OR. implicDiv2DFlow.EQ.zeroRL ) THEN
77			DO j=1-OLy,sNy+OLy
78			DO i=1-OLx,sNx+OLx
79			wSurfP2d(i,j) = 0. _d 0
80			ENDDO
81			ENDDO
82			ELSE
83			DO j=1-OLy,sNy+OLy
84			DO i=1-OLx,sNx+OLx
85			wSurfP2d(i,j) = ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) )
86			& / ( implicDiv2DFlow*deltaTFreeSurf )
87			ENDDO
88			ENDDO
89			ENDIF
90
91			C-- Add EmPmR contribution to top level cg3d_b or to wSurfP2d:
92	jmc	1.1	C (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
93	jmc	1.6	IF ( useRealFreshWaterFlux.AND.fluidIsWater ) THEN
94			IF ( oldFreeSurfTerm .OR. usingPCoords ) THEN
95			tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow/deltaTMom
96			ks = 1
97			IF ( usingPCoords ) ks = Nr
98			DO j=1,sNy
99			DO i=1,sNx
100			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
101			& + tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)
102			& *maskInC(i,j,bi,bj)
103			ENDDO
104			ENDDO
105			ELSE
106			DO j=1-OLy,sNy+OLy
107			DO i=1-OLx,sNx+OLx
108			wSurfP2d(i,j) = wSurfP2d(i,j)
109			& + EmPmR(i,j,bi,bj)*mass2rUnit
110			& *maskInC(i,j,bi,bj)
111			ENDDO
112	jmc	1.1	ENDDO
113	jmc	1.6	ENDIF
114			ENDIF
115	jmc	1.1
116			C-- Update or Add free-surface contribution to cg3d_b:
117	jmc	1.6	surfFac = 0.
118			IF ( selectNHfreeSurf.GE.1 ) THEN
119			tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
120	jmc	1.1	DO j=1,sNy
121			DO i=1,sNx
122	jmc	1.6	locGamma = drC(1)*recip_Bo(i,j,bi,bj)
123			& /( deltaTMom*deltaTFreeSurf
124			& implicitNHPressimplicDiv2DFlow )
125			ks = 1
126			c ks = kSurfC(i,j,bi,bj)
127			c IF ( ks.LE.Nr ) THEN
128	jmc	1.1	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
129	jmc	1.6	& + tmpFac*( wSurfP2d(i,j)
130			& + locGamma*wVel(i,j,ks,bi,bj) )
131			& /( 1. _d 0 + locGamma )
132	jmc	1.1	& _rA(i,j,bi,bj)deepFac2F(ks)
133	jmc	1.6	c ENDIF
134			C- Save wSurfP2d (used in POST_CG3D) into dPhiNH :
135			dPhiNH(i,j,bi,bj) = wSurfP2d(i,j)
136	jmc	1.1	ENDDO
137			ENDDO
138	jmc	1.6	ELSEIF ( .NOT.oldFreeSurfTerm ) THEN
139			tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
140	jmc	1.1	DO j=1,sNy
141			DO i=1,sNx
142	jmc	1.4	ks = kSurfC(i,j,bi,bj)
143	jmc	1.6	IF ( ks.LE.Nr ) THEN
144			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
145			& + tmpFac*wSurfP2d(i,j)
146			& _rA(i,j,bi,bj)deepFac2F(ks)
147			ENDIF
148	jmc	1.1	ENDDO
149			ENDDO
150	jmc	1.5	ELSEIF ( uniformFreeSurfLev ) THEN
151	jmc	1.1	C- Z coordinate: assume surface @ level k=1
152	jmc	1.6	surfFac = freeSurfFac*deepFac2F(1)
153	jmc	1.1	ELSE
154			C- Other than Z coordinate: no assumption on surface level index
155			DO j=1,sNy
156			DO i=1,sNx
157	jmc	1.4	ks = kSurfC(i,j,bi,bj)
158	jmc	1.1	IF ( ks.LE.Nr ) THEN
159			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
160	jmc	1.6	& +freeSurfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
161			& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTMom
162	jmc	1.1	ENDIF
163			ENDDO
164			ENDDO
165			ENDIF
166
167			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
168
169			C-- Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
170			C 2) add vertical velocity contribution.
171			DO j=1,sNy+1
172			DO i=1,sNx+1
173			uf(i,j) = -_recip_dxC(i,j,bi,bj)
174			& * implicSurfPress*implicDiv2DFlow
175			& *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
176	jmc	1.3	#ifdef ALLOW_OBCS
177			& maskInC(i,j,bi,bj)maskInC(i-1,j,bi,bj)
178			#endif
179	jmc	1.1	vf(i,j) = -_recip_dyC(i,j,bi,bj)
180			& * implicSurfPress*implicDiv2DFlow
181			& *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
182	jmc	1.3	#ifdef ALLOW_OBCS
183			& maskInC(i,j,bi,bj)maskInC(i,j-1,bi,bj)
184			#endif
185	jmc	1.1	ENDDO
186			ENDDO
187
188			C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
189			C uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
190			C always recomputed from continuity eq (like eta when exactConserv=T)
191			k=1
192			kp1 = MIN(k+1,Nr)
193			wFacKp = implicDiv2DFlowdeepFac2F(kp1)rhoFacF(kp1)
194			IF (k.GE.Nr) wFacKp = 0.
195			DO j=1,sNy
196			DO i=1,sNx
197			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
198			& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
199			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
200			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
201			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
202	jmc	1.6	& +( surfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
203	jmc	1.1	& -wVel(i,j,kp1,bi,bj)*wFacKp
204	jmc	1.6	& )*_rA(i,j,bi,bj)/deltaTMom
205	jmc	1.1	ENDDO
206			ENDDO
207			DO k=2,Nr
208			kp1 = MIN(k+1,Nr)
209			C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
210			C both appear in wVel term, but at 2 different levels
211			wFacKm = implicDiv2DFlowdeepFac2F( k )rhoFacF( k )
212			wFacKp = implicDiv2DFlowdeepFac2F(kp1)rhoFacF(kp1)
213			IF (k.GE.Nr) wFacKp = 0.
214			DO j=1,sNy
215			DO i=1,sNx
216			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
217			& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
218			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
219			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
220			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
221			& +( wVel(i,j, k ,bi,bj)wFacKmmaskC(i,j,k-1,bi,bj)
222			& -wVel(i,j,kp1,bi,bj)*wFacKp
223	jmc	1.6	& )*_rA(i,j,bi,bj)/deltaTMom
224	jmc	1.1	ENDDO
225			ENDDO
226			ENDDO
227
228			#ifdef ALLOW_OBCS
229	jmc	1.3	C- Note: solver matrix is trivial outside OB region (main diagonal only)
230			C => no real need to reset RHS (=cg3d_b) & cg3d_x, except that:
231			C a) normalisation is fct of Max(RHS), which can be large ouside OB region
232			C (would be different if we were solving for increment of phi_nh
233			C instead of directly for phi_nh).
234			C => need to reset RHS to ensure that interior solution does not depend
235			C on ouside OB region.
236			C b) provide directly the trivial solution cg3d_x == 0 for outside OB region
237			C (=> no residual => no effect on solver convergence and interior solution)
238	jmc	1.1	IF (useOBCS) THEN
239			DO k=1,Nr
240			DO j=1,sNy
241	jmc	1.3	DO i=1,sNx
242			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
243			& *maskInC(i,j,bi,bj)
244			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
245			& *maskInC(i,j,bi,bj)
246			ENDDO
247	jmc	1.1	ENDDO
248			ENDDO
249			ENDIF
250			#endif /* ALLOW_OBCS */
251
252			C- end bi,bj loops
253			ENDDO
254			ENDDO
255
256			c ENDIF
257			#endif /* ALLOW_NONHYDROSTATIC */
258
259			RETURN
260			END