/[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.21 - (hide annotations) (download)
Tue Oct 17 18:52:34 2006 UTC (17 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58q_post, checkpoint58r_post
Changes since 1.20: +2 -2 lines
clean-up multi-threaded problems (reported by debugger tcheck on ACES).

1 jmc 1.21 C $Header: /u/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.20 2006/02/23 20:55:49 jmc 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 jmc 1.18 C !ROUTINE: INI_CG3D
9 cnh 1.13 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.19 #include "SURFACE.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    
42 cnh 1.13 C !LOCAL VARIABLES:
43 adcroft 1.1 C === Local variables ===
44     C bi,bj - Loop counters
45 jmc 1.18 C I,J,K - Loop counters
46 adcroft 1.1 C faceArea - Temporary used to hold cell face areas.
47     C myNorm - Work variable used in clculating normalisation factor
48     CHARACTER*(MAX_LEN_MBUF) msgBuf
49     INTEGER bi, bj
50 jmc 1.20 INTEGER I, J, K, ks
51 adcroft 1.1 _RL faceArea
52     _RS myNorm
53     _RL theRecip_Dr
54     _RL aU, aL, aW, aE, aN, aS, aC
55 jmc 1.20 _RL tmpFac, nh_Fac, igwFac
56 cnh 1.13 CEOP
57 adcroft 1.1
58     CcnhDebugStarts
59 jmc 1.18 c _RL phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
60 adcroft 1.1 CcnhDebugEnds
61    
62 jmc 1.18 C-- Initialise to zero over the full range of indices
63     DO bj=myByLo(myThid),myByHi(myThid)
64     DO bi=myBxLo(myThid),myBxHi(myThid)
65     DO K=1,Nr
66     C- From common bloc CG3D_R:
67     DO J=1-OLy,sNy+OLy
68     DO I=1-OLx,sNx+OLx
69     aW3d(I,J,K,bi,bj) = 0.
70     aS3d(I,J,K,bi,bj) = 0.
71     aV3d(I,J,K,bi,bj) = 0.
72 jmc 1.19 aC3d(I,J,K,bi,bj) = 0.
73 jmc 1.18 zMC (I,J,K,bi,bj) = 0.
74     zML (I,J,K,bi,bj) = 0.
75     zMU (I,J,K,bi,bj) = 0.
76     ENDDO
77     ENDDO
78     C- From common bloc CG3D_WK_R:
79     DO J=0,sNy+1
80     DO I=0,sNx+1
81     cg3d_q(I,J,K,bi,bj) = 0.
82     cg3d_r(I,J,K,bi,bj) = 0.
83     cg3d_s(I,J,K,bi,bj) = 0.
84     ENDDO
85     ENDDO
86     C- From common bloc SFP3D_COMMON_R:
87     DO J=1-OLy,sNy+OLy
88     DO I=1-OLx,sNx+OLx
89     cg3d_b(I,J,K,bi,bj) = 0.
90     ENDDO
91     ENDDO
92     ENDDO
93     ENDDO
94     ENDDO
95    
96 jmc 1.20 nh_Fac = 0.
97     igwFac = 0.
98     IF ( nonHydrostatic
99     & .AND. nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2
100     IF ( implicitIntGravWave ) igwFac = 1. _d 0
101    
102     IF ( use3Dsolver ) THEN
103 adcroft 1.1 C-- Initialise laplace operator
104     C aW3d: Ax/dX
105     C aS3d: Ay/dY
106     C aV3d: Ar/dR
107     myNorm = 0. _d 0
108     DO bj=myByLo(myThid),myByHi(myThid)
109     DO bi=myBxLo(myThid),myBxHi(myThid)
110     DO K=1,Nr
111     DO J=1,sNy
112 jmc 1.19 DO I=1,sNx+1
113 adcroft 1.1 faceArea = _dyG(I,J,bi,bj)*drF(K)
114     & *_hFacW(I,J,K,bi,bj)
115 jmc 1.10 aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
116 jmc 1.19 myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
117     ENDDO
118     ENDDO
119     DO J=1,sNy+1
120     DO I=1,sNx
121 adcroft 1.1 faceArea = _dxG(I,J,bi,bj)*drF(K)
122     & *_hFacS(I,J,K,bi,bj)
123 jmc 1.10 aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
124 adcroft 1.1 myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
125     ENDDO
126     ENDDO
127     ENDDO
128     DO K=1,1
129     DO J=1,sNy
130     DO I=1,sNx
131     aV3d(I,J,K,bi,bj) = 0.
132     myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
133     ENDDO
134     ENDDO
135     ENDDO
136     DO K=2,Nr
137 jmc 1.20 tmpFac = nh_Fac*recip_horiVertRatio*recip_horiVertRatio
138     & + igwFac*dBdrRef(k)*deltaTMom*dTtracerLev(k)
139     IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
140 adcroft 1.1 DO J=1,sNy
141     DO I=1,sNx
142 jmc 1.18 faceArea = _rA(I,J,bi,bj)*maskC(I,J, K ,bi,bj)
143     & *maskC(I,J,K-1,bi,bj)
144     theRecip_Dr = recip_drC(K)
145     c theRecip_Dr =
146 adcroft 1.1 caja & drF(K )*_hFacC(i,j,k ,bi,bj)*0.5
147     caja & +drF(K-1)*_hFacC(i,j,k-1,bi,bj)*0.5
148 jmc 1.18 c IF ( theRecip_Dr .NE. 0. )
149     c & theRecip_Dr = 1. _d 0/theRecip_Dr
150 jmc 1.20 aV3d(I,J,K,bi,bj) = faceArea*theRecip_Dr*tmpFac
151 adcroft 1.1 myNorm = MAX(ABS(aV3d(I,J,K,bi,bj)),myNorm)
152     ENDDO
153     ENDDO
154 adcroft 1.5 ENDDO
155 adcroft 1.4 #ifdef ALLOW_OBCS
156 jmc 1.18 IF ( useOBCS ) THEN
157     DO K=1,Nr
158 adcroft 1.1 DO I=1,sNx
159     IF (OB_Jn(I,bi,bj).NE.0) THEN
160     aS3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
161     aS3d(I,OB_Jn(I,bi,bj)+1,K,bi,bj)=0.
162     aW3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
163     aW3d(I+1,OB_Jn(I,bi,bj),K,bi,bj)=0.
164 adcroft 1.7 aV3d(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
165 adcroft 1.1 ENDIF
166     IF (OB_Js(I,bi,bj).NE.0) THEN
167     aS3d(I,OB_Js(I,bi,bj)+1,K,bi,bj)=0.
168     aS3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
169     aW3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
170     aW3d(I+1,OB_Js(I,bi,bj),K,bi,bj)=0.
171 adcroft 1.7 aV3d(I,OB_Js(I,bi,bj),K,bi,bj)=0.
172 adcroft 1.1 ENDIF
173     ENDDO
174     DO J=1,sNy
175     IF (OB_Ie(J,bi,bj).NE.0) THEN
176     aW3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
177     aW3d(OB_Ie(J,bi,bj)+1,J,K,bi,bj)=0.
178     aS3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
179     aS3d(OB_Ie(J,bi,bj),J+1,K,bi,bj)=0.
180 adcroft 1.8 aV3d(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
181 adcroft 1.1 ENDIF
182     IF (OB_Iw(J,bi,bj).NE.0) THEN
183     aW3d(OB_Iw(J,bi,bj)+1,J,K,bi,bj)=0.
184     aW3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
185     aS3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
186     aS3d(OB_Iw(J,bi,bj),J+1,K,bi,bj)=0.
187 adcroft 1.8 aV3d(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
188 adcroft 1.1 ENDIF
189     ENDDO
190 jmc 1.18 ENDDO
191     ENDIF
192 adcroft 1.4 #endif
193 adcroft 1.1 ENDDO
194     ENDDO
195 adcroft 1.2 _GLOBAL_MAX_R4( myNorm, myThid )
196     IF ( myNorm .NE. 0. _d 0 ) THEN
197     myNorm = 1. _d 0/myNorm
198 adcroft 1.1 ELSE
199     myNorm = 1. _d 0
200     ENDIF
201 jmc 1.21 _BEGIN_MASTER( myThid )
202 adcroft 1.1 cg3dNorm = myNorm
203     CcnhDebugStarts
204     WRITE(msgBuf,'(A,E40.25)') '// CG3D normalisation factor = '
205     & , cg3dNorm
206     CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
207     WRITE(msgBuf,*) ' '
208     CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
209     CcnhDebugEnds
210     _END_MASTER( myThid )
211     DO bj=myByLo(myThid),myByHi(myThid)
212     DO bi=myBxLo(myThid),myBxHi(myThid)
213 jmc 1.19 C- Set solver main diagonal term
214     DO K=1,Nr
215     DO J=1,sNy
216     DO I=1,sNx
217     aW = aW3d(I ,J,K,bi,bj)
218     aE = aW3d(I+1,J,K,bi,bj)
219     aN = aS3d(I,J+1,K,bi,bj)
220     aS = aS3d(I,J ,K,bi,bj)
221     aU = aV3d(I,J,K,bi,bj)
222     IF ( K .NE. Nr ) THEN
223     aL = aV3d(I,J,K+1,bi,bj)
224     ELSE
225     aL = 0.
226     ENDIF
227     aC3d(I,J,K,bi,bj) = -aW-aE-aN-aS-aU-aL
228     ENDDO
229     ENDDO
230     ENDDO
231     C- Add free-surface source term
232     DO J=1,sNy
233     DO I=1,sNx
234     ks = ksurfC(I,J,bi,bj)
235     IF ( ks.LE.Nr ) THEN
236     aC3d(I,J,ks,bi,bj) = aC3d(I,J,ks,bi,bj)
237     & -freeSurfFac*recip_Bo(I,J,bi,bj)
238     & *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
239     ENDIF
240     ENDDO
241     ENDDO
242     C- Matrix solver normalisation
243 adcroft 1.1 DO K=1,Nr
244     DO J=1,sNy
245     DO I=1,sNx
246     aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
247     aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
248     aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
249 jmc 1.19 aC3d(I,J,K,bi,bj) = aC3d(I,J,K,bi,bj)*myNorm
250 adcroft 1.1 ENDDO
251     ENDDO
252     ENDDO
253     ENDDO
254     ENDDO
255 jmc 1.19
256 adcroft 1.1 C-- Update overlap regions
257 jmc 1.18 c _EXCH_XYZ_R4(aW3d, myThid)
258     c _EXCH_XYZ_R4(aS3d, myThid)
259     CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
260 adcroft 1.1 _EXCH_XYZ_R4(aV3d, myThid)
261 jmc 1.19 _EXCH_XYZ_R4(aC3d, myThid)
262 adcroft 1.1 CcnhDebugStarts
263     C CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
264     C CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
265     CcnhDebugEnds
266    
267     C-- Initialise preconditioner
268     C For now PC is just the identity. Change to
269     C be LU factorization of d2/dz2 later. Note
270     C check for consistency with S/R CG3D before
271     C assuming zML is lower and zMU is upper!
272     DO bj=myByLo(myThid),myByHi(myThid)
273     DO bi=myBxLo(myThid),myBxHi(myThid)
274 jmc 1.19 DO K=1,Nr
275 adcroft 1.1 DO J=1,sNy
276     DO I=1,sNx
277 jmc 1.19 IF ( aC3d(I,J,K,bi,bj) .NE. 0. ) THEN
278     zMC(i,j,k,bi,bj) = aC3d(I,J,K,bi,bj)
279     zML(i,j,k,bi,bj) = aV3d(I,J,K,bi,bj)
280     zMU(i,j,k,bi,bj) = 0.
281     IF ( K.NE.Nr )
282     & zMU(i,j,k,bi,bj) = aV3d(I,J,K+1,bi,bj)
283 adcroft 1.1 CcnhDebugStarts
284     C zMC(i,j,k,bi,bj) = 1.
285     C zMU(i,j,k,bi,bj) = 0.
286     C zML(i,j,k,bi,bj) = 0.
287     CcnhDebugEnds
288     ELSE
289     zMC(i,j,k,bi,bj) = 1. _d 0
290     zMU(i,j,k,bi,bj) = 0.
291     zML(i,j,k,bi,bj) = 0.
292     ENDIF
293     ENDDO
294     ENDDO
295     ENDDO
296     DO J=1,sNy
297     DO I=1,sNx
298     zMC(i,j,1,bi,bj)=
299     & 1. _d 0 / zMC(i,j,1,bi,bj)
300     zMU(i,j,1,bi,bj)=
301     & zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
302     ENDDO
303     ENDDO
304 jmc 1.19 DO K=2,Nr
305 adcroft 1.1 DO J=1,sNy
306     DO I=1,sNx
307     zMC(i,j,k,bi,bj) = 1. _d 0 /
308     & (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
309     zMU(i,j,k,bi,bj)=zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
310     ENDDO
311     ENDDO
312     ENDDO
313 jmc 1.19 DO K=1,Nr
314 adcroft 1.1 DO J=1,sNy
315     DO I=1,sNx
316     aW = aW3d(I ,J,K,bi,bj)
317     aE = aW3d(I+1,J,K,bi,bj)
318     aN = aS3d(I,J+1,K,bi,bj)
319     aS = aS3d(I,J ,K,bi,bj)
320     IF ( K .NE. 1 ) THEN
321     aU = aV3d(I,J,K-1,bi,bj)
322     ELSE
323 jmc 1.19 aU = 0.
324 adcroft 1.1 ENDIF
325     IF ( K .NE. Nr ) THEN
326     aL = aV3d(I,J,K+1,bi,bj)
327     ELSE
328 jmc 1.19 aL = 0.
329 adcroft 1.1 ENDIF
330     aC = -aW-aE-aN-aS-aU-aL
331     IF ( aC .EQ. 0. ) THEN
332 adcroft 1.5 zMC(i,j,k,bi,bj) = 1.
333 adcroft 1.1 zML(i,j,k,bi,bj) = 0.
334     zMU(i,j,k,bi,bj) = 0.
335     CcnhDebugStarts
336     C ELSE
337     C zMC(i,j,k,bi,bj) = 1.
338     C zML(i,j,k,bi,bj) = 0.
339     C zMU(i,j,k,bi,bj) = 0.
340     CcnhDEbugEnds
341     ENDIF
342     ENDDO
343     ENDDO
344     ENDDO
345     ENDDO
346     ENDDO
347     C-- Update overlap regions
348     _EXCH_XYZ_R4(zMC, myThid)
349     _EXCH_XYZ_R4(zML, myThid)
350     _EXCH_XYZ_R4(zMU, myThid)
351    
352     CcnhDebugStarts
353 jmc 1.18 c DO k=1,Nr
354     c DO j=1-OLy,sNy+OLy
355     c DO i=1-OLx,sNx+OLx
356     c phi(i,j,1,1) = zMc(i,j,k,1,1)
357     c ENDDO
358     c ENDDO
359 adcroft 1.1 C CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
360 jmc 1.18 c ENDDO
361 adcroft 1.1 C CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
362     C CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
363     CcnhDebugEnds
364    
365 jmc 1.20 C-- end if (use3Dsolver)
366     ENDIF
367    
368 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
369    
370     RETURN
371     END

  ViewVC Help
Powered by ViewVC 1.1.22