/[MITgcm]/MITgcm/model/src/ini_cg2d.F
ViewVC logotype

Annotation of /MITgcm/model/src/ini_cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.14 - (hide annotations) (download)
Mon Jun 15 05:13:56 1998 UTC (25 years, 11 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint7, checkpoint9, checkpoint8
Branch point for: checkpoint7-4degree-ref
Changes since 1.13: +5 -5 lines 
Fairly coplete 4 degree global intercomparison
setup.
 Includes changes to make convective adjustment and hydrostatic
pressure correct as well as IO for climatological datasets
1 cnh 1.14 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_cg2d.F,v 1.13 1998/06/09 15:58:36 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    
15     C === Global variables ===
16     #include "SIZE.h"
17     #include "EEPARAMS.h"
18     #include "PARAMS.h"
19     #include "GRID.h"
20 cnh 1.12 #include "DYNVARS.h"
21 cnh 1.1 #include "CG2D.h"
22    
23     C === Routine arguments ===
24     C myThid - Thread no. that called this routine.
25     INTEGER myThid
26     CEndOfInterface
27    
28     C === Local variables ===
29     C xG, yG - Global coordinate location.
30     C zG
31     C iG, jG - Global coordinate index
32     C bi,bj - Loop counters
33     C faceArea - Temporary used to hold cell face areas.
34     C I,J,K
35     C myNorm - Work variable used in clculating normalisation factor
36     CHARACTER*(MAX_LEN_MBUF) msgBuf
37     INTEGER bi, bj
38     INTEGER I, J, K
39 cnh 1.7 _RL faceArea
40 cnh 1.1 _RL myNorm
41 cnh 1.4 _RL aC, aCw, aCs
42 cnh 1.1
43     C-- Initialise laplace operator
44     C aW2d: integral in Z Ax/dX
45     C aS2d: integral in Z Ay/dY
46     myNorm = 0. _d 0
47     DO bj=myByLo(myThid),myByHi(myThid)
48     DO bi=myBxLo(myThid),myBxHi(myThid)
49     DO J=1,sNy
50     DO I=1,sNx
51     aW2d(I,J,bi,bj) = 0. _d 0
52     aS2d(I,J,bi,bj) = 0. _d 0
53     ENDDO
54     ENDDO
55     DO K=1,Nz
56     DO J=1,sNy
57     DO I=1,sNx
58 cnh 1.10 faceArea = _dyG(I,J,bi,bj)*dzF(K)*_hFacW(I,J,K,bi,bj)
59 cnh 1.1 aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
60 cnh 1.10 & gBaro*faceArea*_rdxC(I,J,bi,bj)
61     faceArea = _dxG(I,J,bi,bj)*dzF(K)*_hFacS(I,J,K,bi,bj)
62 cnh 1.1 aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
63 cnh 1.11 & gBaro*faceArea*_rdyC(I,J,bi,bj)
64 cnh 1.1 ENDDO
65     ENDDO
66     ENDDO
67     DO J=1,sNy
68     DO I=1,sNx
69     myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
70     myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
71     ENDDO
72     ENDDO
73     ENDDO
74     ENDDO
75     cg2dNbuf(1,myThid) = myNorm
76     _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )
77     IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
78     myNorm = 1. _d 0/cg2dNbuf(1,1)
79     ELSE
80     myNorm = 1. _d 0
81     ENDIF
82 cnh 1.5 cg2dNorm = myNorm
83 cnh 1.1 _BEGIN_MASTER( myThid )
84     CcnhDebugStarts
85 cnh 1.3 WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm
86     CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
87     WRITE(msgBuf,*) ' '
88 cnh 1.1 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
89     CcnhDebugEnds
90     _END_MASTER( myThid )
91     DO bj=myByLo(myThid),myByHi(myThid)
92     DO bi=myBxLo(myThid),myBxHi(myThid)
93     DO J=1,sNy
94     DO I=1,sNx
95     aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
96     aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
97     ENDDO
98     ENDDO
99     ENDDO
100     ENDDO
101    
102     C-- Update overlap regions
103     CcnhDebugStarts
104 cnh 1.14 C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
105     C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
106 cnh 1.1 CcnhDebugEnds
107     _EXCH_XY_R4(aW2d, myThid)
108     _EXCH_XY_R4(aS2d, myThid)
109     CcnhDebugStarts
110 cnh 1.14 C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
111     C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
112 cnh 1.1 CcnhDebugEnds
113    
114     C-- Initialise preconditioner
115 cnh 1.4 C Note. 20th May 1998
116     C I made a weird discovery! In the model paper we argue
117     C for the form of the preconditioner used here ( see
118     C A Finite-volume, Incompressible Navier-Stokes Model
119     C ...., Marshall et. al ). The algebra gives a simple
120     C 0.5 factor for the averaging of ac and aCw to get a
121     C symmettric pre-conditioner. By using a factor of 0.51
122     C i.e. scaling the off-diagonal terms in the
123     C preconditioner down slightly I managed to get the
124     C number of iterations for convergence in a test case to
125     C drop form 192 -> 134! Need to investigate this further!
126     C For now I have introduced a parameter cg2dpcOffDFac which
127     C defaults to 0.51 but can be set at runtime.
128 cnh 1.1 DO bj=myByLo(myThid),myByHi(myThid)
129     DO bi=myBxLo(myThid),myBxHi(myThid)
130     DO J=1,sNy
131     DO I=1,sNx
132     pC(I,J,bi,bj) = 1. _d 0
133 cnh 1.4 aC = -(
134     & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
135     & +aS2d(I,J,bi,bj) + aS2D(I ,J+1,bi,bj)
136 cnh 1.5 & +freeSurfFac*myNorm*
137     & zA(I,J,bi,bj)/deltaTMom/deltaTMom
138 cnh 1.4 & )
139     aCs = -(
140     & aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
141     & +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
142 cnh 1.5 & +freeSurfFac*myNorm*
143     & zA(I,J-1,bi,bj)/deltaTMom/deltaTMom
144 cnh 1.4 & )
145     aCw = -(
146     & aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
147     & +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
148 cnh 1.5 & +freeSurfFac*myNorm*
149     & zA(I-1,J,bi,bj)/deltaTMom/deltaTMom
150 cnh 1.4 & )
151     IF ( aC .EQ. 0. ) THEN
152     pC(I,J,bi,bj) = 0. _d 0
153     ELSE
154     pC(I,J,bi,bj) = 1. _d 0 / aC
155     ENDIF
156     IF ( aC + aCw .EQ. 0. ) THEN
157     pW(I,J,bi,bj) = 0.
158     ELSE
159     pW(I,J,bi,bj) =
160     & -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
161     ENDIF
162     IF ( aC + aCs .EQ. 0. ) THEN
163     pS(I,J,bi,bj) = 0.
164     ELSE
165     pS(I,J,bi,bj) =
166     & -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
167     ENDIF
168 cnh 1.6 C pC(I,J,bi,bj) = 1.
169     C pW(I,J,bi,bj) = 0.
170     C pS(I,J,bi,bj) = 0.
171 cnh 1.1 ENDDO
172     ENDDO
173     ENDDO
174     ENDDO
175     C-- Update overlap regions
176     _EXCH_XY_R4(pC, myThid)
177     _EXCH_XY_R4(pW, myThid)
178     _EXCH_XY_R4(pS, myThid)
179    
180 cnh 1.5 C-- Set default values for initial guess and RHS
181 cnh 1.4 IF ( startTime .EQ. 0 ) THEN
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     cg2d_x(I,J,bi,bj) = 0. _d 0
187 cnh 1.5 cg2d_b(I,J,bi,bj) = 0. _d 0
188 cnh 1.12 #ifdef ALLOW_CD
189     cg2d_xNM1(I,J,bi,bj) = 0. _d 0
190     #endif
191 cnh 1.4 ENDDO
192 cnh 1.1 ENDDO
193     ENDDO
194     ENDDO
195 cnh 1.4 C-- Update overlap regions
196     _EXCH_XY_R8(cg2d_x, myThid)
197 cnh 1.5 _EXCH_XY_R8(cg2d_b, myThid)
198 cnh 1.12 #ifdef ALLOW_CD
199     _EXCH_XY_R8(cg2d_xNM1, myThid)
200     #endif
201 cnh 1.4 ENDIF
202 cnh 1.1
203     RETURN
204     END
  ViewVC Help
Powered by ViewVC 1.1.22