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

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

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


Revision 1.37 - (hide annotations) (download)
Tue May 13 13:30:05 2003 UTC (21 years ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint51, checkpoint50d_post, checkpoint51b_pre, checkpoint50f_post, checkpoint50f_pre, checkpoint51b_post, checkpoint50g_post, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, checkpoint50e_post, checkpoint51a_post
Changes since 1.36: +12 -9 lines
o reduce the output frequency of cg3d-related stuff to the monitor frequency,
  analogous to the cg2d-related output.

1 mlosch 1.37 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.36 2002/10/07 16:20:39 jmc Exp $
2 heimbach 1.21 C $Name: $
3 cnh 1.1
4 adcroft 1.5 #include "CPP_OPTIONS.h"
5 cnh 1.1
6 cnh 1.27 CBOP
7     C !ROUTINE: SOLVE_FOR_PRESSURE
8     C !INTERFACE:
9 jmc 1.29 SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
10 cnh 1.27
11     C !DESCRIPTION: \bv
12     C *==========================================================*
13     C | SUBROUTINE SOLVE_FOR_PRESSURE
14     C | o Controls inversion of two and/or three-dimensional
15     C | elliptic problems for the pressure field.
16     C *==========================================================*
17     C \ev
18    
19     C !USES:
20 adcroft 1.8 IMPLICIT NONE
21 cnh 1.4 C == Global variables
22     #include "SIZE.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "DYNVARS.h"
26 adcroft 1.12 #include "GRID.h"
27 jmc 1.17 #include "SURFACE.h"
28 jmc 1.28 #include "FFIELDS.h"
29 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
30 adcroft 1.25 #include "SOLVE_FOR_PRESSURE3D.h"
31 adcroft 1.9 #include "GW.h"
32 adcroft 1.12 #endif
33 adcroft 1.11 #ifdef ALLOW_OBCS
34 adcroft 1.9 #include "OBCS.h"
35 adcroft 1.11 #endif
36 adcroft 1.22 #include "SOLVE_FOR_PRESSURE.h"
37 cnh 1.4
38 jmc 1.32 C === Functions ====
39     LOGICAL DIFFERENT_MULTIPLE
40     EXTERNAL DIFFERENT_MULTIPLE
41    
42 cnh 1.27 C !INPUT/OUTPUT PARAMETERS:
43 cnh 1.1 C == Routine arguments ==
44 jmc 1.28 C myTime - Current time in simulation
45     C myIter - Current iteration number in simulation
46     C myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
47     _RL myTime
48     INTEGER myIter
49 jmc 1.29 INTEGER myThid
50 cnh 1.4
51 cnh 1.27 C !LOCAL VARIABLES:
52 adcroft 1.22 C == Local variables ==
53 cnh 1.6 INTEGER i,j,k,bi,bj
54 adcroft 1.9 _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55     _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 adcroft 1.22 _RL firstResidual,lastResidual
57 jmc 1.36 _RL tmpFac
58 adcroft 1.19 INTEGER numIters
59 adcroft 1.25 CHARACTER*(MAX_LEN_MBUF) msgBuf
60 cnh 1.27 CEOP
61 jmc 1.17
62     C-- Save previous solution & Initialise Vector solution and source term :
63     DO bj=myByLo(myThid),myByHi(myThid)
64     DO bi=myBxLo(myThid),myBxHi(myThid)
65     DO j=1-OLy,sNy+OLy
66     DO i=1-OLx,sNx+OLx
67 jmc 1.26 #ifdef INCLUDE_CD_CODE
68 jmc 1.17 etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
69 jmc 1.26 #endif
70 jmc 1.18 cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
71 jmc 1.17 cg2d_b(i,j,bi,bj) = 0.
72     ENDDO
73     ENDDO
74 jmc 1.29 IF (useRealFreshWaterFlux) THEN
75 jmc 1.36 tmpFac = freeSurfFac*convertEmP2rUnit
76 mlosch 1.35 IF (exactConserv)
77 jmc 1.36 & tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
78 jmc 1.29 DO j=1,sNy
79     DO i=1,sNx
80     cg2d_b(i,j,bi,bj) =
81     & tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
82     ENDDO
83     ENDDO
84     ENDIF
85 jmc 1.17 ENDDO
86     ENDDO
87 adcroft 1.12
88     DO bj=myByLo(myThid),myByHi(myThid)
89     DO bi=myBxLo(myThid),myBxHi(myThid)
90     DO K=Nr,1,-1
91     DO j=1,sNy+1
92     DO i=1,sNx+1
93     uf(i,j) = _dyG(i,j,bi,bj)
94     & *drF(k)*_hFacW(i,j,k,bi,bj)
95     vf(i,j) = _dxG(i,j,bi,bj)
96     & *drF(k)*_hFacS(i,j,k,bi,bj)
97     ENDDO
98     ENDDO
99     CALL CALC_DIV_GHAT(
100     I bi,bj,1,sNx,1,sNy,K,
101     I uf,vf,
102 jmc 1.17 U cg2d_b,
103 adcroft 1.12 I myThid)
104     ENDDO
105     ENDDO
106     ENDDO
107 cnh 1.4
108 adcroft 1.12 C-- Add source term arising from w=d/dt (p_s + p_nh)
109     DO bj=myByLo(myThid),myByHi(myThid)
110     DO bi=myBxLo(myThid),myBxHi(myThid)
111 adcroft 1.13 #ifdef ALLOW_NONHYDROSTATIC
112 jmc 1.28 IF ( nonHydrostatic ) THEN
113     DO j=1,sNy
114     DO i=1,sNx
115     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
116 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
117 jmc 1.28 & *( etaN(i,j,bi,bj)
118     & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
119     cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
120 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
121 jmc 1.28 & *( etaN(i,j,bi,bj)
122     & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
123     ENDDO
124 adcroft 1.12 ENDDO
125 jmc 1.28 ELSEIF ( exactConserv ) THEN
126 adcroft 1.13 #else
127 jmc 1.26 IF ( exactConserv ) THEN
128 jmc 1.28 #endif
129 jmc 1.26 DO j=1,sNy
130     DO i=1,sNx
131     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
132 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
133 jmc 1.26 & * etaH(i,j,bi,bj)
134     ENDDO
135     ENDDO
136     ELSE
137     DO j=1,sNy
138     DO i=1,sNx
139     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
140 adcroft 1.33 & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
141 jmc 1.26 & * etaN(i,j,bi,bj)
142     ENDDO
143 adcroft 1.12 ENDDO
144 jmc 1.26 ENDIF
145 adcroft 1.12
146     #ifdef ALLOW_OBCS
147 adcroft 1.14 IF (useOBCS) THEN
148 adcroft 1.12 DO i=1,sNx
149     C Northern boundary
150     IF (OB_Jn(I,bi,bj).NE.0) THEN
151     cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
152 jmc 1.31 cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
153 adcroft 1.12 ENDIF
154     C Southern boundary
155     IF (OB_Js(I,bi,bj).NE.0) THEN
156     cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
157 jmc 1.31 cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
158 adcroft 1.12 ENDIF
159     ENDDO
160     DO j=1,sNy
161     C Eastern boundary
162     IF (OB_Ie(J,bi,bj).NE.0) THEN
163     cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
164 jmc 1.31 cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
165 adcroft 1.12 ENDIF
166     C Western boundary
167     IF (OB_Iw(J,bi,bj).NE.0) THEN
168     cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
169 jmc 1.31 cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
170 adcroft 1.12 ENDIF
171     ENDDO
172     ENDIF
173     #endif
174     ENDDO
175     ENDDO
176    
177 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
178 adcroft 1.24 IF (debugMode) THEN
179 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
180     & myThid)
181 adcroft 1.24 ENDIF
182 adcroft 1.23 #endif
183 adcroft 1.12
184 cnh 1.1 C-- Find the surface pressure using a two-dimensional conjugate
185     C-- gradient solver.
186 adcroft 1.22 C see CG2D.h for the interface to this routine.
187     firstResidual=0.
188     lastResidual=0.
189 adcroft 1.19 numIters=cg2dMaxIters
190 cnh 1.1 CALL CG2D(
191 adcroft 1.22 U cg2d_b,
192 cnh 1.6 U cg2d_x,
193 adcroft 1.22 O firstResidual,
194     O lastResidual,
195 adcroft 1.19 U numIters,
196 cnh 1.1 I myThid )
197 adcroft 1.19 _EXCH_XY_R8(cg2d_x, myThid )
198 adcroft 1.23
199 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
200 adcroft 1.24 IF (debugMode) THEN
201 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
202     & myThid)
203 adcroft 1.24 ENDIF
204 adcroft 1.23 #endif
205 cnh 1.1
206 jmc 1.32 C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
207     IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
208     & myTime-deltaTClock) ) THEN
209     _BEGIN_MASTER( myThid )
210     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
211     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
212     WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
213     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
214     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
215     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
216     _END_MASTER( )
217     ENDIF
218 jmc 1.17
219     C-- Transfert the 2D-solution to "etaN" :
220     DO bj=myByLo(myThid),myByHi(myThid)
221     DO bi=myBxLo(myThid),myBxHi(myThid)
222     DO j=1-OLy,sNy+OLy
223     DO i=1-OLx,sNx+OLx
224 jmc 1.18 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
225 jmc 1.17 ENDDO
226     ENDDO
227     ENDDO
228     ENDDO
229 adcroft 1.10
230 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
231     IF ( nonHydrostatic ) THEN
232    
233     C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
234     C see CG3D.h for the interface to this routine.
235     DO bj=myByLo(myThid),myByHi(myThid)
236     DO bi=myBxLo(myThid),myBxHi(myThid)
237     DO j=1,sNy+1
238     DO i=1,sNx+1
239 jmc 1.18 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
240 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
241 jmc 1.18 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
242 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
243     ENDDO
244     ENDDO
245    
246 adcroft 1.12 #ifdef ALLOW_OBCS
247 adcroft 1.14 IF (useOBCS) THEN
248 adcroft 1.9 DO i=1,sNx+1
249     C Northern boundary
250     IF (OB_Jn(I,bi,bj).NE.0) THEN
251     vf(I,OB_Jn(I,bi,bj))=0.
252     ENDIF
253     C Southern boundary
254     IF (OB_Js(I,bi,bj).NE.0) THEN
255     vf(I,OB_Js(I,bi,bj)+1)=0.
256     ENDIF
257     ENDDO
258     DO j=1,sNy+1
259     C Eastern boundary
260     IF (OB_Ie(J,bi,bj).NE.0) THEN
261     uf(OB_Ie(J,bi,bj),J)=0.
262     ENDIF
263     C Western boundary
264     IF (OB_Iw(J,bi,bj).NE.0) THEN
265     uf(OB_Iw(J,bi,bj)+1,J)=0.
266     ENDIF
267     ENDDO
268     ENDIF
269 adcroft 1.12 #endif
270 adcroft 1.9
271 adcroft 1.12 K=1
272     DO j=1,sNy
273     DO i=1,sNx
274     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
275     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
276     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
277     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
278     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
279 jmc 1.18 & +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
280     & -wVel(i,j,k+1,bi,bj)
281 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
282     ENDDO
283     ENDDO
284     DO K=2,Nr-1
285 adcroft 1.9 DO j=1,sNy
286     DO i=1,sNx
287     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
288     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
289     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
290     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
291     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
292 adcroft 1.12 & +( wVel(i,j,k ,bi,bj)
293     & -wVel(i,j,k+1,bi,bj)
294     & )*_rA(i,j,bi,bj)/deltaTmom
295    
296 adcroft 1.9 ENDDO
297     ENDDO
298     ENDDO
299 adcroft 1.12 K=Nr
300     DO j=1,sNy
301     DO i=1,sNx
302     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
303     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
304     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
305     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
306     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
307     & +( wVel(i,j,k ,bi,bj)
308     & )*_rA(i,j,bi,bj)/deltaTmom
309    
310     ENDDO
311     ENDDO
312    
313     #ifdef ALLOW_OBCS
314 adcroft 1.14 IF (useOBCS) THEN
315 adcroft 1.12 DO K=1,Nr
316     DO i=1,sNx
317     C Northern boundary
318     IF (OB_Jn(I,bi,bj).NE.0) THEN
319     cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
320     ENDIF
321     C Southern boundary
322     IF (OB_Js(I,bi,bj).NE.0) THEN
323     cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
324     ENDIF
325     ENDDO
326     DO j=1,sNy
327     C Eastern boundary
328     IF (OB_Ie(J,bi,bj).NE.0) THEN
329     cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
330     ENDIF
331     C Western boundary
332     IF (OB_Iw(J,bi,bj).NE.0) THEN
333     cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
334     ENDIF
335     ENDDO
336     ENDDO
337     ENDIF
338     #endif
339 adcroft 1.9
340     ENDDO ! bi
341     ENDDO ! bj
342    
343 adcroft 1.25 firstResidual=0.
344     lastResidual=0.
345     numIters=cg2dMaxIters
346     CALL CG3D(
347     U cg3d_b,
348     U phi_nh,
349     O firstResidual,
350     O lastResidual,
351     U numIters,
352     I myThid )
353     _EXCH_XYZ_R8(phi_nh, myThid )
354    
355 mlosch 1.37 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
356     & myTime-deltaTClock) ) THEN
357     _BEGIN_MASTER( myThid )
358     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
359     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
360     WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
361     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
362     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
363     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
364     _END_MASTER( )
365     ENDIF
366 adcroft 1.9
367     ENDIF
368     #endif
369 cnh 1.1
370     RETURN
371     END

  ViewVC Help
Powered by ViewVC 1.1.22