model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.27 2001/09/26 18:09:16 cnh Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: SOLVE_FOR_PRESSURE
C     !INTERFACE:
      SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid) 

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SOLVE_FOR_PRESSURE                             
C     | o Controls inversion of two and/or three-dimensional      
C     |   elliptic problems for the pressure field.               
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"
#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     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime - Current time in simulation
C     myIter - Current iteration number in simulation
C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
      _RL myTime
      INTEGER myIter
      INTEGER myThid 

C     !LOCAL VARIABLES:
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
CEOP

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