/[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.40 - (hide annotations) (download)
Tue Oct 28 22:57:59 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
Changes since 1.39: +2 -2 lines
 o add a "cd_code" package and update all the verification tests
   so that they use the new package instead of "INCLUDE_CD_CODE"

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

  ViewVC Help
Powered by ViewVC 1.1.22