/[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.38 - (hide annotations) (download)
Tue Jul 8 15:00:26 2003 UTC (20 years, 10 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51f_post, checkpoint51d_post, checkpoint51j_post, checkpoint51h_pre, branchpoint-genmake2, checkpoint51c_post, checkpoint51i_pre, checkpoint51e_post, checkpoint51f_pre, checkpoint51g_post
Branch point for: branch-genmake2
Changes since 1.37: +23 -19 lines
o introducing integer flag debugLevel
o introducing pathname variable mdsioLocalDir for mdsio

1 heimbach 1.38 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.27.6.2 2003/06/24 23:05:29 heimbach 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 heimbach 1.38 IF ( debugLevel .GE. debLevB ) 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 heimbach 1.38 IF ( debugLevel .GE. debLevB ) 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 heimbach 1.38 IF ( debugLevel .GE. debLevA ) THEN
210     _BEGIN_MASTER( myThid )
211     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
212     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
213     WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
214     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
215     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
216     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
217     _END_MASTER( )
218     ENDIF
219 jmc 1.32 ENDIF
220 jmc 1.17
221     C-- Transfert the 2D-solution to "etaN" :
222     DO bj=myByLo(myThid),myByHi(myThid)
223     DO bi=myBxLo(myThid),myBxHi(myThid)
224     DO j=1-OLy,sNy+OLy
225     DO i=1-OLx,sNx+OLx
226 jmc 1.18 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
227 jmc 1.17 ENDDO
228     ENDDO
229     ENDDO
230     ENDDO
231 adcroft 1.10
232 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
233     IF ( nonHydrostatic ) THEN
234    
235     C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
236     C see CG3D.h for the interface to this routine.
237     DO bj=myByLo(myThid),myByHi(myThid)
238     DO bi=myBxLo(myThid),myBxHi(myThid)
239     DO j=1,sNy+1
240     DO i=1,sNx+1
241 jmc 1.18 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
242 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
243 jmc 1.18 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
244 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
245     ENDDO
246     ENDDO
247    
248 adcroft 1.12 #ifdef ALLOW_OBCS
249 adcroft 1.14 IF (useOBCS) THEN
250 adcroft 1.9 DO i=1,sNx+1
251     C Northern boundary
252     IF (OB_Jn(I,bi,bj).NE.0) THEN
253     vf(I,OB_Jn(I,bi,bj))=0.
254     ENDIF
255     C Southern boundary
256     IF (OB_Js(I,bi,bj).NE.0) THEN
257     vf(I,OB_Js(I,bi,bj)+1)=0.
258     ENDIF
259     ENDDO
260     DO j=1,sNy+1
261     C Eastern boundary
262     IF (OB_Ie(J,bi,bj).NE.0) THEN
263     uf(OB_Ie(J,bi,bj),J)=0.
264     ENDIF
265     C Western boundary
266     IF (OB_Iw(J,bi,bj).NE.0) THEN
267     uf(OB_Iw(J,bi,bj)+1,J)=0.
268     ENDIF
269     ENDDO
270     ENDIF
271 adcroft 1.12 #endif
272 adcroft 1.9
273 adcroft 1.12 K=1
274     DO j=1,sNy
275     DO i=1,sNx
276     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
277     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
278     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
279     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
280     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
281 jmc 1.18 & +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
282     & -wVel(i,j,k+1,bi,bj)
283 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
284     ENDDO
285     ENDDO
286     DO K=2,Nr-1
287 adcroft 1.9 DO j=1,sNy
288     DO i=1,sNx
289     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
290     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
291     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
292     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
293     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
294 adcroft 1.12 & +( wVel(i,j,k ,bi,bj)
295     & -wVel(i,j,k+1,bi,bj)
296     & )*_rA(i,j,bi,bj)/deltaTmom
297    
298 adcroft 1.9 ENDDO
299     ENDDO
300     ENDDO
301 adcroft 1.12 K=Nr
302     DO j=1,sNy
303     DO i=1,sNx
304     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
305     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
306     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
307     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
308     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
309     & +( wVel(i,j,k ,bi,bj)
310     & )*_rA(i,j,bi,bj)/deltaTmom
311    
312     ENDDO
313     ENDDO
314    
315     #ifdef ALLOW_OBCS
316 adcroft 1.14 IF (useOBCS) THEN
317 adcroft 1.12 DO K=1,Nr
318     DO i=1,sNx
319     C Northern boundary
320     IF (OB_Jn(I,bi,bj).NE.0) THEN
321     cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
322     ENDIF
323     C Southern boundary
324     IF (OB_Js(I,bi,bj).NE.0) THEN
325     cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
326     ENDIF
327     ENDDO
328     DO j=1,sNy
329     C Eastern boundary
330     IF (OB_Ie(J,bi,bj).NE.0) THEN
331     cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
332     ENDIF
333     C Western boundary
334     IF (OB_Iw(J,bi,bj).NE.0) THEN
335     cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
336     ENDIF
337     ENDDO
338     ENDDO
339     ENDIF
340     #endif
341 adcroft 1.9
342     ENDDO ! bi
343     ENDDO ! bj
344    
345 adcroft 1.25 firstResidual=0.
346     lastResidual=0.
347     numIters=cg2dMaxIters
348     CALL CG3D(
349     U cg3d_b,
350     U phi_nh,
351     O firstResidual,
352     O lastResidual,
353     U numIters,
354     I myThid )
355     _EXCH_XYZ_R8(phi_nh, myThid )
356    
357 mlosch 1.37 IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
358     & myTime-deltaTClock) ) THEN
359 heimbach 1.38 IF ( debugLevel .GE. debLevA ) THEN
360     _BEGIN_MASTER( myThid )
361     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
362     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
363     WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
364     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
365     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
366     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
367     _END_MASTER( )
368     ENDIF
369 mlosch 1.37 ENDIF
370 adcroft 1.9
371     ENDIF
372     #endif
373 cnh 1.1
374     RETURN
375     END

  ViewVC Help
Powered by ViewVC 1.1.22