model/src/ini_cg3d.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.22 2006/12/05 05:25:08 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_CG3D
C     !INTERFACE:
      SUBROUTINE INI_CG3D( myThid )
C     !DESCRIPTION: \bv
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     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "CG3D.h"
#include "SOLVE_FOR_PRESSURE3D.h"
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC

C     !LOCAL VARIABLES:
C     === Local variables ===
C     bi,bj  - Loop counters
C     I,J,K  - Loop counters
C     faceArea - Temporary used to hold cell face areas.
C     myNorm - Work variable used in clculating normalisation factor
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER bi, bj
      INTEGER I, J, K, ks
      _RL     faceArea
      _RS     myNorm
      _RL     theRecip_Dr
      _RL     aU, aL, aW, aE, aN, aS, aC
      _RL     tmpFac, nh_Fac, igwFac
CEOP

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

C--   Initialise to zero over the full range of indices
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=1,Nr
C-       From common bloc CG3D_R:
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           aW3d(I,J,K,bi,bj) = 0.
           aS3d(I,J,K,bi,bj) = 0.
           aV3d(I,J,K,bi,bj) = 0.
           aC3d(I,J,K,bi,bj) = 0.
           zMC (I,J,K,bi,bj) = 0.
           zML (I,J,K,bi,bj) = 0.
           zMU (I,J,K,bi,bj) = 0.
          ENDDO
         ENDDO
C-       From common bloc CG3D_WK_R:
         DO J=0,sNy+1
          DO I=0,sNx+1
           cg3d_q(I,J,K,bi,bj) = 0.
           cg3d_r(I,J,K,bi,bj) = 0.
           cg3d_s(I,J,K,bi,bj) = 0.
          ENDDO
         ENDDO
C-       From common bloc SFP3D_COMMON_R:
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           cg3d_b(I,J,K,bi,bj) = 0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      nh_Fac = 0.
      igwFac = 0.
      IF ( nonHydrostatic
     &      .AND. nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2
      IF ( implicitIntGravWave ) igwFac = 1. _d 0

      IF ( use3Dsolver ) THEN
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+1
           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)
           myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
          ENDDO
         ENDDO
C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
         DO J=1,sNy+1
          DO I=1,sNx
           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(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
         tmpFac = nh_Fac*rVel2wUnit(k)*rVel2wUnit(k)
     &          + igwFac*dBdrRef(k)*deltaTMom*dTtracerLev(k)
         IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
         DO J=1,sNy
          DO I=1,sNx
           faceArea = _rA(I,J,bi,bj)*maskC(I,J, K ,bi,bj)
     &                              *maskC(I,J,K-1,bi,bj)
     &                              *deepFac2F(k)
           theRecip_Dr = recip_drC(K)
c          theRecip_Dr =
caja &      drF(K  )*_hFacC(i,j,k  ,bi,bj)*0.5
caja &     +drF(K-1)*_hFacC(i,j,k-1,bi,bj)*0.5
c          IF ( theRecip_Dr .NE. 0. )
c    &      theRecip_Dr = 1. _d 0/theRecip_Dr
           aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr*tmpFac
           myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_OBCS
        IF ( useOBCS ) THEN
         DO K=1,Nr
          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
         ENDDO
        ENDIF
#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

      _BEGIN_MASTER( myThid )
C-- set global parameter in common block:
       cg3dNorm = myNorm
CcnhDebugStarts
       WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG3D: ',
     &      '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)
C-    Set solver main diagonal term
        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)
           aU = aV3d(I,J,K,bi,bj)
           IF ( K .NE. Nr  ) THEN
            aL = aV3d(I,J,K+1,bi,bj)
           ELSE
            aL = 0.
           ENDIF
           aC3d(I,J,K,bi,bj) = -aW-aE-aN-aS-aU-aL
          ENDDO
         ENDDO
        ENDDO
C-    Add free-surface source term
        DO J=1,sNy
         DO I=1,sNx
           ks = ksurfC(I,J,bi,bj)
           IF ( ks.LE.Nr ) THEN
             aC3d(I,J,ks,bi,bj) = aC3d(I,J,ks,bi,bj)
     &         -freeSurfFac*recip_Bo(I,J,bi,bj)
     &          *rA(I,J,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTfreesurf
           ENDIF
         ENDDO
        ENDDO
C-    Matrix solver normalisation
        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
           aC3d(I,J,K,bi,bj) = aC3d(I,J,K,bi,bj)*myNorm
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Update overlap regions
c     _EXCH_XYZ_R4(aW3d, myThid)
c     _EXCH_XYZ_R4(aS3d, myThid)
      CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
      _EXCH_XYZ_R4(aV3d, myThid)
      _EXCH_XYZ_R4(aC3d, 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
           IF ( aC3d(I,J,K,bi,bj) .NE. 0. ) THEN
            zMC(i,j,k,bi,bj) = aC3d(I,J,K,bi,bj)
            zML(i,j,k,bi,bj) = aV3d(I,J,K,bi,bj)
            zMU(i,j,k,bi,bj) = 0.
            IF ( K.NE.Nr )
     &      zMU(i,j,k,bi,bj) = aV3d(I,J,K+1,bi,bj)
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
c     DO k=1,Nr
c     DO j=1-OLy,sNy+OLy
c     DO i=1-OLx,sNx+OLx
c      phi(i,j,1,1) = zMc(i,j,k,1,1)
c     ENDDO
c     ENDDO
C     CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
c     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--   end if (use3Dsolver)
      ENDIF

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	jmc	1.23	C $Header: /u/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.22 2006/12/05 05:25:08 jmc Exp $
2	jmc	1.10	C $Name: $
3	adcroft	1.1
4	edhill	1.15	#include "PACKAGES_CONFIG.h"
5	adcroft	1.1	#include "CPP_OPTIONS.h"
6
7	cnh	1.13	CBOP
8	jmc	1.18	C !ROUTINE: INI_CG3D
9	cnh	1.13	C !INTERFACE:
10	adcroft	1.1	SUBROUTINE INI_CG3D( myThid )
11	cnh	1.13	C !DESCRIPTION: \bv
12			C ==========================================================
13			C \| SUBROUTINE INI_CG3D
14			C \| o Initialise 3d conjugate gradient solver operators.
15			C ==========================================================
16			C \| These arrays are purely a function of the basin geom.
17			C \| We set then here once and them use then repeatedly.
18			C ==========================================================
19			C \ev
20
21			C !USES:
22	adcroft	1.3	IMPLICIT NONE
23	adcroft	1.1	C === Global variables ===
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "GRID.h"
28	jmc	1.19	#include "SURFACE.h"
29	adcroft	1.1	#include "CG3D.h"
30	adcroft	1.12	#include "SOLVE_FOR_PRESSURE3D.h"
31	adcroft	1.4	#ifdef ALLOW_OBCS
32	adcroft	1.1	#include "OBCS.h"
33	adcroft	1.4	#endif
34	adcroft	1.1
35	cnh	1.13	C !INPUT/OUTPUT PARAMETERS:
36	adcroft	1.1	C === Routine arguments ===
37			C myThid - Thread no. that called this routine.
38			INTEGER myThid
39
40	adcroft	1.6	#ifdef ALLOW_NONHYDROSTATIC
41
42	cnh	1.13	C !LOCAL VARIABLES:
43	adcroft	1.1	C === Local variables ===
44			C bi,bj - Loop counters
45	jmc	1.18	C I,J,K - Loop counters
46	adcroft	1.1	C faceArea - Temporary used to hold cell face areas.
47			C myNorm - Work variable used in clculating normalisation factor
48			CHARACTER*(MAX_LEN_MBUF) msgBuf
49			INTEGER bi, bj
50	jmc	1.20	INTEGER I, J, K, ks
51	adcroft	1.1	_RL faceArea
52			_RS myNorm
53			_RL theRecip_Dr
54			_RL aU, aL, aW, aE, aN, aS, aC
55	jmc	1.20	_RL tmpFac, nh_Fac, igwFac
56	cnh	1.13	CEOP
57	adcroft	1.1
58			CcnhDebugStarts
59	jmc	1.18	c _RL phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
60	adcroft	1.1	CcnhDebugEnds
61
62	jmc	1.18	C-- Initialise to zero over the full range of indices
63			DO bj=myByLo(myThid),myByHi(myThid)
64			DO bi=myBxLo(myThid),myBxHi(myThid)
65			DO K=1,Nr
66			C- From common bloc CG3D_R:
67			DO J=1-OLy,sNy+OLy
68			DO I=1-OLx,sNx+OLx
69			aW3d(I,J,K,bi,bj) = 0.
70			aS3d(I,J,K,bi,bj) = 0.
71			aV3d(I,J,K,bi,bj) = 0.
72	jmc	1.19	aC3d(I,J,K,bi,bj) = 0.
73	jmc	1.18	zMC (I,J,K,bi,bj) = 0.
74			zML (I,J,K,bi,bj) = 0.
75			zMU (I,J,K,bi,bj) = 0.
76			ENDDO
77			ENDDO
78			C- From common bloc CG3D_WK_R:
79			DO J=0,sNy+1
80			DO I=0,sNx+1
81			cg3d_q(I,J,K,bi,bj) = 0.
82			cg3d_r(I,J,K,bi,bj) = 0.
83			cg3d_s(I,J,K,bi,bj) = 0.
84			ENDDO
85			ENDDO
86			C- From common bloc SFP3D_COMMON_R:
87			DO J=1-OLy,sNy+OLy
88			DO I=1-OLx,sNx+OLx
89			cg3d_b(I,J,K,bi,bj) = 0.
90			ENDDO
91			ENDDO
92			ENDDO
93			ENDDO
94			ENDDO
95
96	jmc	1.20	nh_Fac = 0.
97			igwFac = 0.
98			IF ( nonHydrostatic
99			& .AND. nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2
100			IF ( implicitIntGravWave ) igwFac = 1. _d 0
101
102			IF ( use3Dsolver ) THEN
103	adcroft	1.1	C-- Initialise laplace operator
104			C aW3d: Ax/dX
105			C aS3d: Ay/dY
106			C aV3d: Ar/dR
107			myNorm = 0. _d 0
108			DO bj=myByLo(myThid),myByHi(myThid)
109			DO bi=myBxLo(myThid),myBxHi(myThid)
110			DO K=1,Nr
111			DO J=1,sNy
112	jmc	1.19	DO I=1,sNx+1
113	adcroft	1.1	faceArea = _dyG(I,J,bi,bj)*drF(K)
114			& *_hFacW(I,J,K,bi,bj)
115	jmc	1.10	aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
116	jmc	1.19	myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
117			ENDDO
118			ENDDO
119	jmc	1.22	C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
120	jmc	1.19	DO J=1,sNy+1
121			DO I=1,sNx
122	adcroft	1.1	faceArea = _dxG(I,J,bi,bj)*drF(K)
123			& *_hFacS(I,J,K,bi,bj)
124	jmc	1.10	aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
125	adcroft	1.1	myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
126			ENDDO
127			ENDDO
128			ENDDO
129			DO K=1,1
130			DO J=1,sNy
131			DO I=1,sNx
132			aV3d(I,J,K,bi,bj) = 0.
133			myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
134			ENDDO
135			ENDDO
136			ENDDO
137			DO K=2,Nr
138	jmc	1.23	tmpFac = nh_FacrVel2wUnit(k)rVel2wUnit(k)
139	jmc	1.20	& + igwFacdBdrRef(k)deltaTMom*dTtracerLev(k)
140			IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
141	adcroft	1.1	DO J=1,sNy
142			DO I=1,sNx
143	jmc	1.18	faceArea = _rA(I,J,bi,bj)*maskC(I,J, K ,bi,bj)
144			& *maskC(I,J,K-1,bi,bj)
145	jmc	1.22	& *deepFac2F(k)
146	jmc	1.18	theRecip_Dr = recip_drC(K)
147	jmc	1.22	c theRecip_Dr =
148	adcroft	1.1	caja & drF(K )_hFacC(i,j,k ,bi,bj)0.5
149			caja & +drF(K-1)_hFacC(i,j,k-1,bi,bj)0.5
150	jmc	1.22	c IF ( theRecip_Dr .NE. 0. )
151	jmc	1.18	c & theRecip_Dr = 1. _d 0/theRecip_Dr
152	jmc	1.20	aV3d(I,J,K,bi,bj) = faceAreatheRecip_DrtmpFac
153	adcroft	1.1	myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
154			ENDDO
155			ENDDO
156	adcroft	1.5	ENDDO
157	adcroft	1.4	#ifdef ALLOW_OBCS
158	jmc	1.18	IF ( useOBCS ) THEN
159			DO K=1,Nr
160	adcroft	1.1	DO I=1,sNx
161			IF (OB_Jn(I,bi,bj).NE.0) THEN
162			aS3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
163			aS3d(I,OB_Jn(I,bi,bj)+1,K,bi,bj)=0.
164			aW3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
165			aW3d(I+1,OB_Jn(I,bi,bj),K,bi,bj)=0.
166	adcroft	1.7	aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
167	adcroft	1.1	ENDIF
168			IF (OB_Js(I,bi,bj).NE.0) THEN
169			aS3d(I,OB_Js(I,bi,bj)+1,K,bi,bj)=0.
170			aS3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
171			aW3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
172			aW3d(I+1,OB_Js(I,bi,bj),K,bi,bj)=0.
173	adcroft	1.7	aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
174	adcroft	1.1	ENDIF
175			ENDDO
176			DO J=1,sNy
177			IF (OB_Ie(J,bi,bj).NE.0) THEN
178			aW3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
179			aW3d(OB_Ie(J,bi,bj)+1,J,K,bi,bj)=0.
180			aS3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
181			aS3d(OB_Ie(J,bi,bj),J+1,K,bi,bj)=0.
182	adcroft	1.8	aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
183	adcroft	1.1	ENDIF
184			IF (OB_Iw(J,bi,bj).NE.0) THEN
185			aW3d(OB_Iw(J,bi,bj)+1,J,K,bi,bj)=0.
186			aW3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
187			aS3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
188			aS3d(OB_Iw(J,bi,bj),J+1,K,bi,bj)=0.
189	adcroft	1.8	aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
190	adcroft	1.1	ENDIF
191			ENDDO
192	jmc	1.18	ENDDO
193			ENDIF
194	adcroft	1.4	#endif
195	adcroft	1.1	ENDDO
196			ENDDO
197	adcroft	1.2	_GLOBAL_MAX_R4( myNorm, myThid )
198			IF ( myNorm .NE. 0. _d 0 ) THEN
199			myNorm = 1. _d 0/myNorm
200	adcroft	1.1	ELSE
201			myNorm = 1. _d 0
202			ENDIF
203	jmc	1.22
204	jmc	1.21	_BEGIN_MASTER( myThid )
205	jmc	1.22	C-- set global parameter in common block:
206	adcroft	1.1	cg3dNorm = myNorm
207			CcnhDebugStarts
208	jmc	1.22	WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG3D: ',
209			& 'CG3D normalisation factor = ', cg3dNorm
210	adcroft	1.1	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
211			WRITE(msgBuf,*) ' '
212			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
213			CcnhDebugEnds
214			_END_MASTER( myThid )
215	jmc	1.22
216	adcroft	1.1	DO bj=myByLo(myThid),myByHi(myThid)
217			DO bi=myBxLo(myThid),myBxHi(myThid)
218	jmc	1.19	C- Set solver main diagonal term
219			DO K=1,Nr
220			DO J=1,sNy
221			DO I=1,sNx
222			aW = aW3d(I ,J,K,bi,bj)
223			aE = aW3d(I+1,J,K,bi,bj)
224			aN = aS3d(I,J+1,K,bi,bj)
225			aS = aS3d(I,J ,K,bi,bj)
226			aU = aV3d(I,J,K,bi,bj)
227			IF ( K .NE. Nr ) THEN
228			aL = aV3d(I,J,K+1,bi,bj)
229			ELSE
230			aL = 0.
231			ENDIF
232			aC3d(I,J,K,bi,bj) = -aW-aE-aN-aS-aU-aL
233			ENDDO
234			ENDDO
235			ENDDO
236			C- Add free-surface source term
237			DO J=1,sNy
238			DO I=1,sNx
239			ks = ksurfC(I,J,bi,bj)
240			IF ( ks.LE.Nr ) THEN
241			aC3d(I,J,ks,bi,bj) = aC3d(I,J,ks,bi,bj)
242			& -freeSurfFac*recip_Bo(I,J,bi,bj)
243	jmc	1.22	& rA(I,J,bi,bj)deepFac2F(ks)/deltaTMom/deltaTfreesurf
244	jmc	1.19	ENDIF
245			ENDDO
246			ENDDO
247			C- Matrix solver normalisation
248	adcroft	1.1	DO K=1,Nr
249			DO J=1,sNy
250			DO I=1,sNx
251			aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
252			aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
253			aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
254	jmc	1.19	aC3d(I,J,K,bi,bj) = aC3d(I,J,K,bi,bj)*myNorm
255	adcroft	1.1	ENDDO
256			ENDDO
257			ENDDO
258			ENDDO
259			ENDDO
260	jmc	1.19
261	adcroft	1.1	C-- Update overlap regions
262	jmc	1.18	c _EXCH_XYZ_R4(aW3d, myThid)
263			c _EXCH_XYZ_R4(aS3d, myThid)
264			CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
265	adcroft	1.1	_EXCH_XYZ_R4(aV3d, myThid)
266	jmc	1.19	_EXCH_XYZ_R4(aC3d, myThid)
267	adcroft	1.1	CcnhDebugStarts
268			C CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
269			C CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
270			CcnhDebugEnds
271
272			C-- Initialise preconditioner
273			C For now PC is just the identity. Change to
274			C be LU factorization of d2/dz2 later. Note
275			C check for consistency with S/R CG3D before
276			C assuming zML is lower and zMU is upper!
277			DO bj=myByLo(myThid),myByHi(myThid)
278			DO bi=myBxLo(myThid),myBxHi(myThid)
279	jmc	1.19	DO K=1,Nr
280	adcroft	1.1	DO J=1,sNy
281			DO I=1,sNx
282	jmc	1.19	IF ( aC3d(I,J,K,bi,bj) .NE. 0. ) THEN
283			zMC(i,j,k,bi,bj) = aC3d(I,J,K,bi,bj)
284			zML(i,j,k,bi,bj) = aV3d(I,J,K,bi,bj)
285			zMU(i,j,k,bi,bj) = 0.
286			IF ( K.NE.Nr )
287			& zMU(i,j,k,bi,bj) = aV3d(I,J,K+1,bi,bj)
288	adcroft	1.1	CcnhDebugStarts
289			C zMC(i,j,k,bi,bj) = 1.
290			C zMU(i,j,k,bi,bj) = 0.
291			C zML(i,j,k,bi,bj) = 0.
292			CcnhDebugEnds
293			ELSE
294			zMC(i,j,k,bi,bj) = 1. _d 0
295			zMU(i,j,k,bi,bj) = 0.
296			zML(i,j,k,bi,bj) = 0.
297			ENDIF
298			ENDDO
299			ENDDO
300			ENDDO
301			DO J=1,sNy
302			DO I=1,sNx
303			zMC(i,j,1,bi,bj)=
304			& 1. _d 0 / zMC(i,j,1,bi,bj)
305			zMU(i,j,1,bi,bj)=
306			& zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
307			ENDDO
308			ENDDO
309	jmc	1.19	DO K=2,Nr
310	adcroft	1.1	DO J=1,sNy
311			DO I=1,sNx
312			zMC(i,j,k,bi,bj) = 1. _d 0 /
313			& (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
314			zMU(i,j,k,bi,bj)=zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
315			ENDDO
316			ENDDO
317			ENDDO
318	jmc	1.19	DO K=1,Nr
319	adcroft	1.1	DO J=1,sNy
320			DO I=1,sNx
321			aW = aW3d(I ,J,K,bi,bj)
322			aE = aW3d(I+1,J,K,bi,bj)
323			aN = aS3d(I,J+1,K,bi,bj)
324			aS = aS3d(I,J ,K,bi,bj)
325			IF ( K .NE. 1 ) THEN
326			aU = aV3d(I,J,K-1,bi,bj)
327			ELSE
328	jmc	1.19	aU = 0.
329	adcroft	1.1	ENDIF
330			IF ( K .NE. Nr ) THEN
331			aL = aV3d(I,J,K+1,bi,bj)
332			ELSE
333	jmc	1.19	aL = 0.
334	adcroft	1.1	ENDIF
335			aC = -aW-aE-aN-aS-aU-aL
336			IF ( aC .EQ. 0. ) THEN
337	adcroft	1.5	zMC(i,j,k,bi,bj) = 1.
338	adcroft	1.1	zML(i,j,k,bi,bj) = 0.
339			zMU(i,j,k,bi,bj) = 0.
340			CcnhDebugStarts
341			C ELSE
342			C zMC(i,j,k,bi,bj) = 1.
343			C zML(i,j,k,bi,bj) = 0.
344			C zMU(i,j,k,bi,bj) = 0.
345			CcnhDEbugEnds
346			ENDIF
347			ENDDO
348			ENDDO
349			ENDDO
350			ENDDO
351			ENDDO
352			C-- Update overlap regions
353			_EXCH_XYZ_R4(zMC, myThid)
354			_EXCH_XYZ_R4(zML, myThid)
355			_EXCH_XYZ_R4(zMU, myThid)
356
357			CcnhDebugStarts
358	jmc	1.18	c DO k=1,Nr
359			c DO j=1-OLy,sNy+OLy
360			c DO i=1-OLx,sNx+OLx
361			c phi(i,j,1,1) = zMc(i,j,k,1,1)
362			c ENDDO
363			c ENDDO
364	adcroft	1.1	C CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
365	jmc	1.18	c ENDDO
366	adcroft	1.1	C CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
367			C CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
368			CcnhDebugEnds
369
370	jmc	1.20	C-- end if (use3Dsolver)
371			ENDIF
372
373	adcroft	1.1	#endif /* ALLOW_NONHYDROSTATIC */
374
375			RETURN
376			END