model/src/ini_cg2d.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.30 2001/02/20 15:06:21 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE INI_CG2D( myThid )
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     \==========================================================/
      IMPLICIT NONE

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

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

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     aC, aCw, aCs

C-- Initialise -Boyancy at surface level : Bo_surf
C   Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)
C     with rtoz = conversion factor from r-unit to z-unit  (=horiVertRatio)
C   an acurate formulation will include P_surf and T,S_surf effects on rho_surf:
C    z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0 
C    p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
C--
      IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
             Bo_surf(I,J,bi,bj) = recip_rhoConst
             recip_Bo(I,J,bi,bj) = rhoConst
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ELSE
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
             Bo_surf(I,J,bi,bj) = gravity
             recip_Bo(I,J,bi,bj) = recip_Gravity
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ENDIF

C--   Update overlap regions
      _EXCH_XY_R8(Bo_surf, myThid)
      _EXCH_XY_R8(recip_Bo, myThid)

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*
     &      gBaro*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*
     &      gBaro*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
      _EXCH_XY_R4(aW2d, myThid)
      _EXCH_XY_R4(aS2d, 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--   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* horiVertRatio*
     &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
     &    )
          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* horiVertRatio*
     &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
     &    )
          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* horiVertRatio*
     &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
     &    )
          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)
      _EXCH_XY_R4(pW, myThid)
      _EXCH_XY_R4(pS, 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.31	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.30 2001/02/20 15:06:21 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			CStartOfInterface
7			SUBROUTINE INI_CG2D( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE INI_CG2D \|
10			C \| o Initialise 2d 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.23	IMPLICIT NONE
16	cnh	1.1
17			C === Global variables ===
18			#include "SIZE.h"
19			#include "EEPARAMS.h"
20			#include "PARAMS.h"
21			#include "GRID.h"
22	cnh	1.12	#include "DYNVARS.h"
23	jmc	1.31	#include "SURFACE.h"
24			#include "CG2D_INTERNAL.h"
25	adcroft	1.26	#ifdef ALLOW_OBCS
26	adcroft	1.22	#include "OBCS.h"
27	adcroft	1.26	#endif
28	cnh	1.1
29			C === Routine arguments ===
30			C myThid - Thread no. that called this routine.
31			INTEGER myThid
32			CEndOfInterface
33
34			C === Local variables ===
35			C xG, yG - Global coordinate location.
36			C zG
37			C iG, jG - Global coordinate index
38			C bi,bj - Loop counters
39			C faceArea - Temporary used to hold cell face areas.
40			C I,J,K
41			C myNorm - Work variable used in clculating normalisation factor
42			CHARACTER*(MAX_LEN_MBUF) msgBuf
43			INTEGER bi, bj
44			INTEGER I, J, K
45	cnh	1.7	_RL faceArea
46	cnh	1.15	_RS myNorm
47	cnh	1.4	_RL aC, aCw, aCs
48	cnh	1.1
49	jmc	1.31	C-- Initialise -Boyancy at surface level : Bo_surf
50			C Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)
51			C with rtoz = conversion factor from r-unit to z-unit (=horiVertRatio)
52			C an acurate formulation will include P_surf and T,S_surf effects on rho_surf:
53			C z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0
54			C p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
55			C--
56			IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
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			Bo_surf(I,J,bi,bj) = recip_rhoConst
62			recip_Bo(I,J,bi,bj) = rhoConst
63			ENDDO
64			ENDDO
65			ENDDO
66			ENDDO
67			ELSE
68			DO bj=myByLo(myThid),myByHi(myThid)
69			DO bi=myBxLo(myThid),myBxHi(myThid)
70			DO J=1-Oly,sNy+Oly
71			DO I=1-Olx,sNx+Olx
72			Bo_surf(I,J,bi,bj) = gravity
73			recip_Bo(I,J,bi,bj) = recip_Gravity
74			ENDDO
75			ENDDO
76			ENDDO
77			ENDDO
78			ENDIF
79
80			C-- Update overlap regions
81			_EXCH_XY_R8(Bo_surf, myThid)
82			_EXCH_XY_R8(recip_Bo, myThid)
83
84	cnh	1.1	C-- Initialise laplace operator
85			C aW2d: integral in Z Ax/dX
86			C aS2d: integral in Z Ay/dY
87			myNorm = 0. _d 0
88			DO bj=myByLo(myThid),myByHi(myThid)
89			DO bi=myBxLo(myThid),myBxHi(myThid)
90			DO J=1,sNy
91			DO I=1,sNx
92			aW2d(I,J,bi,bj) = 0. _d 0
93			aS2d(I,J,bi,bj) = 0. _d 0
94			ENDDO
95			ENDDO
96	cnh	1.17	DO K=1,Nr
97	cnh	1.1	DO J=1,sNy
98			DO I=1,sNx
99	cnh	1.20	faceArea = _dyG(I,J,bi,bj)*drF(K)
100			& *_hFacW(I,J,K,bi,bj)
101	cnh	1.1	aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
102	jmc	1.30	& implicSurfPressimplicDiv2DFlow
103	cnh	1.17	& gBarofaceArearecip_dxC(I,J,bi,bj)
104	cnh	1.20	faceArea = _dxG(I,J,bi,bj)*drF(K)
105			& *_hFacS(I,J,K,bi,bj)
106	cnh	1.1	aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
107	jmc	1.30	& implicSurfPressimplicDiv2DFlow
108	cnh	1.17	& gBarofaceArearecip_dyC(I,J,bi,bj)
109	cnh	1.1	ENDDO
110			ENDDO
111			ENDDO
112	adcroft	1.26	#ifdef ALLOW_OBCS
113	adcroft	1.28	IF (useOBCS) THEN
114	adcroft	1.22	DO I=1,sNx
115			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
116			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
117			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
118			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
119			ENDDO
120			DO J=1,sNy
121			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
122			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
123			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
124			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
125			ENDDO
126			ENDIF
127	adcroft	1.26	#endif
128	cnh	1.1	DO J=1,sNy
129			DO I=1,sNx
130			myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
131			myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
132			ENDDO
133			ENDDO
134			ENDDO
135			ENDDO
136	adcroft	1.25	_GLOBAL_MAX_R4( myNorm, myThid )
137			IF ( myNorm .NE. 0. _d 0 ) THEN
138			myNorm = 1. _d 0/myNorm
139	cnh	1.1	ELSE
140			myNorm = 1. _d 0
141			ENDIF
142	cnh	1.5	cg2dNorm = myNorm
143	cnh	1.1	_BEGIN_MASTER( myThid )
144			CcnhDebugStarts
145	cnh	1.20	WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
146			& cg2dNorm
147	cnh	1.3	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
148			WRITE(msgBuf,*) ' '
149	cnh	1.1	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
150			CcnhDebugEnds
151			_END_MASTER( myThid )
152			DO bj=myByLo(myThid),myByHi(myThid)
153			DO bi=myBxLo(myThid),myBxHi(myThid)
154			DO J=1,sNy
155			DO I=1,sNx
156			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
157			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
158			ENDDO
159			ENDDO
160			ENDDO
161			ENDDO
162
163			C-- Update overlap regions
164			CcnhDebugStarts
165	cnh	1.14	C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
166			C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
167	cnh	1.1	CcnhDebugEnds
168			_EXCH_XY_R4(aW2d, myThid)
169			_EXCH_XY_R4(aS2d, myThid)
170			CcnhDebugStarts
171	adcroft	1.24	C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
172			C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
173	cnh	1.1	CcnhDebugEnds
174
175			C-- Initialise preconditioner
176	cnh	1.4	C Note. 20th May 1998
177			C I made a weird discovery! In the model paper we argue
178			C for the form of the preconditioner used here ( see
179			C A Finite-volume, Incompressible Navier-Stokes Model
180			C ...., Marshall et. al ). The algebra gives a simple
181			C 0.5 factor for the averaging of ac and aCw to get a
182			C symmettric pre-conditioner. By using a factor of 0.51
183			C i.e. scaling the off-diagonal terms in the
184			C preconditioner down slightly I managed to get the
185			C number of iterations for convergence in a test case to
186			C drop form 192 -> 134! Need to investigate this further!
187			C For now I have introduced a parameter cg2dpcOffDFac which
188			C defaults to 0.51 but can be set at runtime.
189	cnh	1.1	DO bj=myByLo(myThid),myByHi(myThid)
190			DO bi=myBxLo(myThid),myBxHi(myThid)
191			DO J=1,sNy
192			DO I=1,sNx
193			pC(I,J,bi,bj) = 1. _d 0
194	cnh	1.4	aC = -(
195			& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
196	adcroft	1.24	& +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
197	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
198	cnh	1.17	& rA(I,J,bi,bj)/deltaTMom/deltaTMom
199	cnh	1.4	& )
200			aCs = -(
201			& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
202			& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
203	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
204	cnh	1.17	& rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
205	cnh	1.4	& )
206			aCw = -(
207			& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
208			& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
209	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
210	cnh	1.17	& rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
211	cnh	1.4	& )
212			IF ( aC .EQ. 0. ) THEN
213	adcroft	1.27	pC(I,J,bi,bj) = 1. _d 0
214	cnh	1.4	ELSE
215			pC(I,J,bi,bj) = 1. _d 0 / aC
216			ENDIF
217			IF ( aC + aCw .EQ. 0. ) THEN
218			pW(I,J,bi,bj) = 0.
219			ELSE
220			pW(I,J,bi,bj) =
221			& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCw+aC))*2 )
222			ENDIF
223			IF ( aC + aCs .EQ. 0. ) THEN
224			pS(I,J,bi,bj) = 0.
225			ELSE
226			pS(I,J,bi,bj) =
227			& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCs+aC))*2 )
228			ENDIF
229	cnh	1.6	C pC(I,J,bi,bj) = 1.
230			C pW(I,J,bi,bj) = 0.
231			C pS(I,J,bi,bj) = 0.
232	cnh	1.1	ENDDO
233			ENDDO
234			ENDDO
235			ENDDO
236			C-- Update overlap regions
237			_EXCH_XY_R4(pC, myThid)
238			_EXCH_XY_R4(pW, myThid)
239			_EXCH_XY_R4(pS, myThid)
240	cnh	1.18	CcnhDebugStarts
241	adcroft	1.24	C CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid )
242			C CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid )
243			C CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid )
244	cnh	1.18	CcnhDebugEnds
245	cnh	1.1
246			RETURN
247			END