/[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.47 - (show annotations) (download)
Thu Jul 21 19:47:53 2005 UTC (18 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57m_post, checkpoint57s_post, checkpoint57r_post, checkpoint57n_post, checkpoint57t_post, checkpoint57v_post, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, checkpoint57o_post, checkpoint57w_post
Changes since 1.46: +34 -1 lines
add global-mean PmE to the initial cg2d X-vector ; turned off for now.

1 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.46 2005/05/15 03:02:08 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 "DYNVARS.h"
27 #ifdef ALLOW_CD_CODE
28 #include "CD_CODE_VARS.h"
29 #endif
30 #include "GRID.h"
31 #include "SURFACE.h"
32 #include "FFIELDS.h"
33 #ifdef ALLOW_NONHYDROSTATIC
34 #include "SOLVE_FOR_PRESSURE3D.h"
35 #include "GW.h"
36 #endif
37 #ifdef ALLOW_OBCS
38 #include "OBCS.h"
39 #endif
40 #include "SOLVE_FOR_PRESSURE.h"
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 _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59 _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 _RL firstResidual,lastResidual
61 _RL tmpFac
62 _RL sumEmP, tileEmP
63 LOGICAL putPmEinXvector
64 INTEGER numIters
65 CHARACTER*(MAX_LEN_MBUF) msgBuf
66 CEOP
67
68 #ifdef TIME_PER_TIMESTEP
69 CCE107 common block for per timestep timing
70 C !TIMING VARIABLES
71 C == Timing variables ==
72 REAL*8 utnew, utold, stnew, stold, wtnew, wtold
73 COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
74 #endif
75
76 C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
77 C instead of simply etaN ; This can speed-up the solver convergence in
78 C the case where |Global_mean_PmE| is large.
79 putPmEinXvector = .FALSE.
80 c putPmEinXvector = useRealFreshWaterFlux
81
82 C-- Save previous solution & Initialise Vector solution and source term :
83 sumEmP = 0.
84 DO bj=myByLo(myThid),myByHi(myThid)
85 DO bi=myBxLo(myThid),myBxHi(myThid)
86 DO j=1-OLy,sNy+OLy
87 DO i=1-OLx,sNx+OLx
88 #ifdef ALLOW_CD_CODE
89 etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
90 #endif
91 cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
92 cg2d_b(i,j,bi,bj) = 0.
93 ENDDO
94 ENDDO
95 IF (useRealFreshWaterFlux) THEN
96 tmpFac = freeSurfFac*convertEmP2rUnit
97 IF (exactConserv)
98 & tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
99 DO j=1,sNy
100 DO i=1,sNx
101 cg2d_b(i,j,bi,bj) =
102 & tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
103 ENDDO
104 ENDDO
105 ENDIF
106 IF ( putPmEinXvector ) THEN
107 tileEmP = 0.
108 DO j=1,sNy
109 DO i=1,sNx
110 tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
111 & *maskH(i,j,bi,bj)
112 ENDDO
113 ENDDO
114 sumEmP = sumEmP + tileEmP
115 ENDIF
116 ENDDO
117 ENDDO
118 IF ( putPmEinXvector ) THEN
119 _GLOBAL_SUM_R8( sumEmP, myThid )
120 ENDIF
121
122 DO bj=myByLo(myThid),myByHi(myThid)
123 DO bi=myBxLo(myThid),myBxHi(myThid)
124 IF ( putPmEinXvector ) THEN
125 tmpFac = 0.
126 IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
127 & *convertEmP2rUnit*sumEmP/globalArea
128 DO j=1,sNy
129 DO i=1,sNx
130 cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
131 & - tmpFac*Bo_surf(i,j,bi,bj)
132 ENDDO
133 ENDDO
134 ENDIF
135 DO K=Nr,1,-1
136 DO j=1,sNy+1
137 DO i=1,sNx+1
138 uf(i,j) = _dyG(i,j,bi,bj)
139 & *drF(k)*_hFacW(i,j,k,bi,bj)
140 vf(i,j) = _dxG(i,j,bi,bj)
141 & *drF(k)*_hFacS(i,j,k,bi,bj)
142 ENDDO
143 ENDDO
144 CALL CALC_DIV_GHAT(
145 I bi,bj,1,sNx,1,sNy,K,
146 I uf,vf,
147 U cg2d_b,
148 I myThid)
149 ENDDO
150 ENDDO
151 ENDDO
152
153 C-- Add source term arising from w=d/dt (p_s + p_nh)
154 DO bj=myByLo(myThid),myByHi(myThid)
155 DO bi=myBxLo(myThid),myBxHi(myThid)
156 #ifdef ALLOW_NONHYDROSTATIC
157 IF ( nonHydrostatic ) THEN
158 DO j=1,sNy
159 DO i=1,sNx
160 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
161 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
162 & *( etaN(i,j,bi,bj)
163 & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
164 cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
165 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
166 & *( etaN(i,j,bi,bj)
167 & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
168 ENDDO
169 ENDDO
170 ELSEIF ( exactConserv ) THEN
171 #else
172 IF ( exactConserv ) THEN
173 #endif /* ALLOW_NONHYDROSTATIC */
174 DO j=1,sNy
175 DO i=1,sNx
176 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
177 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
178 & * etaH(i,j,bi,bj)
179 ENDDO
180 ENDDO
181 ELSE
182 DO j=1,sNy
183 DO i=1,sNx
184 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
185 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
186 & * etaN(i,j,bi,bj)
187 ENDDO
188 ENDDO
189 ENDIF
190
191 #ifdef ALLOW_OBCS
192 IF (useOBCS) THEN
193 DO i=1,sNx
194 C Northern boundary
195 IF (OB_Jn(I,bi,bj).NE.0) THEN
196 cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
197 cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
198 ENDIF
199 C Southern boundary
200 IF (OB_Js(I,bi,bj).NE.0) THEN
201 cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
202 cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
203 ENDIF
204 ENDDO
205 DO j=1,sNy
206 C Eastern boundary
207 IF (OB_Ie(J,bi,bj).NE.0) THEN
208 cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
209 cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
210 ENDIF
211 C Western boundary
212 IF (OB_Iw(J,bi,bj).NE.0) THEN
213 cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
214 cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
215 ENDIF
216 ENDDO
217 ENDIF
218 #endif
219 ENDDO
220 ENDDO
221
222 #ifdef ALLOW_DEBUG
223 IF ( debugLevel .GE. debLevB ) THEN
224 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
225 & myThid)
226 ENDIF
227 #endif
228
229 C-- Find the surface pressure using a two-dimensional conjugate
230 C-- gradient solver.
231 C see CG2D.h for the interface to this routine.
232 firstResidual=0.
233 lastResidual=0.
234 numIters=cg2dMaxIters
235 CALL CG2D(
236 U cg2d_b,
237 U cg2d_x,
238 O firstResidual,
239 O lastResidual,
240 U numIters,
241 I myThid )
242 _EXCH_XY_R8(cg2d_x, myThid )
243
244 #ifdef ALLOW_DEBUG
245 IF ( debugLevel .GE. debLevB ) THEN
246 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
247 & myThid)
248 ENDIF
249 #endif
250
251 C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
252 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
253 & ) THEN
254 IF ( debugLevel .GE. debLevA ) THEN
255 _BEGIN_MASTER( myThid )
256 WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
257 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
258 WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
259 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
260 WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
261 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
262 _END_MASTER( myThid )
263 ENDIF
264 ENDIF
265
266 C-- Transfert the 2D-solution to "etaN" :
267 DO bj=myByLo(myThid),myByHi(myThid)
268 DO bi=myBxLo(myThid),myBxHi(myThid)
269 DO j=1-OLy,sNy+OLy
270 DO i=1-OLx,sNx+OLx
271 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
272 ENDDO
273 ENDDO
274 ENDDO
275 ENDDO
276
277 #ifdef ALLOW_NONHYDROSTATIC
278 IF ( nonHydrostatic ) THEN
279
280 C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
281 C see CG3D.h for the interface to this routine.
282 DO bj=myByLo(myThid),myByHi(myThid)
283 DO bi=myBxLo(myThid),myBxHi(myThid)
284 DO j=1,sNy+1
285 DO i=1,sNx+1
286 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
287 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
288 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
289 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
290 ENDDO
291 ENDDO
292
293 #ifdef ALLOW_OBCS
294 IF (useOBCS) THEN
295 DO i=1,sNx+1
296 C Northern boundary
297 IF (OB_Jn(I,bi,bj).NE.0) THEN
298 vf(I,OB_Jn(I,bi,bj))=0.
299 ENDIF
300 C Southern boundary
301 IF (OB_Js(I,bi,bj).NE.0) THEN
302 vf(I,OB_Js(I,bi,bj)+1)=0.
303 ENDIF
304 ENDDO
305 DO j=1,sNy+1
306 C Eastern boundary
307 IF (OB_Ie(J,bi,bj).NE.0) THEN
308 uf(OB_Ie(J,bi,bj),J)=0.
309 ENDIF
310 C Western boundary
311 IF (OB_Iw(J,bi,bj).NE.0) THEN
312 uf(OB_Iw(J,bi,bj)+1,J)=0.
313 ENDIF
314 ENDDO
315 ENDIF
316 #endif
317
318 K=1
319 DO j=1,sNy
320 DO i=1,sNx
321 cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
322 & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
323 & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
324 & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
325 & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
326 & +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
327 & -wVel(i,j,k+1,bi,bj)
328 & )*_rA(i,j,bi,bj)/deltaTmom
329 ENDDO
330 ENDDO
331 DO K=2,Nr-1
332 DO j=1,sNy
333 DO i=1,sNx
334 cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
335 & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
336 & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
337 & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
338 & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
339 & +( wVel(i,j,k ,bi,bj)
340 & -wVel(i,j,k+1,bi,bj)
341 & )*_rA(i,j,bi,bj)/deltaTmom
342
343 ENDDO
344 ENDDO
345 ENDDO
346 K=Nr
347 DO j=1,sNy
348 DO i=1,sNx
349 cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
350 & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
351 & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
352 & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
353 & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
354 & +( wVel(i,j,k ,bi,bj)
355 & )*_rA(i,j,bi,bj)/deltaTmom
356
357 ENDDO
358 ENDDO
359
360 #ifdef ALLOW_OBCS
361 IF (useOBCS) THEN
362 DO K=1,Nr
363 DO i=1,sNx
364 C Northern boundary
365 IF (OB_Jn(I,bi,bj).NE.0) THEN
366 cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
367 ENDIF
368 C Southern boundary
369 IF (OB_Js(I,bi,bj).NE.0) THEN
370 cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
371 ENDIF
372 ENDDO
373 DO j=1,sNy
374 C Eastern boundary
375 IF (OB_Ie(J,bi,bj).NE.0) THEN
376 cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
377 ENDIF
378 C Western boundary
379 IF (OB_Iw(J,bi,bj).NE.0) THEN
380 cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
381 ENDIF
382 ENDDO
383 ENDDO
384 ENDIF
385 #endif
386
387 ENDDO ! bi
388 ENDDO ! bj
389
390 firstResidual=0.
391 lastResidual=0.
392 numIters=cg2dMaxIters
393 CALL CG3D(
394 U cg3d_b,
395 U phi_nh,
396 O firstResidual,
397 O lastResidual,
398 U numIters,
399 I myThid )
400 _EXCH_XYZ_R8(phi_nh, myThid )
401
402 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
403 & ) THEN
404 IF ( debugLevel .GE. debLevA ) THEN
405 _BEGIN_MASTER( myThid )
406 WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
407 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
408 WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
409 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
410 WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
411 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
412 _END_MASTER( myThid )
413 ENDIF
414 ENDIF
415
416 ENDIF
417 #endif
418
419 #ifdef TIME_PER_TIMESTEP
420 CCE107 Time per timestep information
421 _BEGIN_MASTER( myThid )
422 CALL TIMER_GET_TIME( utnew, stnew, wtnew )
423 C Only output timing information after the 1st timestep
424 IF ( wtold .NE. 0.0D0 ) THEN
425 WRITE(msgBuf,'(A34,3F10.6)')
426 $ 'User, system and wallclock time:', utnew - utold,
427 $ stnew - stold, wtnew - wtold
428 CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
429 ENDIF
430 utold = utnew
431 stold = stnew
432 wtold = wtnew
433 _END_MASTER( myThid )
434 #endif
435
436 RETURN
437 END
438
439 #ifdef TIME_PER_TIMESTEP
440 CCE107 Initialization of common block for per timestep timing
441 BLOCK DATA settimers
442 C !TIMING VARIABLES
443 C == Timing variables ==
444 REAL*8 utnew, utold, stnew, stold, wtnew, wtold
445 COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
446 DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
447 END
448 #endif

  ViewVC Help
Powered by ViewVC 1.1.22