model/src/ini_cg2d.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.3 1998/04/26 23:41:54 cnh Exp $

#include "CPP_EEOPTIONS.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     \==========================================================/

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "CG2D.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
      real    faceArea
      _RL     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,Nz
         DO J=1,sNy
          DO I=1,sNx
           faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)
           aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) + 
     &      gravity*faceArea*rDxC(I,J,bi,bj)
           faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)
           aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) + 
     &      gravity*faceArea*rDyC(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        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_R8( 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
      _BEGIN_MASTER( myThid )
       cg2dNorm = myNorm
CcnhDebugStarts
       WRITE(msgBuf,'(A,F)') '// 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_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
CcnhDebugEnds
      _EXCH_XY_R4(aW2d, myThid)
      _EXCH_XY_R4(aS2d, myThid)
CcnhDebugStarts
C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
C     CALL PLOT_FIELD_XYR8( 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)
     &    )
          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)
     &    )
          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)
     &    )
          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
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C--   Update overlap regions
      _EXCH_XY_R4(pC, myThid)
      _EXCH_XY_R4(pW, myThid)
      _EXCH_XY_R4(pS, myThid)

C--   Set default values for initial guess
      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
          ENDDO
         ENDDO
        ENDDO
       ENDDO
C--    Update overlap regions
       _EXCH_XY_R8(cg2d_x, myThid)
      ENDIF

      RETURN
      END
1	cnh	1.4	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.3 1998/04/26 23:41:54 cnh Exp $
2	cnh	1.1
3			#include "CPP_EEOPTIONS.h"
4
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
15			C === Global variables ===
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "GRID.h"
20			#include "CG2D.h"
21
22			C === Routine arguments ===
23			C myThid - Thread no. that called this routine.
24			INTEGER myThid
25			CEndOfInterface
26
27			C === Local variables ===
28			C xG, yG - Global coordinate location.
29			C zG
30			C iG, jG - Global coordinate index
31			C bi,bj - Loop counters
32			C faceArea - Temporary used to hold cell face areas.
33			C I,J,K
34			C myNorm - Work variable used in clculating normalisation factor
35			CHARACTER*(MAX_LEN_MBUF) msgBuf
36			INTEGER bi, bj
37			INTEGER I, J, K
38			real faceArea
39			_RL myNorm
40	cnh	1.4	_RL aC, aCw, aCs
41	cnh	1.1
42			C-- Initialise laplace operator
43			C aW2d: integral in Z Ax/dX
44			C aS2d: integral in Z Ay/dY
45			myNorm = 0. _d 0
46			DO bj=myByLo(myThid),myByHi(myThid)
47			DO bi=myBxLo(myThid),myBxHi(myThid)
48			DO J=1,sNy
49			DO I=1,sNx
50			aW2d(I,J,bi,bj) = 0. _d 0
51			aS2d(I,J,bi,bj) = 0. _d 0
52			ENDDO
53			ENDDO
54			DO K=1,Nz
55			DO J=1,sNy
56			DO I=1,sNx
57			faceArea = dyG(I,J,bi,bj)dzF(K)HFacW(I,J,K,bi,bj)
58			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
59			& gravityfaceArearDxC(I,J,bi,bj)
60			faceArea = dxG(I,J,bi,bj)dzF(K)HFacS(I,J,K,bi,bj)
61			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
62			& gravityfaceArearDyC(I,J,bi,bj)
63			ENDDO
64			ENDDO
65			ENDDO
66			DO J=1,sNy
67			DO I=1,sNx
68			myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
69			myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
70			ENDDO
71			ENDDO
72			ENDDO
73			ENDDO
74			cg2dNbuf(1,myThid) = myNorm
75			_GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )
76			IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
77			myNorm = 1. _d 0/cg2dNbuf(1,1)
78			ELSE
79			myNorm = 1. _d 0
80			ENDIF
81			_BEGIN_MASTER( myThid )
82			cg2dNorm = myNorm
83			CcnhDebugStarts
84	cnh	1.3	WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm
85			CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
86			WRITE(msgBuf,*) ' '
87	cnh	1.1	CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
88			CcnhDebugEnds
89			_END_MASTER( myThid )
90			DO bj=myByLo(myThid),myByHi(myThid)
91			DO bi=myBxLo(myThid),myBxHi(myThid)
92			DO J=1,sNy
93			DO I=1,sNx
94			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
95			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
96			ENDDO
97			ENDDO
98			ENDDO
99			ENDDO
100
101			C-- Update overlap regions
102			CcnhDebugStarts
103			C CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
104			C CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
105			CcnhDebugEnds
106			_EXCH_XY_R4(aW2d, myThid)
107			_EXCH_XY_R4(aS2d, myThid)
108			CcnhDebugStarts
109			C CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
110			C CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
111			CcnhDebugEnds
112
113			C-- Initialise preconditioner
114	cnh	1.4	C Note. 20th May 1998
115			C I made a weird discovery! In the model paper we argue
116			C for the form of the preconditioner used here ( see
117			C A Finite-volume, Incompressible Navier-Stokes Model
118			C ...., Marshall et. al ). The algebra gives a simple
119			C 0.5 factor for the averaging of ac and aCw to get a
120			C symmettric pre-conditioner. By using a factor of 0.51
121			C i.e. scaling the off-diagonal terms in the
122			C preconditioner down slightly I managed to get the
123			C number of iterations for convergence in a test case to
124			C drop form 192 -> 134! Need to investigate this further!
125			C For now I have introduced a parameter cg2dpcOffDFac which
126			C defaults to 0.51 but can be set at runtime.
127	cnh	1.1	DO bj=myByLo(myThid),myByHi(myThid)
128			DO bi=myBxLo(myThid),myBxHi(myThid)
129			DO J=1,sNy
130			DO I=1,sNx
131			pC(I,J,bi,bj) = 1. _d 0
132	cnh	1.4	aC = -(
133			& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
134			& +aS2d(I,J,bi,bj) + aS2D(I ,J+1,bi,bj)
135			& )
136			aCs = -(
137			& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
138			& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
139			& )
140			aCw = -(
141			& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
142			& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
143			& )
144			IF ( aC .EQ. 0. ) THEN
145			pC(I,J,bi,bj) = 0. _d 0
146			ELSE
147			pC(I,J,bi,bj) = 1. _d 0 / aC
148			ENDIF
149			IF ( aC + aCw .EQ. 0. ) THEN
150			pW(I,J,bi,bj) = 0.
151			ELSE
152			pW(I,J,bi,bj) =
153			& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCw+aC))*2 )
154			ENDIF
155			IF ( aC + aCs .EQ. 0. ) THEN
156			pS(I,J,bi,bj) = 0.
157			ELSE
158			pS(I,J,bi,bj) =
159			& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac (aCs+aC))*2 )
160			ENDIF
161	cnh	1.1	ENDDO
162			ENDDO
163			ENDDO
164			ENDDO
165			C-- Update overlap regions
166			_EXCH_XY_R4(pC, myThid)
167			_EXCH_XY_R4(pW, myThid)
168			_EXCH_XY_R4(pS, myThid)
169
170			C-- Set default values for initial guess
171	cnh	1.4	IF ( startTime .EQ. 0 ) THEN
172			DO bj=myByLo(myThid),myByHi(myThid)
173			DO bi=myBxLo(myThid),myBxHi(myThid)
174			DO J=1,sNy
175			DO I=1,sNx
176			cg2d_x(I,J,bi,bj) = 0. _d 0
177			ENDDO
178	cnh	1.1	ENDDO
179			ENDDO
180			ENDDO
181	cnh	1.4	C-- Update overlap regions
182			_EXCH_XY_R8(cg2d_x, myThid)
183			ENDIF
184	cnh	1.1
185			RETURN
186			END