model/src/solve_for_pressure.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.19 2001/03/09 20:45:09 adcroft Exp $
C $Name: checkpoint37 $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SOLVE_FOR_PRESSURE( myThid )
C     /==========================================================\
C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |
C     | o Controls inversion of two and/or three-dimensional     |
C     |   elliptic problems for the pressure field.              |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "CG3D.h"
#include "GW.h"
#endif
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif

C     == Routine arguments ==
C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
      INTEGER myThid
CEndOfInterface

C     Local variables
C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector
C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
      INTEGER i,j,k,bi,bj
      _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL tolerance,residual
      INTEGER numIters

C--   Save previous solution & Initialise Vector solution and source term :
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#ifdef INCLUDE_CD_CODE
          etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
#endif
          cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
          cg2d_b(i,j,bi,bj) = 0.
#ifdef USE_NATURAL_BCS
     &     + freeSurfFac*_rA(i,j,bi,bj)*
     &       EmPmR(I,J,bi,bj)/deltaTMom
#endif
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=Nr,1,-1
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j) = _dyG(i,j,bi,bj)
     &      *drF(k)*_hFacW(i,j,k,bi,bj)
           vf(i,j) = _dxG(i,j,bi,bj)
     &      *drF(k)*_hFacS(i,j,k,bi,bj)
          ENDDO
         ENDDO
         CALL CALC_DIV_GHAT(
     I       bi,bj,1,sNx,1,sNy,K,
     I       uf,vf,
     U       cg2d_b,
     I       myThid)
        ENDDO
       ENDDO
      ENDDO

C--   Add source term arising from w=d/dt (p_s + p_nh)
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_NONHYDROSTATIC
        DO j=1,sNy
         DO i=1,sNx
          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
     &        *( etaN(i,j,bi,bj)
     &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
          cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
     &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
     &        *( etaN(i,j,bi,bj)
     &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
C-jmc
c    &      -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
c    &        *( cg2d_x(i,j,bi,bj) + cg3d_x(i,j,1,bi,bj) )
c    &        /deltaTMom/deltaTMom
C-jmc
         ENDDO
        ENDDO
#else
        DO j=1,sNy
         DO i=1,sNx
          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
     &        * etaN(i,j,bi,bj)
         ENDDO
        ENDDO
#endif

#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
         DO i=1,sNx
C Northern boundary
          IF (OB_Jn(I,bi,bj).NE.0) THEN
           cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
          ENDIF
C Southern boundary
          IF (OB_Js(I,bi,bj).NE.0) THEN
           cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
          ENDIF
         ENDDO
         DO j=1,sNy
C Eastern boundary
          IF (OB_Ie(J,bi,bj).NE.0) THEN
           cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(J,bi,bj).NE.0) THEN
           cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
          ENDIF
         ENDDO
        ENDIF
#endif
       ENDDO
      ENDDO


C--   Find the surface pressure using a two-dimensional conjugate
C--   gradient solver.
C     see CG2D_INTERNAL.h for the interface to this routine.
      tolerance=cg2dTargetResidual
      residual=0.
      numIters=cg2dMaxIters
      CALL CG2D(
     I           cg2d_b,
     U           cg2d_x,
     U           tolerance,
     U           residual,
     U           numIters,
     I           myThid )
      _EXCH_XY_R8(cg2d_x, myThid )

      _BEGIN_MASTER( myThid )
#ifndef ALLOW_AUTODIFF_TAMC
       WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
     &                         0, tolerance
       WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
     &                         numIters, residual
#endif
      _END_MASTER( )

C--   Transfert the 2D-solution to "etaN" :
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_NONHYDROSTATIC
      IF ( nonHydrostatic ) 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)
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j)=-_recip_dxC(i,j,bi,bj)*
     &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
           vf(i,j)=-_recip_dyC(i,j,bi,bj)*
     &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
          ENDDO
         ENDDO

#ifdef ALLOW_OBCS
         IF (useOBCS) THEN
          DO i=1,sNx+1
C Northern boundary
          IF (OB_Jn(I,bi,bj).NE.0) THEN
           vf(I,OB_Jn(I,bi,bj))=0.
          ENDIF
C Southern boundary
          IF (OB_Js(I,bi,bj).NE.0) THEN
           vf(I,OB_Js(I,bi,bj)+1)=0.
          ENDIF
          ENDDO
          DO j=1,sNy+1
C Eastern boundary
          IF (OB_Ie(J,bi,bj).NE.0) THEN
           uf(OB_Ie(J,bi,bj),J)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(J,bi,bj).NE.0) THEN
           uf(OB_Iw(J,bi,bj)+1,J)=0.
          ENDIF
          ENDDO
         ENDIF
#endif

         K=1
         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 )
     &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
     &          -wVel(i,j,k+1,bi,bj)
     &        )*_rA(i,j,bi,bj)/deltaTmom
          ENDDO
         ENDDO
         DO K=2,Nr-1
          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)
     &         -wVel(i,j,k+1,bi,bj)
     &        )*_rA(i,j,bi,bj)/deltaTmom

           ENDDO
          ENDDO
         ENDDO
         K=Nr
         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)
     &        )*_rA(i,j,bi,bj)/deltaTmom

          ENDDO
         ENDDO

#ifdef ALLOW_OBCS
         IF (useOBCS) THEN
          DO K=1,Nr
          DO i=1,sNx
C Northern boundary
          IF (OB_Jn(I,bi,bj).NE.0) THEN
           cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
          ENDIF
C Southern boundary
          IF (OB_Js(I,bi,bj).NE.0) THEN
           cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
          ENDIF
          ENDDO
          DO j=1,sNy
C Eastern boundary
          IF (OB_Ie(J,bi,bj).NE.0) THEN
           cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(J,bi,bj).NE.0) THEN
           cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
          ENDIF
          ENDDO
          ENDDO
         ENDIF
#endif

        ENDDO ! bi
       ENDDO ! bj

       CALL CG3D( myThid )
       _EXCH_XYZ_R8(cg3d_x, myThid )

      ENDIF
#endif

      RETURN
      END
1	heimbach	1.20	C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.19 2001/03/09 20:45:09 adcroft Exp $
2			C $Name: checkpoint37 $
3	cnh	1.1
4	adcroft	1.5	#include "CPP_OPTIONS.h"
5	cnh	1.1
6			CStartOfInterface
7			SUBROUTINE SOLVE_FOR_PRESSURE( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE SOLVE_FOR_PRESSURE \|
10			C \| o Controls inversion of two and/or three-dimensional \|
11			C \| elliptic problems for the pressure field. \|
12			C \==========================================================/
13	adcroft	1.8	IMPLICIT NONE
14	cnh	1.1
15	cnh	1.4	C == Global variables
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "DYNVARS.h"
20	adcroft	1.12	#include "GRID.h"
21	jmc	1.17	#include "SURFACE.h"
22	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
23			#include "CG3D.h"
24			#include "GW.h"
25	adcroft	1.12	#endif
26	adcroft	1.11	#ifdef ALLOW_OBCS
27	adcroft	1.9	#include "OBCS.h"
28	adcroft	1.11	#endif
29	cnh	1.4
30	cnh	1.1	C == Routine arguments ==
31			C myThid - Number of this instance of SOLVE_FOR_PRESSURE
32			INTEGER myThid
33			CEndOfInterface
34	cnh	1.4
35			C Local variables
36	jmc	1.17	C cg2d_x - Conjugate Gradient 2-D solver : Solution vector
37			C cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
38	cnh	1.6	INTEGER i,j,k,bi,bj
39	adcroft	1.9	_RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40			_RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41	jmc	1.17	_RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42			_RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43	adcroft	1.19	_RL tolerance,residual
44			INTEGER numIters
45	jmc	1.17
46			C-- Save previous solution & Initialise Vector solution and source term :
47			DO bj=myByLo(myThid),myByHi(myThid)
48			DO bi=myBxLo(myThid),myBxHi(myThid)
49			DO j=1-OLy,sNy+OLy
50			DO i=1-OLx,sNx+OLx
51			#ifdef INCLUDE_CD_CODE
52			etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
53			#endif
54	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
55	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
56			#ifdef USE_NATURAL_BCS
57	jmc	1.18	& + freeSurfFac_rA(i,j,bi,bj)
58	jmc	1.17	& EmPmR(I,J,bi,bj)/deltaTMom
59			#endif
60			ENDDO
61			ENDDO
62			ENDDO
63			ENDDO
64	adcroft	1.12
65			DO bj=myByLo(myThid),myByHi(myThid)
66			DO bi=myBxLo(myThid),myBxHi(myThid)
67			DO K=Nr,1,-1
68			DO j=1,sNy+1
69			DO i=1,sNx+1
70			uf(i,j) = _dyG(i,j,bi,bj)
71			& drF(k)_hFacW(i,j,k,bi,bj)
72			vf(i,j) = _dxG(i,j,bi,bj)
73			& drF(k)_hFacS(i,j,k,bi,bj)
74			ENDDO
75			ENDDO
76			CALL CALC_DIV_GHAT(
77			I bi,bj,1,sNx,1,sNy,K,
78			I uf,vf,
79	jmc	1.17	U cg2d_b,
80	adcroft	1.12	I myThid)
81			ENDDO
82			ENDDO
83			ENDDO
84	cnh	1.4
85	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
86			DO bj=myByLo(myThid),myByHi(myThid)
87			DO bi=myBxLo(myThid),myBxHi(myThid)
88	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
89	adcroft	1.12	DO j=1,sNy
90			DO i=1,sNx
91			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
92	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
93			& *( etaN(i,j,bi,bj)
94			& +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
95	adcroft	1.13	cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
96	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
97			& *( etaN(i,j,bi,bj)
98			& +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
99			C-jmc
100			c & -freeSurfFac_rA(i,j,bi,bj)recip_Bo(i,j,bi,bj)
101			c & *( cg2d_x(i,j,bi,bj) + cg3d_x(i,j,1,bi,bj) )
102			c & /deltaTMom/deltaTMom
103			C-jmc
104	adcroft	1.12	ENDDO
105			ENDDO
106	adcroft	1.13	#else
107	adcroft	1.12	DO j=1,sNy
108			DO i=1,sNx
109	adcroft	1.13	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
110	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
111			& * etaN(i,j,bi,bj)
112	adcroft	1.12	ENDDO
113			ENDDO
114			#endif
115
116			#ifdef ALLOW_OBCS
117	adcroft	1.14	IF (useOBCS) THEN
118	adcroft	1.12	DO i=1,sNx
119			C Northern boundary
120			IF (OB_Jn(I,bi,bj).NE.0) THEN
121			cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
122			ENDIF
123			C Southern boundary
124			IF (OB_Js(I,bi,bj).NE.0) THEN
125			cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
126			ENDIF
127			ENDDO
128			DO j=1,sNy
129			C Eastern boundary
130			IF (OB_Ie(J,bi,bj).NE.0) THEN
131			cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
132			ENDIF
133			C Western boundary
134			IF (OB_Iw(J,bi,bj).NE.0) THEN
135			cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
136			ENDIF
137			ENDDO
138			ENDIF
139			#endif
140			ENDDO
141			ENDDO
142
143
144	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
145			C-- gradient solver.
146	jmc	1.17	C see CG2D_INTERNAL.h for the interface to this routine.
147	adcroft	1.19	tolerance=cg2dTargetResidual
148			residual=0.
149			numIters=cg2dMaxIters
150	cnh	1.1	CALL CG2D(
151	cnh	1.6	I cg2d_b,
152			U cg2d_x,
153	adcroft	1.19	U tolerance,
154			U residual,
155			U numIters,
156	cnh	1.1	I myThid )
157	adcroft	1.19	_EXCH_XY_R8(cg2d_x, myThid )
158	cnh	1.1
159	adcroft	1.19	_BEGIN_MASTER( myThid )
160	heimbach	1.20	#ifndef ALLOW_AUTODIFF_TAMC
161	adcroft	1.19	WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
162			& 0, tolerance
163			WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
164			& numIters, residual
165	heimbach	1.20	#endif
166	adcroft	1.19	_END_MASTER( )
167	jmc	1.17
168			C-- Transfert the 2D-solution to "etaN" :
169			DO bj=myByLo(myThid),myByHi(myThid)
170			DO bi=myBxLo(myThid),myBxHi(myThid)
171			DO j=1-OLy,sNy+OLy
172			DO i=1-OLx,sNx+OLx
173	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
174	jmc	1.17	ENDDO
175			ENDDO
176			ENDDO
177			ENDDO
178	adcroft	1.10
179	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
180			IF ( nonHydrostatic ) THEN
181
182			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
183			C see CG3D.h for the interface to this routine.
184			DO bj=myByLo(myThid),myByHi(myThid)
185			DO bi=myBxLo(myThid),myBxHi(myThid)
186			DO j=1,sNy+1
187			DO i=1,sNx+1
188	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
189	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
190	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
191	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
192			ENDDO
193			ENDDO
194
195	adcroft	1.12	#ifdef ALLOW_OBCS
196	adcroft	1.14	IF (useOBCS) THEN
197	adcroft	1.9	DO i=1,sNx+1
198			C Northern boundary
199			IF (OB_Jn(I,bi,bj).NE.0) THEN
200			vf(I,OB_Jn(I,bi,bj))=0.
201			ENDIF
202			C Southern boundary
203			IF (OB_Js(I,bi,bj).NE.0) THEN
204			vf(I,OB_Js(I,bi,bj)+1)=0.
205			ENDIF
206			ENDDO
207			DO j=1,sNy+1
208			C Eastern boundary
209			IF (OB_Ie(J,bi,bj).NE.0) THEN
210			uf(OB_Ie(J,bi,bj),J)=0.
211			ENDIF
212			C Western boundary
213			IF (OB_Iw(J,bi,bj).NE.0) THEN
214			uf(OB_Iw(J,bi,bj)+1,J)=0.
215			ENDIF
216			ENDDO
217			ENDIF
218	adcroft	1.12	#endif
219	adcroft	1.9
220	adcroft	1.12	K=1
221			DO j=1,sNy
222			DO i=1,sNx
223			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
224			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
225			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
226			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
227			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
228	jmc	1.18	& +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
229			& -wVel(i,j,k+1,bi,bj)
230	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
231			ENDDO
232			ENDDO
233			DO K=2,Nr-1
234	adcroft	1.9	DO j=1,sNy
235			DO i=1,sNx
236			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
237			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
238			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
239			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
240			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
241	adcroft	1.12	& +( wVel(i,j,k ,bi,bj)
242			& -wVel(i,j,k+1,bi,bj)
243			& )*_rA(i,j,bi,bj)/deltaTmom
244
245	adcroft	1.9	ENDDO
246			ENDDO
247			ENDDO
248	adcroft	1.12	K=Nr
249			DO j=1,sNy
250			DO i=1,sNx
251			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
252			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
253			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
254			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
255			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
256			& +( wVel(i,j,k ,bi,bj)
257			& )*_rA(i,j,bi,bj)/deltaTmom
258
259			ENDDO
260			ENDDO
261
262			#ifdef ALLOW_OBCS
263	adcroft	1.14	IF (useOBCS) THEN
264	adcroft	1.12	DO K=1,Nr
265			DO i=1,sNx
266			C Northern boundary
267			IF (OB_Jn(I,bi,bj).NE.0) THEN
268			cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
269			ENDIF
270			C Southern boundary
271			IF (OB_Js(I,bi,bj).NE.0) THEN
272			cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
273			ENDIF
274			ENDDO
275			DO j=1,sNy
276			C Eastern boundary
277			IF (OB_Ie(J,bi,bj).NE.0) THEN
278			cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
279			ENDIF
280			C Western boundary
281			IF (OB_Iw(J,bi,bj).NE.0) THEN
282			cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
283			ENDIF
284			ENDDO
285			ENDDO
286			ENDIF
287			#endif
288	adcroft	1.9
289			ENDDO ! bi
290			ENDDO ! bj
291
292			CALL CG3D( myThid )
293	adcroft	1.10	_EXCH_XYZ_R8(cg3d_x, myThid )
294	adcroft	1.9
295			ENDIF
296			#endif
297	cnh	1.1
298			RETURN
299			END