model/src/solve_for_pressure.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.24 2001/06/06 15:14:06 adcroft Exp $
C $Name:  $

#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 "SOLVE_FOR_PRESSURE3D.h"
#include "GW.h"
#endif
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif
#include "SOLVE_FOR_PRESSURE.h"

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

C     == Local variables ==
      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 firstResidual,lastResidual
      INTEGER numIters
      CHARACTER*(MAX_LEN_MBUF) msgBuf

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
          etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
          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)
     &          +phi_nh(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)
     &          +phi_nh(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) + phi_nh(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

#ifndef EXCLUDE_DEBUGMODE
      IF (debugMode) THEN
       CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif

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

#ifndef EXCLUDE_DEBUGMODE
      IF (debugMode) THEN
       CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif

      _BEGIN_MASTER( myThid )
      WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      _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

      firstResidual=0.
      lastResidual=0.
      numIters=cg2dMaxIters
      CALL CG3D(
     U           cg3d_b,
     U           phi_nh,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      _EXCH_XYZ_R8(phi_nh, myThid )

      _BEGIN_MASTER( myThid )
      WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      _END_MASTER( )

      ENDIF
#endif

      RETURN
      END
1	adcroft	1.25	C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.24 2001/06/06 15:14:06 adcroft Exp $
2	heimbach	1.21	C $Name: $
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	adcroft	1.25	#include "SOLVE_FOR_PRESSURE3D.h"
24	adcroft	1.9	#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	adcroft	1.22	#include "SOLVE_FOR_PRESSURE.h"
30	cnh	1.4
31	cnh	1.1	C == Routine arguments ==
32			C myThid - Number of this instance of SOLVE_FOR_PRESSURE
33			INTEGER myThid
34			CEndOfInterface
35	cnh	1.4
36	adcroft	1.22	C == Local variables ==
37	cnh	1.6	INTEGER i,j,k,bi,bj
38	adcroft	1.9	_RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39			_RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40	adcroft	1.22	_RL firstResidual,lastResidual
41	adcroft	1.19	INTEGER numIters
42	adcroft	1.25	CHARACTER*(MAX_LEN_MBUF) msgBuf
43	jmc	1.17
44			C-- Save previous solution & Initialise Vector solution and source term :
45			DO bj=myByLo(myThid),myByHi(myThid)
46			DO bi=myBxLo(myThid),myBxHi(myThid)
47			DO j=1-OLy,sNy+OLy
48			DO i=1-OLx,sNx+OLx
49			etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
50	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
51	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
52			#ifdef USE_NATURAL_BCS
53	jmc	1.18	& + freeSurfFac_rA(i,j,bi,bj)
54	jmc	1.17	& EmPmR(I,J,bi,bj)/deltaTMom
55			#endif
56			ENDDO
57			ENDDO
58			ENDDO
59			ENDDO
60	adcroft	1.12
61			DO bj=myByLo(myThid),myByHi(myThid)
62			DO bi=myBxLo(myThid),myBxHi(myThid)
63			DO K=Nr,1,-1
64			DO j=1,sNy+1
65			DO i=1,sNx+1
66			uf(i,j) = _dyG(i,j,bi,bj)
67			& drF(k)_hFacW(i,j,k,bi,bj)
68			vf(i,j) = _dxG(i,j,bi,bj)
69			& drF(k)_hFacS(i,j,k,bi,bj)
70			ENDDO
71			ENDDO
72			CALL CALC_DIV_GHAT(
73			I bi,bj,1,sNx,1,sNy,K,
74			I uf,vf,
75	jmc	1.17	U cg2d_b,
76	adcroft	1.12	I myThid)
77			ENDDO
78			ENDDO
79			ENDDO
80	cnh	1.4
81	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
82			DO bj=myByLo(myThid),myByHi(myThid)
83			DO bi=myBxLo(myThid),myBxHi(myThid)
84	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
85	adcroft	1.12	DO j=1,sNy
86			DO i=1,sNx
87			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
88	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
89			& *( etaN(i,j,bi,bj)
90	adcroft	1.25	& +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
91	adcroft	1.13	cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
92	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
93			& *( etaN(i,j,bi,bj)
94	adcroft	1.25	& +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
95	jmc	1.18	C-jmc
96			c & -freeSurfFac_rA(i,j,bi,bj)recip_Bo(i,j,bi,bj)
97	adcroft	1.25	c & *( cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj) )
98	jmc	1.18	c & /deltaTMom/deltaTMom
99			C-jmc
100	adcroft	1.12	ENDDO
101			ENDDO
102	adcroft	1.13	#else
103	adcroft	1.12	DO j=1,sNy
104			DO i=1,sNx
105	adcroft	1.13	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
106	jmc	1.18	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
107			& * etaN(i,j,bi,bj)
108	adcroft	1.12	ENDDO
109			ENDDO
110			#endif
111
112			#ifdef ALLOW_OBCS
113	adcroft	1.14	IF (useOBCS) THEN
114	adcroft	1.12	DO i=1,sNx
115			C Northern boundary
116			IF (OB_Jn(I,bi,bj).NE.0) THEN
117			cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
118			ENDIF
119			C Southern boundary
120			IF (OB_Js(I,bi,bj).NE.0) THEN
121			cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
122			ENDIF
123			ENDDO
124			DO j=1,sNy
125			C Eastern boundary
126			IF (OB_Ie(J,bi,bj).NE.0) THEN
127			cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
128			ENDIF
129			C Western boundary
130			IF (OB_Iw(J,bi,bj).NE.0) THEN
131			cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
132			ENDIF
133			ENDDO
134			ENDIF
135			#endif
136			ENDDO
137			ENDDO
138
139	adcroft	1.23	#ifndef EXCLUDE_DEBUGMODE
140	adcroft	1.24	IF (debugMode) THEN
141	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
142			& myThid)
143	adcroft	1.24	ENDIF
144	adcroft	1.23	#endif
145	adcroft	1.12
146	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
147			C-- gradient solver.
148	adcroft	1.22	C see CG2D.h for the interface to this routine.
149			firstResidual=0.
150			lastResidual=0.
151	adcroft	1.19	numIters=cg2dMaxIters
152	cnh	1.1	CALL CG2D(
153	adcroft	1.22	U cg2d_b,
154	cnh	1.6	U cg2d_x,
155	adcroft	1.22	O firstResidual,
156			O lastResidual,
157	adcroft	1.19	U numIters,
158	cnh	1.1	I myThid )
159	adcroft	1.19	_EXCH_XY_R8(cg2d_x, myThid )
160	adcroft	1.23
161			#ifndef EXCLUDE_DEBUGMODE
162	adcroft	1.24	IF (debugMode) THEN
163	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
164			& myThid)
165	adcroft	1.24	ENDIF
166	adcroft	1.23	#endif
167	cnh	1.1
168	adcroft	1.19	_BEGIN_MASTER( myThid )
169	adcroft	1.25	WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
170			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
171			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
172			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
173			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
174			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
175	adcroft	1.19	_END_MASTER( )
176	jmc	1.17
177			C-- Transfert the 2D-solution to "etaN" :
178			DO bj=myByLo(myThid),myByHi(myThid)
179			DO bi=myBxLo(myThid),myBxHi(myThid)
180			DO j=1-OLy,sNy+OLy
181			DO i=1-OLx,sNx+OLx
182	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
183	jmc	1.17	ENDDO
184			ENDDO
185			ENDDO
186			ENDDO
187	adcroft	1.10
188	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
189			IF ( nonHydrostatic ) THEN
190
191			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
192			C see CG3D.h for the interface to this routine.
193			DO bj=myByLo(myThid),myByHi(myThid)
194			DO bi=myBxLo(myThid),myBxHi(myThid)
195			DO j=1,sNy+1
196			DO i=1,sNx+1
197	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
198	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
199	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
200	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
201			ENDDO
202			ENDDO
203
204	adcroft	1.12	#ifdef ALLOW_OBCS
205	adcroft	1.14	IF (useOBCS) THEN
206	adcroft	1.9	DO i=1,sNx+1
207			C Northern boundary
208			IF (OB_Jn(I,bi,bj).NE.0) THEN
209			vf(I,OB_Jn(I,bi,bj))=0.
210			ENDIF
211			C Southern boundary
212			IF (OB_Js(I,bi,bj).NE.0) THEN
213			vf(I,OB_Js(I,bi,bj)+1)=0.
214			ENDIF
215			ENDDO
216			DO j=1,sNy+1
217			C Eastern boundary
218			IF (OB_Ie(J,bi,bj).NE.0) THEN
219			uf(OB_Ie(J,bi,bj),J)=0.
220			ENDIF
221			C Western boundary
222			IF (OB_Iw(J,bi,bj).NE.0) THEN
223			uf(OB_Iw(J,bi,bj)+1,J)=0.
224			ENDIF
225			ENDDO
226			ENDIF
227	adcroft	1.12	#endif
228	adcroft	1.9
229	adcroft	1.12	K=1
230			DO j=1,sNy
231			DO i=1,sNx
232			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
233			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
234			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
235			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
236			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
237	jmc	1.18	& +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
238			& -wVel(i,j,k+1,bi,bj)
239	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
240			ENDDO
241			ENDDO
242			DO K=2,Nr-1
243	adcroft	1.9	DO j=1,sNy
244			DO i=1,sNx
245			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
246			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
247			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
248			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
249			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
250	adcroft	1.12	& +( wVel(i,j,k ,bi,bj)
251			& -wVel(i,j,k+1,bi,bj)
252			& )*_rA(i,j,bi,bj)/deltaTmom
253
254	adcroft	1.9	ENDDO
255			ENDDO
256			ENDDO
257	adcroft	1.12	K=Nr
258			DO j=1,sNy
259			DO i=1,sNx
260			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
261			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
262			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
263			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
264			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
265			& +( wVel(i,j,k ,bi,bj)
266			& )*_rA(i,j,bi,bj)/deltaTmom
267
268			ENDDO
269			ENDDO
270
271			#ifdef ALLOW_OBCS
272	adcroft	1.14	IF (useOBCS) THEN
273	adcroft	1.12	DO K=1,Nr
274			DO i=1,sNx
275			C Northern boundary
276			IF (OB_Jn(I,bi,bj).NE.0) THEN
277			cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
278			ENDIF
279			C Southern boundary
280			IF (OB_Js(I,bi,bj).NE.0) THEN
281			cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
282			ENDIF
283			ENDDO
284			DO j=1,sNy
285			C Eastern boundary
286			IF (OB_Ie(J,bi,bj).NE.0) THEN
287			cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
288			ENDIF
289			C Western boundary
290			IF (OB_Iw(J,bi,bj).NE.0) THEN
291			cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
292			ENDIF
293			ENDDO
294			ENDDO
295			ENDIF
296			#endif
297	adcroft	1.9
298			ENDDO ! bi
299			ENDDO ! bj
300
301	adcroft	1.25	firstResidual=0.
302			lastResidual=0.
303			numIters=cg2dMaxIters
304			CALL CG3D(
305			U cg3d_b,
306			U phi_nh,
307			O firstResidual,
308			O lastResidual,
309			U numIters,
310			I myThid )
311			_EXCH_XYZ_R8(phi_nh, myThid )
312
313			_BEGIN_MASTER( myThid )
314			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
315			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
316			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
317			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
318			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
319			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320			_END_MASTER( )
321	adcroft	1.9
322			ENDIF
323			#endif
324	cnh	1.1
325			RETURN
326			END