model/src/ini_cg3d.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg3d.F,v 1.9 2001/02/04 14:38:47 cnh Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE INI_CG3D( myThid )
C     /==========================================================\
C     | SUBROUTINE INI_CG3D                                      |
C     | o Initialise 3d conjugate gradient solver operators.     |
C     |==========================================================|
C     | These arrays are purely a function of the basin geom.    |
C     | We set then here once and them use then repeatedly.      |
C     \==========================================================/
      IMPLICIT NONE

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

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface

#ifdef ALLOW_NONHYDROSTATIC

C     === Local variables ===
C     xG, yG - Global coordinate location.
C     zG
C     iG, jG - Global coordinate index
C     bi,bj  - Loop counters
C     faceArea - Temporary used to hold cell face areas.
C     I,J,K
C     myNorm - Work variable used in clculating normalisation factor
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER bi, bj
      INTEGER I,  J, K
      _RL     faceArea
      _RS     myNorm
      _RL     aCw, aCs
      _RL     theRecip_Dr
      _RL     aU, aL, aW, aE, aN, aS, aC

CcnhDebugStarts
      _RL    phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
CcnhDebugEnds

C--   Initialise laplace operator
C     aW3d: Ax/dX
C     aS3d: Ay/dY
C     aV3d: Ar/dR
      myNorm = 0. _d 0
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=1,Nr
         DO J=1,sNy
          DO I=1,sNx
           faceArea = _dyG(I,J,bi,bj)*drF(K)
     &                *_hFacW(I,J,K,bi,bj)
           aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
           faceArea = _dxG(I,J,bi,bj)*drF(K)
     &                *_hFacS(I,J,K,bi,bj)
           aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
           myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
           myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
          ENDDO
         ENDDO
        ENDDO
        DO K=1,1
         DO J=1,sNy
          DO I=1,sNx
           aV3d(I,J,K,bi,bj) =  0.
           myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
          ENDDO
         ENDDO
        ENDDO
        DO K=2,Nr
         DO J=1,sNy
          DO I=1,sNx
           IF ( _hFacC(i,j,k,bi,bj) .EQ. 0. ) THEN
            faceArea = 0.
           ELSE
            faceArea = _rA(I,J,bi,bj)
           ENDIF
           theRecip_Dr = 
     &      drC(K)
caja &      drF(K  )*_hFacC(i,j,k  ,bi,bj)*0.5
caja &     +drF(K-1)*_hFacC(i,j,k-1,bi,bj)*0.5
           IF ( theRecip_Dr .NE. 0. ) 
     &      theRecip_Dr = 1. _d 0/theRecip_Dr
           aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr
     &                 *horiVertRatio*horiVertRatio 
           myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_OBCS
        DO K=1,Nr
         IF (useOBCS) THEN
          DO I=1,sNx
           IF (OB_Jn(I,bi,bj).NE.0) THEN
            aS3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
            aS3d(I,OB_Jn(I,bi,bj)+1,K,bi,bj)=0.
            aW3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
            aW3d(I+1,OB_Jn(I,bi,bj),K,bi,bj)=0.
            aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
           ENDIF
           IF (OB_Js(I,bi,bj).NE.0) THEN
            aS3d(I,OB_Js(I,bi,bj)+1,K,bi,bj)=0.
            aS3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
            aW3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
            aW3d(I+1,OB_Js(I,bi,bj),K,bi,bj)=0.
            aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
           ENDIF
          ENDDO
          DO J=1,sNy
           IF (OB_Ie(J,bi,bj).NE.0) THEN
            aW3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
            aW3d(OB_Ie(J,bi,bj)+1,J,K,bi,bj)=0.
            aS3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
            aS3d(OB_Ie(J,bi,bj),J+1,K,bi,bj)=0.
            aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
           ENDIF
           IF (OB_Iw(J,bi,bj).NE.0) THEN
            aW3d(OB_Iw(J,bi,bj)+1,J,K,bi,bj)=0.
            aW3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
            aS3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
            aS3d(OB_Iw(J,bi,bj),J+1,K,bi,bj)=0.
            aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
           ENDIF
          ENDDO
         ENDIF
        ENDDO
#endif
       ENDDO
      ENDDO
      _GLOBAL_MAX_R4( myNorm, myThid )
      IF ( myNorm .NE. 0. _d 0 ) THEN
       myNorm = 1. _d 0/myNorm
      ELSE
       myNorm = 1. _d 0
      ENDIF
       cg3dNorm = myNorm
      _BEGIN_MASTER( myThid )
CcnhDebugStarts
       WRITE(msgBuf,'(A,E40.25)') '// CG3D normalisation factor = '
     &                , cg3dNorm
       CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
       WRITE(msgBuf,*) '                               '
       CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
CcnhDebugEnds
      _END_MASTER( myThid )
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=1,Nr
         DO J=1,sNy
          DO I=1,sNx
           aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
           aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
           aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      
C--   Update overlap regions
      _EXCH_XYZ_R4(aW3d, myThid)
      _EXCH_XYZ_R4(aS3d, myThid)
      _EXCH_XYZ_R4(aV3d, myThid)
CcnhDebugStarts
C     CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
C     CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
CcnhDebugEnds

C--   Initialise preconditioner
C     For now PC is just the identity. Change to 
C     be LU factorization of d2/dz2 later. Note
C     check for consistency with S/R CG3D before
C     assuming zML is lower and zMU is upper!
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=1,Nr  
         DO J=1,sNy
          DO I=1,sNx
           aW = aW3d(I  ,J,K,bi,bj)
           aE = aW3d(I+1,J,K,bi,bj)
           aN = aS3d(I,J+1,K,bi,bj)
           aS = aS3d(I,J  ,K,bi,bj)
           IF ( K .NE. 1 ) THEN
            aU = aV3d(I,J,K,bi,bj)
           ELSE
            aU = 0
           ENDIF
           IF ( K .NE. Nr ) THEN
            aL = aV3d(I,J,K+1,bi,bj)
           ELSE
            aL = 0
           ENDIF
           aC = -aW-aE-aN-aS-aU-aL
           IF ( K .EQ. 1 .AND. aC .NE. 0. ) THEN
            aC = aC
     &         -freeSurfFac*myNorm*(horiVertRatio/gravity)*
     &          rA(I,J,bi,bj)/deltaTMom/deltaTMom
           ENDIF
           IF ( aC .NE. 0. ) THEN
            zMC(i,j,k,bi,bj) = aC
            zMU(i,j,k,bi,bj) = aL
            zML(i,j,k,bi,bj) = aU
CcnhDebugStarts
C           zMC(i,j,k,bi,bj) = 1.
C           zMU(i,j,k,bi,bj) = 0.
C           zML(i,j,k,bi,bj) = 0.
CcnhDebugEnds
           ELSE
            zMC(i,j,k,bi,bj) = 1. _d 0
            zMU(i,j,k,bi,bj) = 0.
            zML(i,j,k,bi,bj) = 0.
           ENDIF
          ENDDO
         ENDDO
        ENDDO
        DO J=1,sNy
         DO I=1,sNx
          zMC(i,j,1,bi,bj)=
     &     1. _d 0 / zMC(i,j,1,bi,bj)
          zMU(i,j,1,bi,bj)=
     &     zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
         ENDDO
        ENDDO
        DO K=2,Nr  
         DO J=1,sNy
          DO I=1,sNx
           zMC(i,j,k,bi,bj) = 1. _d 0 /
     &     (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
           zMU(i,j,k,bi,bj)=zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        DO K=1,Nr  
         DO J=1,sNy
          DO I=1,sNx
           aW = aW3d(I  ,J,K,bi,bj)
           aE = aW3d(I+1,J,K,bi,bj)
           aN = aS3d(I,J+1,K,bi,bj)
           aS = aS3d(I,J  ,K,bi,bj)
           IF ( K .NE. 1 ) THEN
            aU = aV3d(I,J,K-1,bi,bj)
           ELSE
            aU = 0
           ENDIF
           IF ( K .NE. Nr ) THEN
            aL = aV3d(I,J,K+1,bi,bj)
           ELSE
            aL = 0
           ENDIF
           aC = -aW-aE-aN-aS-aU-aL
           IF ( aC .EQ. 0. ) THEN
            zMC(i,j,k,bi,bj) = 1.
            zML(i,j,k,bi,bj) = 0.
            zMU(i,j,k,bi,bj) = 0.
CcnhDebugStarts
C          ELSE
C           zMC(i,j,k,bi,bj) = 1.
C           zML(i,j,k,bi,bj) = 0.
C           zMU(i,j,k,bi,bj) = 0.
CcnhDEbugEnds
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C--   Update overlap regions
      _EXCH_XYZ_R4(zMC, myThid)
      _EXCH_XYZ_R4(zML, myThid)
      _EXCH_XYZ_R4(zMU, myThid)

CcnhDebugStarts
      DO k=1,nr
      DO j=1-OLy,sNy+OLy
      DO i=1-OLx,sNx+OLx
       phi(i,j,1,1) = zMc(i,j,k,1,1)
      ENDDO
      ENDDO
C     CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
      ENDDO
C     CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
C     CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
CcnhDebugEnds

C--   Set default values for initial guess and RHS
      IF ( startTime .EQ. 0 ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO K=1,Nr
          DO J=1-OLy,sNy+OLy
           DO I=1-OLx,sNx+OLx
            cg3d_x(I,J,K,bi,bj) = 0. _d 0
            cg3d_b(I,J,K,bi,bj) = 0. _d 0
#ifdef INCLUDE_CD_CODE
            cg3d_xNM1(I,J,K,bi,bj) = 0. _d 0
#endif
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
C--    Update overlap regions
       _EXCH_XYZ_R8(cg3d_x, myThid)
       _EXCH_XYZ_R8(cg3d_b, myThid)
#ifdef INCLUDE_CD_CODE
       _EXCH_XYZ_R8(cg3d_xNM1, myThid)
#endif
      ENDIF

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	jmc	1.10	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg3d.F,v 1.9 2001/02/04 14:38:47 cnh Exp $
2			C $Name: $
3	adcroft	1.1
4			#include "CPP_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE INI_CG3D( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE INI_CG3D \|
10			C \| o Initialise 3d conjugate gradient solver operators. \|
11			C \|==========================================================\|
12			C \| These arrays are purely a function of the basin geom. \|
13			C \| We set then here once and them use then repeatedly. \|
14			C \==========================================================/
15	adcroft	1.3	IMPLICIT NONE
16	adcroft	1.1
17			C === Global variables ===
18			#include "SIZE.h"
19			#include "EEPARAMS.h"
20			#include "PARAMS.h"
21			#include "GRID.h"
22	jmc	1.10	c_jmc #include "DYNVARS.h"
23	adcroft	1.1	#include "CG3D.h"
24	adcroft	1.4	#ifdef ALLOW_OBCS
25	adcroft	1.1	#include "OBCS.h"
26	adcroft	1.4	#endif
27	adcroft	1.1
28			C === Routine arguments ===
29			C myThid - Thread no. that called this routine.
30			INTEGER myThid
31			CEndOfInterface
32
33	adcroft	1.6	#ifdef ALLOW_NONHYDROSTATIC
34
35	adcroft	1.1	C === Local variables ===
36			C xG, yG - Global coordinate location.
37			C zG
38			C iG, jG - Global coordinate index
39			C bi,bj - Loop counters
40			C faceArea - Temporary used to hold cell face areas.
41			C I,J,K
42			C myNorm - Work variable used in clculating normalisation factor
43			CHARACTER*(MAX_LEN_MBUF) msgBuf
44			INTEGER bi, bj
45			INTEGER I, J, K
46			_RL faceArea
47			_RS myNorm
48			_RL aCw, aCs
49			_RL theRecip_Dr
50			_RL aU, aL, aW, aE, aN, aS, aC
51
52			CcnhDebugStarts
53			_RL phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54			CcnhDebugEnds
55
56			C-- Initialise laplace operator
57			C aW3d: Ax/dX
58			C aS3d: Ay/dY
59			C aV3d: Ar/dR
60			myNorm = 0. _d 0
61			DO bj=myByLo(myThid),myByHi(myThid)
62			DO bi=myBxLo(myThid),myBxHi(myThid)
63			DO K=1,Nr
64			DO J=1,sNy
65			DO I=1,sNx
66			faceArea = _dyG(I,J,bi,bj)*drF(K)
67			& *_hFacW(I,J,K,bi,bj)
68	jmc	1.10	aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
69	adcroft	1.1	faceArea = _dxG(I,J,bi,bj)*drF(K)
70			& *_hFacS(I,J,K,bi,bj)
71	jmc	1.10	aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
72	adcroft	1.1	myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
73			myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
74			ENDDO
75			ENDDO
76			ENDDO
77			DO K=1,1
78			DO J=1,sNy
79			DO I=1,sNx
80			aV3d(I,J,K,bi,bj) = 0.
81			myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
82			ENDDO
83			ENDDO
84			ENDDO
85			DO K=2,Nr
86			DO J=1,sNy
87			DO I=1,sNx
88			IF ( _hFacC(i,j,k,bi,bj) .EQ. 0. ) THEN
89			faceArea = 0.
90			ELSE
91			faceArea = _rA(I,J,bi,bj)
92			ENDIF
93			theRecip_Dr =
94			& drC(K)
95			caja & drF(K )_hFacC(i,j,k ,bi,bj)0.5
96			caja & +drF(K-1)_hFacC(i,j,k-1,bi,bj)0.5
97			IF ( theRecip_Dr .NE. 0. )
98			& theRecip_Dr = 1. _d 0/theRecip_Dr
99	jmc	1.10	aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr
100			& horiVertRatiohoriVertRatio
101	adcroft	1.1	myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
102			ENDDO
103			ENDDO
104	adcroft	1.5	ENDDO
105	adcroft	1.4	#ifdef ALLOW_OBCS
106	adcroft	1.5	DO K=1,Nr
107	adcroft	1.8	IF (useOBCS) THEN
108	adcroft	1.1	DO I=1,sNx
109			IF (OB_Jn(I,bi,bj).NE.0) THEN
110			aS3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
111			aS3d(I,OB_Jn(I,bi,bj)+1,K,bi,bj)=0.
112			aW3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
113			aW3d(I+1,OB_Jn(I,bi,bj),K,bi,bj)=0.
114	adcroft	1.7	aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
115	adcroft	1.1	ENDIF
116			IF (OB_Js(I,bi,bj).NE.0) THEN
117			aS3d(I,OB_Js(I,bi,bj)+1,K,bi,bj)=0.
118			aS3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
119			aW3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
120			aW3d(I+1,OB_Js(I,bi,bj),K,bi,bj)=0.
121	adcroft	1.7	aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
122	adcroft	1.1	ENDIF
123			ENDDO
124			DO J=1,sNy
125			IF (OB_Ie(J,bi,bj).NE.0) THEN
126			aW3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
127			aW3d(OB_Ie(J,bi,bj)+1,J,K,bi,bj)=0.
128			aS3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
129			aS3d(OB_Ie(J,bi,bj),J+1,K,bi,bj)=0.
130	adcroft	1.8	aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
131	adcroft	1.1	ENDIF
132			IF (OB_Iw(J,bi,bj).NE.0) THEN
133			aW3d(OB_Iw(J,bi,bj)+1,J,K,bi,bj)=0.
134			aW3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
135			aS3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
136			aS3d(OB_Iw(J,bi,bj),J+1,K,bi,bj)=0.
137	adcroft	1.8	aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
138	adcroft	1.1	ENDIF
139			ENDDO
140			ENDIF
141	adcroft	1.5	ENDDO
142	adcroft	1.4	#endif
143	adcroft	1.1	ENDDO
144			ENDDO
145	adcroft	1.2	_GLOBAL_MAX_R4( myNorm, myThid )
146			IF ( myNorm .NE. 0. _d 0 ) THEN
147			myNorm = 1. _d 0/myNorm
148	adcroft	1.1	ELSE
149			myNorm = 1. _d 0
150			ENDIF
151			cg3dNorm = myNorm
152			_BEGIN_MASTER( myThid )
153			CcnhDebugStarts
154			WRITE(msgBuf,'(A,E40.25)') '// CG3D normalisation factor = '
155			& , cg3dNorm
156			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
157			WRITE(msgBuf,*) ' '
158			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
159			CcnhDebugEnds
160			_END_MASTER( myThid )
161			DO bj=myByLo(myThid),myByHi(myThid)
162			DO bi=myBxLo(myThid),myBxHi(myThid)
163			DO K=1,Nr
164			DO J=1,sNy
165			DO I=1,sNx
166			aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
167			aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
168			aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
169			ENDDO
170			ENDDO
171			ENDDO
172			ENDDO
173			ENDDO
174
175			C-- Update overlap regions
176			_EXCH_XYZ_R4(aW3d, myThid)
177			_EXCH_XYZ_R4(aS3d, myThid)
178			_EXCH_XYZ_R4(aV3d, myThid)
179			CcnhDebugStarts
180			C CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
181			C CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
182			CcnhDebugEnds
183
184			C-- Initialise preconditioner
185			C For now PC is just the identity. Change to
186			C be LU factorization of d2/dz2 later. Note
187			C check for consistency with S/R CG3D before
188			C assuming zML is lower and zMU is upper!
189			DO bj=myByLo(myThid),myByHi(myThid)
190			DO bi=myBxLo(myThid),myBxHi(myThid)
191			DO K=1,Nr
192			DO J=1,sNy
193			DO I=1,sNx
194			aW = aW3d(I ,J,K,bi,bj)
195			aE = aW3d(I+1,J,K,bi,bj)
196			aN = aS3d(I,J+1,K,bi,bj)
197			aS = aS3d(I,J ,K,bi,bj)
198			IF ( K .NE. 1 ) THEN
199			aU = aV3d(I,J,K,bi,bj)
200			ELSE
201			aU = 0
202			ENDIF
203			IF ( K .NE. Nr ) THEN
204			aL = aV3d(I,J,K+1,bi,bj)
205			ELSE
206			aL = 0
207			ENDIF
208			aC = -aW-aE-aN-aS-aU-aL
209			IF ( K .EQ. 1 .AND. aC .NE. 0. ) THEN
210			aC = aC
211	jmc	1.10	& -freeSurfFacmyNorm(horiVertRatio/gravity)*
212			& rA(I,J,bi,bj)/deltaTMom/deltaTMom
213	adcroft	1.1	ENDIF
214			IF ( aC .NE. 0. ) THEN
215			zMC(i,j,k,bi,bj) = aC
216			zMU(i,j,k,bi,bj) = aL
217			zML(i,j,k,bi,bj) = aU
218			CcnhDebugStarts
219			C zMC(i,j,k,bi,bj) = 1.
220			C zMU(i,j,k,bi,bj) = 0.
221			C zML(i,j,k,bi,bj) = 0.
222			CcnhDebugEnds
223			ELSE
224			zMC(i,j,k,bi,bj) = 1. _d 0
225			zMU(i,j,k,bi,bj) = 0.
226			zML(i,j,k,bi,bj) = 0.
227			ENDIF
228			ENDDO
229			ENDDO
230			ENDDO
231			DO J=1,sNy
232			DO I=1,sNx
233			zMC(i,j,1,bi,bj)=
234			& 1. _d 0 / zMC(i,j,1,bi,bj)
235			zMU(i,j,1,bi,bj)=
236			& zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
237			ENDDO
238			ENDDO
239			DO K=2,Nr
240			DO J=1,sNy
241			DO I=1,sNx
242			zMC(i,j,k,bi,bj) = 1. _d 0 /
243			& (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
244			zMU(i,j,k,bi,bj)=zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
245			ENDDO
246			ENDDO
247			ENDDO
248			DO K=1,Nr
249			DO J=1,sNy
250			DO I=1,sNx
251			aW = aW3d(I ,J,K,bi,bj)
252			aE = aW3d(I+1,J,K,bi,bj)
253			aN = aS3d(I,J+1,K,bi,bj)
254			aS = aS3d(I,J ,K,bi,bj)
255			IF ( K .NE. 1 ) THEN
256			aU = aV3d(I,J,K-1,bi,bj)
257			ELSE
258			aU = 0
259			ENDIF
260			IF ( K .NE. Nr ) THEN
261			aL = aV3d(I,J,K+1,bi,bj)
262			ELSE
263			aL = 0
264			ENDIF
265			aC = -aW-aE-aN-aS-aU-aL
266			IF ( aC .EQ. 0. ) THEN
267	adcroft	1.5	zMC(i,j,k,bi,bj) = 1.
268	adcroft	1.1	zML(i,j,k,bi,bj) = 0.
269			zMU(i,j,k,bi,bj) = 0.
270			CcnhDebugStarts
271			C ELSE
272			C zMC(i,j,k,bi,bj) = 1.
273			C zML(i,j,k,bi,bj) = 0.
274			C zMU(i,j,k,bi,bj) = 0.
275			CcnhDEbugEnds
276			ENDIF
277			ENDDO
278			ENDDO
279			ENDDO
280			ENDDO
281			ENDDO
282			C-- Update overlap regions
283			_EXCH_XYZ_R4(zMC, myThid)
284			_EXCH_XYZ_R4(zML, myThid)
285			_EXCH_XYZ_R4(zMU, myThid)
286
287			CcnhDebugStarts
288			DO k=1,nr
289			DO j=1-OLy,sNy+OLy
290			DO i=1-OLx,sNx+OLx
291			phi(i,j,1,1) = zMc(i,j,k,1,1)
292			ENDDO
293			ENDDO
294			C CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
295			ENDDO
296			C CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
297			C CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
298			CcnhDebugEnds
299
300			C-- Set default values for initial guess and RHS
301			IF ( startTime .EQ. 0 ) THEN
302			DO bj=myByLo(myThid),myByHi(myThid)
303			DO bi=myBxLo(myThid),myBxHi(myThid)
304			DO K=1,Nr
305			DO J=1-OLy,sNy+OLy
306			DO I=1-OLx,sNx+OLx
307			cg3d_x(I,J,K,bi,bj) = 0. _d 0
308			cg3d_b(I,J,K,bi,bj) = 0. _d 0
309	adcroft	1.2	#ifdef INCLUDE_CD_CODE
310	adcroft	1.1	cg3d_xNM1(I,J,K,bi,bj) = 0. _d 0
311			#endif
312			ENDDO
313			ENDDO
314			ENDDO
315			ENDDO
316			ENDDO
317			C-- Update overlap regions
318			_EXCH_XYZ_R8(cg3d_x, myThid)
319			_EXCH_XYZ_R8(cg3d_b, myThid)
320	adcroft	1.2	#ifdef INCLUDE_CD_CODE
321	adcroft	1.1	_EXCH_XYZ_R8(cg3d_xNM1, myThid)
322			#endif
323			ENDIF
324
325			#endif /* ALLOW_NONHYDROSTATIC */
326
327			RETURN
328			END