model/src/ini_cg2d.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cg2d.F,v 1.44 2006/10/17 18:52:34 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: INI_CG2D
C     !INTERFACE:
      SUBROUTINE INI_CG2D( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_CG2D                                       
C     | o Initialise 2d 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 "CG2D.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

C     !LOCAL VARIABLES:
C     === Local variables ===
C     bi,bj  :: tile indices
C     I,J,K  :: Loop counters
C     faceArea - Temporary used to hold cell face areas.
C     myNorm - Work variable used in calculating normalisation factor
C     sumArea - Work variable used to compute the total Domain Area
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER bi, bj
      INTEGER I, J, K, ks
      _RL     faceArea
      _RS     myNorm
      _RS     aC, aCw, aCs
CEOP

C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
C     but safer when EXCH do not fill all the overlap regions.
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          aW2d(I,J,bi,bj) = 0. _d 0
          aS2d(I,J,bi,bj) = 0. _d 0
          aC2d(I,J,bi,bj) = 0. _d 0
          pW(I,J,bi,bj) = 0. _d 0
          pS(I,J,bi,bj) = 0. _d 0
          pC(I,J,bi,bj) = 0. _d 0
          cg2d_q(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
        DO J=1-1,sNy+1
         DO I=1-1,sNx+1
          cg2d_r(I,J,bi,bj) = 0. _d 0
          cg2d_s(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Initialise laplace operator
C     aW2d: integral in Z Ax/dX
C     aS2d: integral in Z Ay/dY
      myNorm = 0. _d 0
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          aW2d(I,J,bi,bj) = 0. _d 0
          aS2d(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
        DO K=1,Nr
         DO J=1,sNy
          DO I=1,sNx
C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
           faceArea = _dyG(I,J,bi,bj)*drF(K)
     &               *_hFacW(I,J,K,bi,bj)
           aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)
     &              + implicSurfPress*implicDiv2DFlow
     &               *faceArea*recip_dxC(I,J,bi,bj)
           faceArea = _dxG(I,J,bi,bj)*drF(K)
     &               *_hFacS(I,J,K,bi,bj)
           aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)
     &              + implicSurfPress*implicDiv2DFlow
     &               *faceArea*recip_dyC(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
         DO I=1,sNx
          IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
          IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
          IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
          IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
         ENDDO
         DO J=1,sNy
          IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
          IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
          IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
          IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
         ENDDO
        ENDIF
#endif
        DO J=1,sNy
         DO I=1,sNx
          myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
          myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
         ENDDO
        ENDDO
       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
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
          aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Update overlap regions
CcnhDebugStarts
C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
CcnhDebugEnds
c     _EXCH_XY_R4(aW2d, myThid)
c     _EXCH_XY_R4(aS2d, myThid)
      CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
CcnhDebugStarts
C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
CcnhDebugEnds

      _BEGIN_MASTER(myThid)
C-- set global parameter in common block:
      cg2dNorm = myNorm
C-- Define the solver tolerance in the appropriate Unit :
      cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
      IF (cg2dNormaliseRHS) THEN
C-  when using a normalisation of RHS, tolerance has no unit => no conversion
        cg2dTolerance = cg2dTargetResidual
      ELSE
C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
        cg2dTolerance = cg2dNorm * cg2dTargetResWunit
     &                           * globalArea / deltaTmom
      ENDIF
      _END_MASTER(myThid)

CcnhDebugStarts
      _BEGIN_MASTER( myThid )
       WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
     &      'CG2D normalisation factor = ', cg2dNorm
       CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
       IF (.NOT.cg2dNormaliseRHS) THEN
        WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
     &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
       ENDIF
       WRITE(msgBuf,*) ' '
       CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      _END_MASTER( myThid )
CcnhDebugEnds

C--   Initialise preconditioner
C     Note. 20th May 1998
C           I made a weird discovery! In the model paper we argue
C           for the form of the preconditioner used here ( see
C           A Finite-volume, Incompressible Navier-Stokes Model
C           ...., Marshall et. al ). The algebra gives a simple
C           0.5 factor for the averaging of ac and aCw to get a
C           symmettric pre-conditioner. By using a factor of 0.51
C           i.e. scaling the off-diagonal terms in the
C           preconditioner down slightly I managed to get the
C           number of iterations for convergence in a test case to
C           drop form 192 -> 134! Need to investigate this further!
C           For now I have introduced a parameter cg2dpcOffDFac which
C           defaults to 0.51 but can be set at runtime.
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          ks = ksurfC(I,J,bi,bj)
          pC(I,J,bi,bj) = 1. _d 0
          aC = -(
     &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
     &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
     &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks)
     &                *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
     &    )
          aCs = -(
     &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
     &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
     &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*deepFac2F(ks)
     &                *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
     &    )
          aCw = -(
     &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
     &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
     &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*deepFac2F(ks)
     &                *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
     &    )
          IF ( aC .EQ. 0. ) THEN
            pC(I,J,bi,bj) = 1. _d 0
          ELSE
           pC(I,J,bi,bj) =  1. _d 0 / aC
          ENDIF
          IF ( aC + aCw .EQ. 0. ) THEN
           pW(I,J,bi,bj) = 0.
          ELSE
           pW(I,J,bi,bj) =
     &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
          ENDIF
          IF ( aC + aCs .EQ. 0. ) THEN
           pS(I,J,bi,bj) = 0.
          ELSE
           pS(I,J,bi,bj) =
     &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
          ENDIF
C-    store solver main diagonal :
          aC2d(I,J,bi,bj) = aC
C         pC(I,J,bi,bj) = 1.
C         pW(I,J,bi,bj) = 0.
C         pS(I,J,bi,bj) = 0.
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C--   Update overlap regions
      _EXCH_XY_R4(pC, myThid)
c     _EXCH_XY_R4(pW, myThid)
c     _EXCH_XY_R4(pS, myThid)
      CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
CcnhDebugStarts
C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
CcnhDebugEnds

      RETURN
      END
1	jmc	1.45	C $Header: /u/gcmpack/MITgcm/model/src/ini_cg2d.F,v 1.44 2006/10/17 18:52:34 jmc Exp $
2	jmc	1.30	C $Name: $
3	cnh	1.1
4	edhill	1.42	#include "PACKAGES_CONFIG.h"
5	adcroft	1.13	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.37	CBOP
8			C !ROUTINE: INI_CG2D
9			C !INTERFACE:
10	cnh	1.1	SUBROUTINE INI_CG2D( myThid )
11	cnh	1.37
12			C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE INI_CG2D
15			C \| o Initialise 2d conjugate gradient solver operators.
16			C ==========================================================
17			C \| These arrays are purely a function of the basin geom.
18			C \| We set then here once and them use then repeatedly.
19			C ==========================================================
20			C \ev
21
22			C !USES:
23	adcroft	1.23	IMPLICIT NONE
24	cnh	1.1	C === Global variables ===
25			#include "SIZE.h"
26			#include "EEPARAMS.h"
27			#include "PARAMS.h"
28			#include "GRID.h"
29	jmc	1.31	#include "SURFACE.h"
30	adcroft	1.34	#include "CG2D.h"
31	adcroft	1.26	#ifdef ALLOW_OBCS
32	adcroft	1.22	#include "OBCS.h"
33	adcroft	1.26	#endif
34	cnh	1.1
35	cnh	1.37	C !INPUT/OUTPUT PARAMETERS:
36	cnh	1.1	C === Routine arguments ===
37			C myThid - Thread no. that called this routine.
38			INTEGER myThid
39
40	cnh	1.37	C !LOCAL VARIABLES:
41	cnh	1.1	C === Local variables ===
42	jmc	1.45	C bi,bj :: tile indices
43			C I,J,K :: Loop counters
44	cnh	1.1	C faceArea - Temporary used to hold cell face areas.
45	adcroft	1.34	C myNorm - Work variable used in calculating normalisation factor
46	jmc	1.45	C sumArea - Work variable used to compute the total Domain Area
47	cnh	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
48			INTEGER bi, bj
49	jmc	1.45	INTEGER I, J, K, ks
50	cnh	1.7	_RL faceArea
51	cnh	1.15	_RS myNorm
52	jmc	1.45	_RS aC, aCw, aCs
53	cnh	1.37	CEOP
54	jmc	1.38
55	jmc	1.45	C-- Initialize arrays in common blocs (CG2D.h) ; not really necessary
56	jmc	1.38	C but safer when EXCH do not fill all the overlap regions.
57			DO bj=myByLo(myThid),myByHi(myThid)
58			DO bi=myBxLo(myThid),myBxHi(myThid)
59			DO J=1-OLy,sNy+OLy
60			DO I=1-OLx,sNx+OLx
61			aW2d(I,J,bi,bj) = 0. _d 0
62			aS2d(I,J,bi,bj) = 0. _d 0
63	jmc	1.45	aC2d(I,J,bi,bj) = 0. _d 0
64	jmc	1.38	pW(I,J,bi,bj) = 0. _d 0
65			pS(I,J,bi,bj) = 0. _d 0
66			pC(I,J,bi,bj) = 0. _d 0
67			cg2d_q(I,J,bi,bj) = 0. _d 0
68			ENDDO
69			ENDDO
70			DO J=1-1,sNy+1
71			DO I=1-1,sNx+1
72			cg2d_r(I,J,bi,bj) = 0. _d 0
73			cg2d_s(I,J,bi,bj) = 0. _d 0
74			ENDDO
75			ENDDO
76			ENDDO
77			ENDDO
78	jmc	1.31
79	cnh	1.1	C-- Initialise laplace operator
80			C aW2d: integral in Z Ax/dX
81			C aS2d: integral in Z Ay/dY
82			myNorm = 0. _d 0
83			DO bj=myByLo(myThid),myByHi(myThid)
84			DO bi=myBxLo(myThid),myBxHi(myThid)
85			DO J=1,sNy
86			DO I=1,sNx
87			aW2d(I,J,bi,bj) = 0. _d 0
88			aS2d(I,J,bi,bj) = 0. _d 0
89			ENDDO
90			ENDDO
91	cnh	1.17	DO K=1,Nr
92	cnh	1.1	DO J=1,sNy
93			DO I=1,sNx
94	jmc	1.45	C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
95	cnh	1.20	faceArea = _dyG(I,J,bi,bj)*drF(K)
96			& *_hFacW(I,J,K,bi,bj)
97	jmc	1.45	aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)
98			& + implicSurfPress*implicDiv2DFlow
99			& faceArearecip_dxC(I,J,bi,bj)
100	cnh	1.20	faceArea = _dxG(I,J,bi,bj)*drF(K)
101			& *_hFacS(I,J,K,bi,bj)
102	jmc	1.45	aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)
103			& + implicSurfPress*implicDiv2DFlow
104			& faceArearecip_dyC(I,J,bi,bj)
105	cnh	1.1	ENDDO
106			ENDDO
107			ENDDO
108	adcroft	1.26	#ifdef ALLOW_OBCS
109	adcroft	1.28	IF (useOBCS) THEN
110	adcroft	1.22	DO I=1,sNx
111			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
112			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
113			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
114			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
115			ENDDO
116			DO J=1,sNy
117			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
118			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
119			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
120			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
121			ENDDO
122			ENDIF
123	adcroft	1.26	#endif
124	cnh	1.1	DO J=1,sNy
125			DO I=1,sNx
126			myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
127			myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
128			ENDDO
129			ENDDO
130			ENDDO
131			ENDDO
132	adcroft	1.25	_GLOBAL_MAX_R4( myNorm, myThid )
133			IF ( myNorm .NE. 0. _d 0 ) THEN
134			myNorm = 1. _d 0/myNorm
135	cnh	1.1	ELSE
136			myNorm = 1. _d 0
137			ENDIF
138			DO bj=myByLo(myThid),myByHi(myThid)
139			DO bi=myBxLo(myThid),myBxHi(myThid)
140			DO J=1,sNy
141			DO I=1,sNx
142			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
143			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
144			ENDDO
145			ENDDO
146			ENDDO
147			ENDDO
148	jmc	1.45
149	cnh	1.1	C-- Update overlap regions
150			CcnhDebugStarts
151	cnh	1.14	C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
152			C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
153	cnh	1.1	CcnhDebugEnds
154	adcroft	1.34	c _EXCH_XY_R4(aW2d, myThid)
155			c _EXCH_XY_R4(aS2d, myThid)
156			CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
157	cnh	1.1	CcnhDebugStarts
158	adcroft	1.24	C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
159			C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
160	cnh	1.1	CcnhDebugEnds
161
162	jmc	1.44	_BEGIN_MASTER(myThid)
163			C-- set global parameter in common block:
164			cg2dNorm = myNorm
165	adcroft	1.34	C-- Define the solver tolerance in the appropriate Unit :
166	jmc	1.43	cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
167	adcroft	1.34	IF (cg2dNormaliseRHS) THEN
168			C- when using a normalisation of RHS, tolerance has no unit => no conversion
169			cg2dTolerance = cg2dTargetResidual
170			ELSE
171			C- convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
172			cg2dTolerance = cg2dNorm * cg2dTargetResWunit
173	jmc	1.43	& * globalArea / deltaTmom
174	adcroft	1.34	ENDIF
175	jmc	1.44	_END_MASTER(myThid)
176	jmc	1.43
177			CcnhDebugStarts
178			_BEGIN_MASTER( myThid )
179			WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
180			& 'CG2D normalisation factor = ', cg2dNorm
181			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
182			IF (.NOT.cg2dNormaliseRHS) THEN
183			WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
184			& 'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
185			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
186			ENDIF
187			WRITE(msgBuf,*) ' '
188			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
189			_END_MASTER( myThid )
190			CcnhDebugEnds
191	jmc	1.45
192	cnh	1.1	C-- Initialise preconditioner
193	cnh	1.4	C Note. 20th May 1998
194			C I made a weird discovery! In the model paper we argue
195			C for the form of the preconditioner used here ( see
196			C A Finite-volume, Incompressible Navier-Stokes Model
197			C ...., Marshall et. al ). The algebra gives a simple
198			C 0.5 factor for the averaging of ac and aCw to get a
199			C symmettric pre-conditioner. By using a factor of 0.51
200			C i.e. scaling the off-diagonal terms in the
201			C preconditioner down slightly I managed to get the
202			C number of iterations for convergence in a test case to
203			C drop form 192 -> 134! Need to investigate this further!
204			C For now I have introduced a parameter cg2dpcOffDFac which
205			C defaults to 0.51 but can be set at runtime.
206	cnh	1.1	DO bj=myByLo(myThid),myByHi(myThid)
207			DO bi=myBxLo(myThid),myBxHi(myThid)
208			DO J=1,sNy
209			DO I=1,sNx
210	jmc	1.45	ks = ksurfC(I,J,bi,bj)
211	cnh	1.1	pC(I,J,bi,bj) = 1. _d 0
212	cnh	1.4	aC = -(
213			& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
214	adcroft	1.24	& +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
215	jmc	1.45	& +freeSurfFacmyNormrecip_Bo(I,J,bi,bj)*deepFac2F(ks)
216			& *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
217	cnh	1.4	& )
218			aCs = -(
219			& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
220			& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
221	jmc	1.45	& +freeSurfFacmyNormrecip_Bo(I,J-1,bi,bj)*deepFac2F(ks)
222			& *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
223	cnh	1.4	& )
224			aCw = -(
225			& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
226			& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
227	jmc	1.45	& +freeSurfFacmyNormrecip_Bo(I-1,J,bi,bj)*deepFac2F(ks)
228			& *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
229	cnh	1.4	& )
230			IF ( aC .EQ. 0. ) THEN
231	adcroft	1.27	pC(I,J,bi,bj) = 1. _d 0
232	cnh	1.4	ELSE
233			pC(I,J,bi,bj) = 1. _d 0 / aC
234			ENDIF
235			IF ( aC + aCw .EQ. 0. ) THEN
236			pW(I,J,bi,bj) = 0.
237			ELSE
238	jmc	1.45	pW(I,J,bi,bj) =
239	cnh	1.4	& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCw+aC))*2 )
240			ENDIF
241			IF ( aC + aCs .EQ. 0. ) THEN
242			pS(I,J,bi,bj) = 0.
243			ELSE
244			pS(I,J,bi,bj) =
245			& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCs+aC))*2 )
246			ENDIF
247	jmc	1.45	C- store solver main diagonal :
248			aC2d(I,J,bi,bj) = aC
249	cnh	1.6	C pC(I,J,bi,bj) = 1.
250			C pW(I,J,bi,bj) = 0.
251			C pS(I,J,bi,bj) = 0.
252	cnh	1.1	ENDDO
253			ENDDO
254			ENDDO
255			ENDDO
256			C-- Update overlap regions
257			_EXCH_XY_R4(pC, myThid)
258	adcroft	1.34	c _EXCH_XY_R4(pW, myThid)
259			c _EXCH_XY_R4(pS, myThid)
260			CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
261	cnh	1.18	CcnhDebugStarts
262	adcroft	1.24	C CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid )
263			C CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid )
264			C CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid )
265	cnh	1.18	CcnhDebugEnds
266	cnh	1.1
267			RETURN
268			END