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	C $Header: /u/gcmpack/MITgcm/model/src/pre_cg3d.F,v 1.5 2012/04/11 15:49:38 jmc Exp $
2	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	C wSurfP2d :: surface vertical velocity after 2-D solver
57	INTEGER i,j,k,bi,bj
58	INTEGER ks, kp1
59	c CHARACTER*10 sufx
60	c CHARACTER*(MAX_LEN_MBUF) msgBuf
61	_RL locGamma, surfFac, tmpFac
62	_RL wFacKm, wFacKp
63	_RL uf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64	_RL vf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65	_RL wSurfP2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66	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	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	C (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
93	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	ENDDO
113	ENDIF
114	ENDIF
115
116	C-- Update or Add free-surface contribution to cg3d_b:
117	surfFac = 0.
118	IF ( selectNHfreeSurf.GE.1 ) THEN
119	tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
120	DO j=1,sNy
121	DO i=1,sNx
122	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	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
129	& + tmpFac*( wSurfP2d(i,j)
130	& + locGamma*wVel(i,j,ks,bi,bj) )
131	& /( 1. _d 0 + locGamma )
132	& _rA(i,j,bi,bj)deepFac2F(ks)
133	c ENDIF
134	C- Save wSurfP2d (used in POST_CG3D) into dPhiNH :
135	dPhiNH(i,j,bi,bj) = wSurfP2d(i,j)
136	ENDDO
137	ENDDO
138	ELSEIF ( .NOT.oldFreeSurfTerm ) THEN
139	tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
140	DO j=1,sNy
141	DO i=1,sNx
142	ks = kSurfC(i,j,bi,bj)
143	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	ENDDO
149	ENDDO
150	ELSEIF ( uniformFreeSurfLev ) THEN
151	C- Z coordinate: assume surface @ level k=1
152	surfFac = freeSurfFac*deepFac2F(1)
153	ELSE
154	C- Other than Z coordinate: no assumption on surface level index
155	DO j=1,sNy
156	DO i=1,sNx
157	ks = kSurfC(i,j,bi,bj)
158	IF ( ks.LE.Nr ) THEN
159	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
160	& +freeSurfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
161	& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTMom
162	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	#ifdef ALLOW_OBCS
177	& maskInC(i,j,bi,bj)maskInC(i-1,j,bi,bj)
178	#endif
179	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	#ifdef ALLOW_OBCS
183	& maskInC(i,j,bi,bj)maskInC(i,j-1,bi,bj)
184	#endif
185	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	& +( surfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
203	& -wVel(i,j,kp1,bi,bj)*wFacKp
204	& )*_rA(i,j,bi,bj)/deltaTMom
205	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	& )*_rA(i,j,bi,bj)/deltaTMom
224	ENDDO
225	ENDDO
226	ENDDO
227
228	#ifdef ALLOW_OBCS
229	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	IF (useOBCS) THEN
239	DO k=1,Nr
240	DO j=1,sNy
241	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	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