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

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

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


Revision 1.77 - (show annotations) (download)
Wed Jun 8 01:21:14 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint63, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.76: +7 -7 lines
refine debugLevel criteria when printing messages

1 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.76 2011/05/18 23:41:26 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 CBOP
8 C !ROUTINE: SOLVE_FOR_PRESSURE
9 C !INTERFACE:
10 SUBROUTINE SOLVE_FOR_PRESSURE( myTime, myIter, myThid )
11
12 C !DESCRIPTION: \bv
13 C *==========================================================*
14 C | SUBROUTINE SOLVE_FOR_PRESSURE
15 C | o Controls inversion of two and/or three-dimensional
16 C | elliptic problems for the pressure field.
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 IMPLICIT NONE
22 C == Global variables
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "SURFACE.h"
28 #include "FFIELDS.h"
29 #include "DYNVARS.h"
30 #include "SOLVE_FOR_PRESSURE.h"
31 #ifdef ALLOW_NONHYDROSTATIC
32 #include "SOLVE_FOR_PRESSURE3D.h"
33 #include "NH_VARS.h"
34 #endif
35 #ifdef ALLOW_CD_CODE
36 #include "CD_CODE_VARS.h"
37 #endif
38
39 C === Functions ====
40 LOGICAL DIFFERENT_MULTIPLE
41 EXTERNAL DIFFERENT_MULTIPLE
42
43 C !INPUT/OUTPUT PARAMETERS:
44 C == Routine arguments ==
45 C myTime :: Current time in simulation
46 C myIter :: Current iteration number in simulation
47 C myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
48 _RL myTime
49 INTEGER myIter
50 INTEGER myThid
51
52 C !LOCAL VARIABLES:
53 C == Local variables ==
54 INTEGER i,j,k,bi,bj
55 INTEGER ks
56 INTEGER numIters
57 _RL firstResidual,lastResidual
58 _RL tmpFac
59 _RL sumEmP, tileEmP(nSx,nSy)
60 LOGICAL putPmEinXvector
61 INTEGER ioUnit
62 CHARACTER*10 sufx
63 CHARACTER*(MAX_LEN_MBUF) msgBuf
64 #ifdef ALLOW_NONHYDROSTATIC
65 LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
66 #else
67 _RL cg3d_b(1)
68 #endif
69 CEOP
70
71 #ifdef ALLOW_NONHYDROSTATIC
72 zeroPsNH = .FALSE.
73 c zeroPsNH = use3Dsolver .AND. exactConserv
74 c & .AND. select_rStar.EQ.0
75 zeroMeanPnh = .FALSE.
76 c zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
77 c oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
78 c & .AND. .NOT.zeroPsNH
79 oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
80 #else
81 cg3d_b(1) = 0.
82 #endif
83
84 C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
85 C anelastic (always Z-coordinate):
86 C 1) assume that rhoFacF(1)=1 (and ksurf == 1);
87 C (this reduces the number of lines of code to modify)
88 C 2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
89 C (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
90 C => 2 factors cancel in elliptic eq. for Phi_s ,
91 C but 1rst factor(a) remains in RHS cg2d_b.
92
93 C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
94 C instead of simply etaN ; This can speed-up the solver convergence in
95 C the case where |Global_mean_PmE| is large.
96 putPmEinXvector = .FALSE.
97 c putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
98
99 IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
100 _BEGIN_MASTER( myThid )
101 ioUnit = standardMessageUnit
102 WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
103 & ' putPmEinXvector =', putPmEinXvector
104 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
105 #ifdef ALLOW_NONHYDROSTATIC
106 WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
107 & ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
108 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
109 WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
110 & ' oldFreeSurfTerm =', oldFreeSurfTerm
111 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
112 #endif
113 _END_MASTER( myThid )
114 ENDIF
115
116 C-- Save previous solution & Initialise Vector solution and source term :
117 sumEmP = 0.
118 DO bj=myByLo(myThid),myByHi(myThid)
119 DO bi=myBxLo(myThid),myBxHi(myThid)
120 DO j=1-OLy,sNy+OLy
121 DO i=1-OLx,sNx+OLx
122 #ifdef ALLOW_CD_CODE
123 etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
124 #endif
125 cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
126 cg2d_b(i,j,bi,bj) = 0.
127 ENDDO
128 ENDDO
129 IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
130 tmpFac = freeSurfFac*mass2rUnit
131 IF (exactConserv)
132 & tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
133 DO j=1,sNy
134 DO i=1,sNx
135 cg2d_b(i,j,bi,bj) =
136 & tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
137 ENDDO
138 ENDDO
139 ENDIF
140 IF ( putPmEinXvector ) THEN
141 tileEmP(bi,bj) = 0.
142 DO j=1,sNy
143 DO i=1,sNx
144 tileEmP(bi,bj) = tileEmP(bi,bj)
145 & + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
146 & *maskInC(i,j,bi,bj)
147 ENDDO
148 ENDDO
149 ENDIF
150 ENDDO
151 ENDDO
152 IF ( putPmEinXvector ) THEN
153 CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
154 ENDIF
155
156 DO bj=myByLo(myThid),myByHi(myThid)
157 DO bi=myBxLo(myThid),myBxHi(myThid)
158 IF ( putPmEinXvector ) THEN
159 tmpFac = 0.
160 IF (globalArea.GT.0.) tmpFac =
161 & freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
162 DO j=1,sNy
163 DO i=1,sNx
164 cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
165 & - tmpFac*Bo_surf(i,j,bi,bj)
166 ENDDO
167 ENDDO
168 ENDIF
169 C- RHS: similar to the divergence of the vertically integrated mass transport:
170 C del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] } / deltaT
171 DO k=Nr,1,-1
172 CALL CALC_DIV_GHAT(
173 I bi,bj,k,
174 U cg2d_b, cg3d_b,
175 I myThid )
176 ENDDO
177 ENDDO
178 ENDDO
179
180 DO bj=myByLo(myThid),myByHi(myThid)
181 DO bi=myBxLo(myThid),myBxHi(myThid)
182 #ifdef ALLOW_NONHYDROSTATIC
183 IF ( oldFreeSurfTerm ) THEN
184 C-- Add source term arising from w=d/dt (p_s + p_nh)
185 DO j=1,sNy
186 DO i=1,sNx
187 ks = kSurfC(i,j,bi,bj)
188 IF ( ks.LE.Nr ) THEN
189 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
190 & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
191 & /deltaTMom/deltaTfreesurf
192 & *( etaN(i,j,bi,bj)
193 & +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
194 cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
195 & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
196 & /deltaTMom/deltaTfreesurf
197 & *( etaN(i,j,bi,bj)
198 & +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
199 ENDIF
200 ENDDO
201 ENDDO
202 ELSEIF ( exactConserv ) THEN
203 #else
204 C-- Add source term arising from w=d/dt (p_s)
205 IF ( exactConserv ) THEN
206 #endif /* ALLOW_NONHYDROSTATIC */
207 DO j=1,sNy
208 DO i=1,sNx
209 ks = kSurfC(i,j,bi,bj)
210 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
211 & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
212 & /deltaTMom/deltaTfreesurf
213 & * etaH(i,j,bi,bj)
214 ENDDO
215 ENDDO
216 ELSE
217 DO j=1,sNy
218 DO i=1,sNx
219 ks = kSurfC(i,j,bi,bj)
220 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
221 & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
222 & /deltaTMom/deltaTfreesurf
223 & * etaN(i,j,bi,bj)
224 ENDDO
225 ENDDO
226 ENDIF
227
228 #ifdef ALLOW_OBCS
229 C- Note: solver matrix is trivial outside OB region (main diagonal only)
230 C => no real need to reset RHS (=cg2d_b) & cg2d_x, except that:
231 C a) normalisation is fct of Max(RHS), which can be large ouside OB region
232 C (would be different if we were solving for increment of eta/g
233 C instead of directly for eta/g).
234 C => need to reset RHS to ensure that interior solution does not depend
235 C on ouside OB region.
236 C b) provide directly the trivial solution cg2d_x == 0 for outside OB region
237 C (=> no residual => no effect on solver convergence and interior solution)
238 IF (useOBCS) THEN
239 DO j=1,sNy
240 DO i=1,sNx
241 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*maskInC(i,j,bi,bj)
242 cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*maskInC(i,j,bi,bj)
243 ENDDO
244 ENDDO
245 ENDIF
246 #endif /* ALLOW_OBCS */
247 C- end bi,bj loops
248 ENDDO
249 ENDDO
250
251 #ifdef ALLOW_DEBUG
252 IF ( debugLevel .GE. debLevD ) THEN
253 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
254 & myThid)
255 ENDIF
256 #endif
257 IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
258 WRITE(sufx,'(I10.10)') myIter
259 CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
260 ENDIF
261
262 C-- Find the surface pressure using a two-dimensional conjugate
263 C-- gradient solver.
264 C see CG2D.h for the interface to this routine.
265 firstResidual=0.
266 lastResidual=0.
267 numIters=cg2dMaxIters
268 c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
269 #ifdef ALLOW_CG2D_NSA
270 C-- Call the not-self-adjoint version of cg2d
271 CALL CG2D_NSA(
272 U cg2d_b,
273 U cg2d_x,
274 O firstResidual,
275 O lastResidual,
276 U numIters,
277 I myThid )
278 #else /* not ALLOW_CG2D_NSA = default */
279 #ifdef ALLOW_SRCG
280 IF ( useSRCGSolver ) THEN
281 C-- Call the single reduce CG solver
282 CALL CG2D_SR(
283 U cg2d_b,
284 U cg2d_x,
285 O firstResidual,
286 O lastResidual,
287 U numIters,
288 I myThid )
289 ELSE
290 #else
291 IF (.TRUE.) THEN
292 C-- Call the default CG solver
293 #endif /* ALLOW_SRCG */
294 CALL CG2D(
295 U cg2d_b,
296 U cg2d_x,
297 O firstResidual,
298 O lastResidual,
299 U numIters,
300 I myThid )
301 ENDIF
302 #endif /* ALLOW_CG2D_NSA */
303 _EXCH_XY_RL( cg2d_x, myThid )
304 c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
305
306 #ifdef ALLOW_DEBUG
307 IF ( debugLevel .GE. debLevD ) THEN
308 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
309 & myThid)
310 ENDIF
311 #endif
312
313 C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
314 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
315 & ) THEN
316 IF ( debugLevel .GE. debLevA ) THEN
317 _BEGIN_MASTER( myThid )
318 WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
319 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320 WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
321 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
322 WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
323 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
324 _END_MASTER( myThid )
325 ENDIF
326 ENDIF
327
328 C-- Transfert the 2D-solution to "etaN" :
329 DO bj=myByLo(myThid),myByHi(myThid)
330 DO bi=myBxLo(myThid),myBxHi(myThid)
331 DO j=1-OLy,sNy+OLy
332 DO i=1-OLx,sNx+OLx
333 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
334 ENDDO
335 ENDDO
336 ENDDO
337 ENDDO
338
339 #ifdef ALLOW_NONHYDROSTATIC
340 IF ( use3Dsolver ) THEN
341 IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
342 WRITE(sufx,'(I10.10)') myIter
343 CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
344 ENDIF
345
346 C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
347 C see CG3D.h for the interface to this routine.
348
349 C-- Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
350 C 2) Update or Add free-surface contribution
351 C 3) increment in horiz velocity due to new cg2d_x
352 C 4) add vertical velocity contribution.
353 CALL PRE_CG3D(
354 I oldFreeSurfTerm,
355 I cg2d_x,
356 U cg3d_b,
357 I myTime, myIter, myThid )
358
359 #ifdef ALLOW_DEBUG
360 IF ( debugLevel .GE. debLevD ) THEN
361 CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
362 & myThid)
363 ENDIF
364 #endif
365 IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
366 WRITE(sufx,'(I10.10)') myIter
367 CALL WRITE_FLD_XYZ_RL('cg3d_b.',sufx, cg3d_b, myIter,myThid )
368 ENDIF
369
370 firstResidual=0.
371 lastResidual=0.
372 numIters=cg3dMaxIters
373 CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
374 CALL CG3D(
375 U cg3d_b,
376 U phi_nh,
377 O firstResidual,
378 O lastResidual,
379 U numIters,
380 I myIter, myThid )
381 _EXCH_XYZ_RL( phi_nh, myThid )
382 CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
383
384 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
385 & ) THEN
386 IF ( debugLevel .GE. debLevA ) THEN
387 _BEGIN_MASTER( myThid )
388 WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
389 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
390 WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
391 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
392 WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
393 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
394 _END_MASTER( myThid )
395 ENDIF
396 ENDIF
397
398 C-- Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
399 C from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
400 IF ( nonHydrostatic .AND. exactConserv ) THEN
401 IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
402 WRITE(sufx,'(I10.10)') myIter
403 CALL WRITE_FLD_XYZ_RL('cg3d_x.',sufx, phi_nh, myIter,myThid )
404 ENDIF
405 CALL POST_CG3D(
406 I zeroPsNH, zeroMeanPnh,
407 I myTime, myIter, myThid )
408 ENDIF
409
410 ENDIF
411 #endif /* ALLOW_NONHYDROSTATIC */
412
413 #ifdef ALLOW_SHOWFLOPS
414 CALL SHOWFLOPS_INSOLVE( myThid)
415 #endif
416
417 RETURN
418 END

  ViewVC Help
Powered by ViewVC 1.1.22