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

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

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

Revision 1.15 - (hide annotations) (download)
Fri Oct 31 20:35:32 2003 UTC (20 years, 8 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 edhill 1.15 C $Header: /u/u3/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.14 2003/10/09 04:19:18 edhill Exp $
2 jmc 1.10 C $Name: $
3 adcroft 1.1
4 edhill 1.15 #include "PACKAGES_CONFIG.h"
5 adcroft 1.1 #include "CPP_OPTIONS.h"
6    
7 cnh 1.13 CBOP
8     C !ROUTINE: INI_CG2D
9     C !INTERFACE:
10 adcroft 1.1 SUBROUTINE INI_CG3D( myThid )
11 cnh 1.13 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 adcroft 1.3 IMPLICIT NONE
23 adcroft 1.1 C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28 jmc 1.10 c_jmc #include "DYNVARS.h"
29 adcroft 1.1 #include "CG3D.h"
30 adcroft 1.12 #include "SOLVE_FOR_PRESSURE3D.h"
31 adcroft 1.4 #ifdef ALLOW_OBCS
32 adcroft 1.1 #include "OBCS.h"
33 adcroft 1.4 #endif
34 adcroft 1.1
35 cnh 1.13 C !INPUT/OUTPUT PARAMETERS:
36 adcroft 1.1 C === Routine arguments ===
37     C myThid - Thread no. that called this routine.
38     INTEGER myThid
39    
40 adcroft 1.6 #ifdef ALLOW_NONHYDROSTATIC
41 edhill 1.14 ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
42 adcroft 1.6
43 cnh 1.13 C !LOCAL VARIABLES:
44 adcroft 1.1 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 cnh 1.13 CEOP
61 adcroft 1.1
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 jmc 1.10 aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
79 adcroft 1.1 faceArea = _dxG(I,J,bi,bj)*drF(K)
80     & *_hFacS(I,J,K,bi,bj)
81 jmc 1.10 aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
82 adcroft 1.1 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 jmc 1.10 aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr
110     & *horiVertRatio*horiVertRatio
111 adcroft 1.1 myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
112     ENDDO
113     ENDDO
114 adcroft 1.5 ENDDO
115 adcroft 1.4 #ifdef ALLOW_OBCS
116 edhill 1.14 ceh3 needs an IF ( useOBCS ) THEN
117 adcroft 1.5 DO K=1,Nr
118 adcroft 1.8 IF (useOBCS) THEN
119 adcroft 1.1 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 adcroft 1.7 aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
126 adcroft 1.1 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 adcroft 1.7 aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
133 adcroft 1.1 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 adcroft 1.8 aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
142 adcroft 1.1 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 adcroft 1.8 aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
149 adcroft 1.1 ENDIF
150     ENDDO
151     ENDIF
152 adcroft 1.5 ENDDO
153 adcroft 1.4 #endif
154 adcroft 1.1 ENDDO
155     ENDDO
156 adcroft 1.2 _GLOBAL_MAX_R4( myNorm, myThid )
157     IF ( myNorm .NE. 0. _d 0 ) THEN
158     myNorm = 1. _d 0/myNorm
159 adcroft 1.1 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 jmc 1.10 & -freeSurfFac*myNorm*(horiVertRatio/gravity)*
223     & rA(I,J,bi,bj)/deltaTMom/deltaTMom
224 adcroft 1.1 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 adcroft 1.5 zMC(i,j,k,bi,bj) = 1.
279 adcroft 1.1 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