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

Contents of /MITgcm/model/src/ini_cg3d.F

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

Revision 1.15 - (show annotations) (download)
Fri Oct 31 20:35:32 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint52l_pre, hrcube4, hrcube5, checkpoint57b_post, checkpoint52d_pre, checkpoint56b_post, checkpoint52j_pre, checkpoint54d_post, checkpoint54e_post, checkpoint57d_post, checkpoint52l_post, checkpoint52k_post, checkpoint55, checkpoint54, checkpoint57, checkpoint56, checkpoint53, checkpoint52, checkpoint52f_post, checkpoint54f_post, checkpoint51t_post, checkpoint55i_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55c_post, checkpoint52e_pre, checkpoint52e_post, checkpoint53d_post, checkpoint57a_post, checkpoint52b_pre, checkpoint54b_post, checkpoint52m_post, checkpoint55g_post, checkpoint52b_post, checkpoint52c_post, checkpoint57c_post, checkpoint52f_pre, checkpoint55d_post, checkpoint54a_pre, checkpoint53c_post, checkpoint55d_pre, checkpoint57c_pre, checkpoint55j_post, checkpoint54a_post, checkpoint55h_post, checkpoint51r_post, checkpoint57e_post, checkpoint55b_post, checkpoint53a_post, checkpoint55f_post, checkpoint52d_post, checkpoint53g_post, eckpoint57e_pre, checkpoint52a_pre, checkpoint52i_post, checkpoint52h_pre, checkpoint56a_post, checkpoint53f_post, checkpoint52j_post, checkpoint57f_pre, branch-netcdf, checkpoint52n_post, checkpoint53b_pre, checkpoint56c_post, checkpoint57a_pre, checkpoint55a_post, checkpoint53b_post, checkpoint52a_post, ecco_c52_e35, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint51u_post
Branch point for: branch-nonh, netcdf-sm0
Changes since 1.14: +2 -1 lines 
 o remove all '#include "PACACKAGES_CONFIG.h"' from model/inc/* and cleanup
   the verification tests that this breaks
 o this was confirmed to work for the basic tests ("testreport -ieee") on
   shelley
1 C $Header: /u/u3/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.14 2003/10/09 04:19:18 edhill Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 CBOP
8 C !ROUTINE: INI_CG2D
9 C !INTERFACE:
10 SUBROUTINE INI_CG3D( myThid )
11 C !DESCRIPTION: \bv
12 C *==========================================================*
13 C | SUBROUTINE INI_CG3D
14 C | o Initialise 3d 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 IMPLICIT NONE
23 C === Global variables ===
24 #include "SIZE.h"
25 #include "EEPARAMS.h"
26 #include "PARAMS.h"
27 #include "GRID.h"
28 c_jmc #include "DYNVARS.h"
29 #include "CG3D.h"
30 #include "SOLVE_FOR_PRESSURE3D.h"
31 #ifdef ALLOW_OBCS
32 #include "OBCS.h"
33 #endif
34
35 C !INPUT/OUTPUT PARAMETERS:
36 C === Routine arguments ===
37 C myThid - Thread no. that called this routine.
38 INTEGER myThid
39
40 #ifdef ALLOW_NONHYDROSTATIC
41 ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
42
43 C !LOCAL VARIABLES:
44 C === Local variables ===
45 C xG, yG - Global coordinate location.
46 C zG
47 C iG, jG - Global coordinate index
48 C bi,bj - Loop counters
49 C faceArea - Temporary used to hold cell face areas.
50 C I,J,K
51 C myNorm - Work variable used in clculating normalisation factor
52 CHARACTER*(MAX_LEN_MBUF) msgBuf
53 INTEGER bi, bj
54 INTEGER I, J, K
55 _RL faceArea
56 _RS myNorm
57 _RL aCw, aCs
58 _RL theRecip_Dr
59 _RL aU, aL, aW, aE, aN, aS, aC
60 CEOP
61
62 CcnhDebugStarts
63 _RL phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
64 CcnhDebugEnds
65
66 C-- Initialise laplace operator
67 C aW3d: Ax/dX
68 C aS3d: Ay/dY
69 C aV3d: Ar/dR
70 myNorm = 0. _d 0
71 DO bj=myByLo(myThid),myByHi(myThid)
72 DO bi=myBxLo(myThid),myBxHi(myThid)
73 DO K=1,Nr
74 DO J=1,sNy
75 DO I=1,sNx
76 faceArea = _dyG(I,J,bi,bj)*drF(K)
77 & *_hFacW(I,J,K,bi,bj)
78 aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
79 faceArea = _dxG(I,J,bi,bj)*drF(K)
80 & *_hFacS(I,J,K,bi,bj)
81 aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
82 myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
83 myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
84 ENDDO
85 ENDDO
86 ENDDO
87 DO K=1,1
88 DO J=1,sNy
89 DO I=1,sNx
90 aV3d(I,J,K,bi,bj) = 0.
91 myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
92 ENDDO
93 ENDDO
94 ENDDO
95 DO K=2,Nr
96 DO J=1,sNy
97 DO I=1,sNx
98 IF ( _hFacC(i,j,k,bi,bj) .EQ. 0. ) THEN
99 faceArea = 0.
100 ELSE
101 faceArea = _rA(I,J,bi,bj)
102 ENDIF
103 theRecip_Dr =
104 & drC(K)
105 caja & drF(K )*_hFacC(i,j,k ,bi,bj)*0.5
106 caja & +drF(K-1)*_hFacC(i,j,k-1,bi,bj)*0.5
107 IF ( theRecip_Dr .NE. 0. )
108 & theRecip_Dr = 1. _d 0/theRecip_Dr
109 aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr
110 & *horiVertRatio*horiVertRatio
111 myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
112 ENDDO
113 ENDDO
114 ENDDO
115 #ifdef ALLOW_OBCS
116 ceh3 needs an IF ( useOBCS ) THEN
117 DO K=1,Nr
118 IF (useOBCS) THEN
119 DO I=1,sNx
120 IF (OB_Jn(I,bi,bj).NE.0) THEN
121 aS3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
122 aS3d(I,OB_Jn(I,bi,bj)+1,K,bi,bj)=0.
123 aW3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
124 aW3d(I+1,OB_Jn(I,bi,bj),K,bi,bj)=0.
125 aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
126 ENDIF
127 IF (OB_Js(I,bi,bj).NE.0) THEN
128 aS3d(I,OB_Js(I,bi,bj)+1,K,bi,bj)=0.
129 aS3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
130 aW3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
131 aW3d(I+1,OB_Js(I,bi,bj),K,bi,bj)=0.
132 aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
133 ENDIF
134 ENDDO
135 DO J=1,sNy
136 IF (OB_Ie(J,bi,bj).NE.0) THEN
137 aW3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
138 aW3d(OB_Ie(J,bi,bj)+1,J,K,bi,bj)=0.
139 aS3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
140 aS3d(OB_Ie(J,bi,bj),J+1,K,bi,bj)=0.
141 aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
142 ENDIF
143 IF (OB_Iw(J,bi,bj).NE.0) THEN
144 aW3d(OB_Iw(J,bi,bj)+1,J,K,bi,bj)=0.
145 aW3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
146 aS3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
147 aS3d(OB_Iw(J,bi,bj),J+1,K,bi,bj)=0.
148 aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
149 ENDIF
150 ENDDO
151 ENDIF
152 ENDDO
153 #endif
154 ENDDO
155 ENDDO
156 _GLOBAL_MAX_R4( myNorm, myThid )
157 IF ( myNorm .NE. 0. _d 0 ) THEN
158 myNorm = 1. _d 0/myNorm
159 ELSE
160 myNorm = 1. _d 0
161 ENDIF
162 cg3dNorm = myNorm
163 _BEGIN_MASTER( myThid )
164 CcnhDebugStarts
165 WRITE(msgBuf,'(A,E40.25)') '// CG3D normalisation factor = '
166 & , cg3dNorm
167 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
168 WRITE(msgBuf,*) ' '
169 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
170 CcnhDebugEnds
171 _END_MASTER( myThid )
172 DO bj=myByLo(myThid),myByHi(myThid)
173 DO bi=myBxLo(myThid),myBxHi(myThid)
174 DO K=1,Nr
175 DO J=1,sNy
176 DO I=1,sNx
177 aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
178 aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
179 aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
180 ENDDO
181 ENDDO
182 ENDDO
183 ENDDO
184 ENDDO
185
186 C-- Update overlap regions
187 _EXCH_XYZ_R4(aW3d, myThid)
188 _EXCH_XYZ_R4(aS3d, myThid)
189 _EXCH_XYZ_R4(aV3d, myThid)
190 CcnhDebugStarts
191 C CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
192 C CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
193 CcnhDebugEnds
194
195 C-- Initialise preconditioner
196 C For now PC is just the identity. Change to
197 C be LU factorization of d2/dz2 later. Note
198 C check for consistency with S/R CG3D before
199 C assuming zML is lower and zMU is upper!
200 DO bj=myByLo(myThid),myByHi(myThid)
201 DO bi=myBxLo(myThid),myBxHi(myThid)
202 DO K=1,Nr
203 DO J=1,sNy
204 DO I=1,sNx
205 aW = aW3d(I ,J,K,bi,bj)
206 aE = aW3d(I+1,J,K,bi,bj)
207 aN = aS3d(I,J+1,K,bi,bj)
208 aS = aS3d(I,J ,K,bi,bj)
209 IF ( K .NE. 1 ) THEN
210 aU = aV3d(I,J,K,bi,bj)
211 ELSE
212 aU = 0
213 ENDIF
214 IF ( K .NE. Nr ) THEN
215 aL = aV3d(I,J,K+1,bi,bj)
216 ELSE
217 aL = 0
218 ENDIF
219 aC = -aW-aE-aN-aS-aU-aL
220 IF ( K .EQ. 1 .AND. aC .NE. 0. ) THEN
221 aC = aC
222 & -freeSurfFac*myNorm*(horiVertRatio/gravity)*
223 & rA(I,J,bi,bj)/deltaTMom/deltaTMom
224 ENDIF
225 IF ( aC .NE. 0. ) THEN
226 zMC(i,j,k,bi,bj) = aC
227 zMU(i,j,k,bi,bj) = aL
228 zML(i,j,k,bi,bj) = aU
229 CcnhDebugStarts
230 C zMC(i,j,k,bi,bj) = 1.
231 C zMU(i,j,k,bi,bj) = 0.
232 C zML(i,j,k,bi,bj) = 0.
233 CcnhDebugEnds
234 ELSE
235 zMC(i,j,k,bi,bj) = 1. _d 0
236 zMU(i,j,k,bi,bj) = 0.
237 zML(i,j,k,bi,bj) = 0.
238 ENDIF
239 ENDDO
240 ENDDO
241 ENDDO
242 DO J=1,sNy
243 DO I=1,sNx
244 zMC(i,j,1,bi,bj)=
245 & 1. _d 0 / zMC(i,j,1,bi,bj)
246 zMU(i,j,1,bi,bj)=
247 & zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
248 ENDDO
249 ENDDO
250 DO K=2,Nr
251 DO J=1,sNy
252 DO I=1,sNx
253 zMC(i,j,k,bi,bj) = 1. _d 0 /
254 & (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
255 zMU(i,j,k,bi,bj)=zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
256 ENDDO
257 ENDDO
258 ENDDO
259 DO K=1,Nr
260 DO J=1,sNy
261 DO I=1,sNx
262 aW = aW3d(I ,J,K,bi,bj)
263 aE = aW3d(I+1,J,K,bi,bj)
264 aN = aS3d(I,J+1,K,bi,bj)
265 aS = aS3d(I,J ,K,bi,bj)
266 IF ( K .NE. 1 ) THEN
267 aU = aV3d(I,J,K-1,bi,bj)
268 ELSE
269 aU = 0
270 ENDIF
271 IF ( K .NE. Nr ) THEN
272 aL = aV3d(I,J,K+1,bi,bj)
273 ELSE
274 aL = 0
275 ENDIF
276 aC = -aW-aE-aN-aS-aU-aL
277 IF ( aC .EQ. 0. ) THEN
278 zMC(i,j,k,bi,bj) = 1.
279 zML(i,j,k,bi,bj) = 0.
280 zMU(i,j,k,bi,bj) = 0.
281 CcnhDebugStarts
282 C ELSE
283 C zMC(i,j,k,bi,bj) = 1.
284 C zML(i,j,k,bi,bj) = 0.
285 C zMU(i,j,k,bi,bj) = 0.
286 CcnhDEbugEnds
287 ENDIF
288 ENDDO
289 ENDDO
290 ENDDO
291 ENDDO
292 ENDDO
293 C-- Update overlap regions
294 _EXCH_XYZ_R4(zMC, myThid)
295 _EXCH_XYZ_R4(zML, myThid)
296 _EXCH_XYZ_R4(zMU, myThid)
297
298 CcnhDebugStarts
299 DO k=1,nr
300 DO j=1-OLy,sNy+OLy
301 DO i=1-OLx,sNx+OLx
302 phi(i,j,1,1) = zMc(i,j,k,1,1)
303 ENDDO
304 ENDDO
305 C CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
306 ENDDO
307 C CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
308 C CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
309 CcnhDebugEnds
310
311 C-- Set default values for initial guess and RHS
312 IF ( startTime .EQ. 0 ) THEN
313 DO bj=myByLo(myThid),myByHi(myThid)
314 DO bi=myBxLo(myThid),myBxHi(myThid)
315 DO K=1,Nr
316 DO J=1-OLy,sNy+OLy
317 DO I=1-OLx,sNx+OLx
318 cg3d_b(I,J,K,bi,bj) = 0. _d 0
319 ENDDO
320 ENDDO
321 ENDDO
322 ENDDO
323 ENDDO
324 C-- Update overlap regions
325 _EXCH_XYZ_R8(cg3d_b, myThid)
326 ENDIF
327
328 #endif /* ALLOW_NONHYDROSTATIC */
329
330 RETURN
331 END
  ViewVC Help
Powered by ViewVC 1.1.22