model/src/ini_cg2d.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.22 1998/12/08 19:44:28 adcroft Exp $

#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 "CG2D.h"
#include "OBCS.h"

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 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) + 
     &      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) + 
     &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        IF (openBoundaries) 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
        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
      cg2dNbuf(1,myThid) = myNorm
      _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )
      IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
       myNorm = 1. _d 0/cg2dNbuf(1,1)
      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
      CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
      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) = 0. _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
      CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
      CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
      CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 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 J=1,sNy
          DO I=1,sNx
           cg2d_x(I,J,bi,bj) = 0. _d 0
           cg2d_b(I,J,bi,bj) = 0. _d 0
#ifdef INCLUDE_CD_CODE
           cg2d_xNM1(I,J,bi,bj) = 0. _d 0
#endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
C--    Update overlap regions
       _EXCH_XY_R8(cg2d_x, myThid)
       _EXCH_XY_R8(cg2d_b, myThid)
#ifdef INCLUDE_CD_CODE
       _EXCH_XY_R8(cg2d_xNM1, myThid)
#endif
      ENDIF

      RETURN
      END
1	adcroft	1.23	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.22 1998/12/08 19:44:28 adcroft Exp $
2	cnh	1.1
3	adcroft	1.13	#include "CPP_OPTIONS.h"
4	cnh	1.1
5			CStartOfInterface
6			SUBROUTINE INI_CG2D( myThid )
7			C /==========================================================\
8			C \| SUBROUTINE INI_CG2D \|
9			C \| o Initialise 2d conjugate gradient solver operators. \|
10			C \|==========================================================\|
11			C \| These arrays are purely a function of the basin geom. \|
12			C \| We set then here once and them use then repeatedly. \|
13			C \==========================================================/
14	adcroft	1.23	IMPLICIT NONE
15	cnh	1.1
16			C === Global variables ===
17			#include "SIZE.h"
18			#include "EEPARAMS.h"
19			#include "PARAMS.h"
20			#include "GRID.h"
21	cnh	1.12	#include "DYNVARS.h"
22	cnh	1.1	#include "CG2D.h"
23	adcroft	1.22	#include "OBCS.h"
24	cnh	1.1
25			C === Routine arguments ===
26			C myThid - Thread no. that called this routine.
27			INTEGER myThid
28			CEndOfInterface
29
30			C === Local variables ===
31			C xG, yG - Global coordinate location.
32			C zG
33			C iG, jG - Global coordinate index
34			C bi,bj - Loop counters
35			C faceArea - Temporary used to hold cell face areas.
36			C I,J,K
37			C myNorm - Work variable used in clculating normalisation factor
38			CHARACTER*(MAX_LEN_MBUF) msgBuf
39			INTEGER bi, bj
40			INTEGER I, J, K
41	cnh	1.7	_RL faceArea
42	cnh	1.15	_RS myNorm
43	cnh	1.4	_RL aC, aCw, aCs
44	cnh	1.1
45			C-- Initialise laplace operator
46			C aW2d: integral in Z Ax/dX
47			C aS2d: integral in Z Ay/dY
48			myNorm = 0. _d 0
49			DO bj=myByLo(myThid),myByHi(myThid)
50			DO bi=myBxLo(myThid),myBxHi(myThid)
51			DO J=1,sNy
52			DO I=1,sNx
53			aW2d(I,J,bi,bj) = 0. _d 0
54			aS2d(I,J,bi,bj) = 0. _d 0
55			ENDDO
56			ENDDO
57	cnh	1.17	DO K=1,Nr
58	cnh	1.1	DO J=1,sNy
59			DO I=1,sNx
60	cnh	1.20	faceArea = _dyG(I,J,bi,bj)*drF(K)
61			& *_hFacW(I,J,K,bi,bj)
62	cnh	1.1	aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
63	cnh	1.17	& gBarofaceArearecip_dxC(I,J,bi,bj)
64	cnh	1.20	faceArea = _dxG(I,J,bi,bj)*drF(K)
65			& *_hFacS(I,J,K,bi,bj)
66	cnh	1.1	aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
67	cnh	1.17	& gBarofaceArearecip_dyC(I,J,bi,bj)
68	cnh	1.1	ENDDO
69			ENDDO
70			ENDDO
71	adcroft	1.22	IF (openBoundaries) THEN
72			DO I=1,sNx
73			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
74			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
75			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
76			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
77			ENDDO
78			DO J=1,sNy
79			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
80			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
81			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
82			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
83			ENDDO
84			ENDIF
85	cnh	1.1	DO J=1,sNy
86			DO I=1,sNx
87			myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
88			myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
89			ENDDO
90			ENDDO
91			ENDDO
92			ENDDO
93			cg2dNbuf(1,myThid) = myNorm
94	adcroft	1.16	_GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )
95	cnh	1.1	IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
96			myNorm = 1. _d 0/cg2dNbuf(1,1)
97			ELSE
98			myNorm = 1. _d 0
99			ENDIF
100	cnh	1.5	cg2dNorm = myNorm
101	cnh	1.1	_BEGIN_MASTER( myThid )
102			CcnhDebugStarts
103	cnh	1.20	WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
104			& cg2dNorm
105	cnh	1.3	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
106			WRITE(msgBuf,*) ' '
107	cnh	1.1	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
108			CcnhDebugEnds
109			_END_MASTER( myThid )
110			DO bj=myByLo(myThid),myByHi(myThid)
111			DO bi=myBxLo(myThid),myBxHi(myThid)
112			DO J=1,sNy
113			DO I=1,sNx
114			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
115			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
116			ENDDO
117			ENDDO
118			ENDDO
119			ENDDO
120
121			C-- Update overlap regions
122			CcnhDebugStarts
123	cnh	1.14	C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
124			C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
125	cnh	1.1	CcnhDebugEnds
126			_EXCH_XY_R4(aW2d, myThid)
127			_EXCH_XY_R4(aS2d, myThid)
128			CcnhDebugStarts
129	cnh	1.18	CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
130			CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
131	cnh	1.1	CcnhDebugEnds
132
133			C-- Initialise preconditioner
134	cnh	1.4	C Note. 20th May 1998
135			C I made a weird discovery! In the model paper we argue
136			C for the form of the preconditioner used here ( see
137			C A Finite-volume, Incompressible Navier-Stokes Model
138			C ...., Marshall et. al ). The algebra gives a simple
139			C 0.5 factor for the averaging of ac and aCw to get a
140			C symmettric pre-conditioner. By using a factor of 0.51
141			C i.e. scaling the off-diagonal terms in the
142			C preconditioner down slightly I managed to get the
143			C number of iterations for convergence in a test case to
144			C drop form 192 -> 134! Need to investigate this further!
145			C For now I have introduced a parameter cg2dpcOffDFac which
146			C defaults to 0.51 but can be set at runtime.
147	cnh	1.1	DO bj=myByLo(myThid),myByHi(myThid)
148			DO bi=myBxLo(myThid),myBxHi(myThid)
149			DO J=1,sNy
150			DO I=1,sNx
151			pC(I,J,bi,bj) = 1. _d 0
152	cnh	1.4	aC = -(
153			& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
154			& +aS2d(I,J,bi,bj) + aS2D(I ,J+1,bi,bj)
155	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
156	cnh	1.17	& rA(I,J,bi,bj)/deltaTMom/deltaTMom
157	cnh	1.4	& )
158			aCs = -(
159			& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
160			& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
161	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
162	cnh	1.17	& rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
163	cnh	1.4	& )
164			aCw = -(
165			& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
166			& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
167	cnh	1.19	& +freeSurfFacmyNorm horiVertRatio*
168	cnh	1.17	& rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
169	cnh	1.4	& )
170			IF ( aC .EQ. 0. ) THEN
171			pC(I,J,bi,bj) = 0. _d 0
172			ELSE
173			pC(I,J,bi,bj) = 1. _d 0 / aC
174			ENDIF
175			IF ( aC + aCw .EQ. 0. ) THEN
176			pW(I,J,bi,bj) = 0.
177			ELSE
178			pW(I,J,bi,bj) =
179			& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCw+aC))*2 )
180			ENDIF
181			IF ( aC + aCs .EQ. 0. ) THEN
182			pS(I,J,bi,bj) = 0.
183			ELSE
184			pS(I,J,bi,bj) =
185			& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCs+aC))*2 )
186			ENDIF
187	cnh	1.6	C pC(I,J,bi,bj) = 1.
188			C pW(I,J,bi,bj) = 0.
189			C pS(I,J,bi,bj) = 0.
190	cnh	1.1	ENDDO
191			ENDDO
192			ENDDO
193			ENDDO
194			C-- Update overlap regions
195			_EXCH_XY_R4(pC, myThid)
196			_EXCH_XY_R4(pW, myThid)
197			_EXCH_XY_R4(pS, myThid)
198	cnh	1.18	CcnhDebugStarts
199			CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid )
200			CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid )
201			CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid )
202			CcnhDebugEnds
203	cnh	1.1
204	cnh	1.5	C-- Set default values for initial guess and RHS
205	cnh	1.4	IF ( startTime .EQ. 0 ) THEN
206			DO bj=myByLo(myThid),myByHi(myThid)
207			DO bi=myBxLo(myThid),myBxHi(myThid)
208			DO J=1,sNy
209			DO I=1,sNx
210			cg2d_x(I,J,bi,bj) = 0. _d 0
211	cnh	1.5	cg2d_b(I,J,bi,bj) = 0. _d 0
212	cnh	1.21	#ifdef INCLUDE_CD_CODE
213	cnh	1.12	cg2d_xNM1(I,J,bi,bj) = 0. _d 0
214			#endif
215	cnh	1.4	ENDDO
216	cnh	1.1	ENDDO
217			ENDDO
218			ENDDO
219	cnh	1.4	C-- Update overlap regions
220			_EXCH_XY_R8(cg2d_x, myThid)
221	cnh	1.5	_EXCH_XY_R8(cg2d_b, myThid)
222	cnh	1.21	#ifdef INCLUDE_CD_CODE
223	cnh	1.12	_EXCH_XY_R8(cg2d_xNM1, myThid)
224			#endif
225	cnh	1.4	ENDIF
226	cnh	1.1
227			RETURN
228			END