/[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.74 - (show annotations) (download)
Mon Dec 21 00:24:58 2009 UTC (14 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62, checkpoint62b
Changes since 1.73: +2 -2 lines
- use interior masks (instead of maskH, <- to be remove).

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

  ViewVC Help
Powered by ViewVC 1.1.22