model/src/ini_cg2d.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cg2d.F,v 1.41 2003/01/13 01:31:22 jmc Exp $
C $Name:  $

#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"
c#include "DYNVARS.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     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 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
      _RL     faceArea
      _RS     myNorm
      _RL     sumArea
      _RL     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
          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
           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
       cg2dNorm = myNorm
      _BEGIN_MASTER( myThid )
CcnhDebugStarts
       WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ', 
     &               cg2dNorm
       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 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

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