/[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.32 - (hide annotations) (download)
Sat Jun 15 03:18:07 2002 UTC (21 years, 11 months ago) by jmc
Branch: MAIN
Changes since 1.31: +17 -9 lines
* reduce cg2d_ output Freq to monitorFreq.

1 jmc 1.32 C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.31 2002/03/27 23:13: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.29 _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     tmpFac = freeSurfFac
76     IF (exactConserv) tmpFac = freeSurfFac*implicDiv2DFlow
77     DO j=1,sNy
78     DO i=1,sNx
79     cg2d_b(i,j,bi,bj) =
80     & tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
81     ENDDO
82     ENDDO
83     ENDIF
84 jmc 1.17 ENDDO
85     ENDDO
86 adcroft 1.12
87     DO bj=myByLo(myThid),myByHi(myThid)
88     DO bi=myBxLo(myThid),myBxHi(myThid)
89     DO K=Nr,1,-1
90     DO j=1,sNy+1
91     DO i=1,sNx+1
92     uf(i,j) = _dyG(i,j,bi,bj)
93     & *drF(k)*_hFacW(i,j,k,bi,bj)
94     vf(i,j) = _dxG(i,j,bi,bj)
95     & *drF(k)*_hFacS(i,j,k,bi,bj)
96     ENDDO
97     ENDDO
98     CALL CALC_DIV_GHAT(
99     I bi,bj,1,sNx,1,sNy,K,
100     I uf,vf,
101 jmc 1.17 U cg2d_b,
102 adcroft 1.12 I myThid)
103     ENDDO
104     ENDDO
105     ENDDO
106 cnh 1.4
107 adcroft 1.12 C-- Add source term arising from w=d/dt (p_s + p_nh)
108     DO bj=myByLo(myThid),myByHi(myThid)
109     DO bi=myBxLo(myThid),myBxHi(myThid)
110 adcroft 1.13 #ifdef ALLOW_NONHYDROSTATIC
111 jmc 1.28 IF ( nonHydrostatic ) THEN
112     DO j=1,sNy
113     DO i=1,sNx
114     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
115     & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
116     & *( etaN(i,j,bi,bj)
117     & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
118     cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
119     & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
120     & *( etaN(i,j,bi,bj)
121     & +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
122     ENDDO
123 adcroft 1.12 ENDDO
124 jmc 1.28 ELSEIF ( exactConserv ) THEN
125 adcroft 1.13 #else
126 jmc 1.26 IF ( exactConserv ) THEN
127 jmc 1.28 #endif
128 jmc 1.26 DO j=1,sNy
129     DO i=1,sNx
130     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
131     & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
132     & * etaH(i,j,bi,bj)
133     ENDDO
134     ENDDO
135     ELSE
136     DO j=1,sNy
137     DO i=1,sNx
138     cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
139     & -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
140     & * etaN(i,j,bi,bj)
141     ENDDO
142 adcroft 1.12 ENDDO
143 jmc 1.26 ENDIF
144 adcroft 1.12
145     #ifdef ALLOW_OBCS
146 adcroft 1.14 IF (useOBCS) THEN
147 adcroft 1.12 DO i=1,sNx
148     C Northern boundary
149     IF (OB_Jn(I,bi,bj).NE.0) THEN
150     cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
151 jmc 1.31 cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
152 adcroft 1.12 ENDIF
153     C Southern boundary
154     IF (OB_Js(I,bi,bj).NE.0) THEN
155     cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
156 jmc 1.31 cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
157 adcroft 1.12 ENDIF
158     ENDDO
159     DO j=1,sNy
160     C Eastern boundary
161     IF (OB_Ie(J,bi,bj).NE.0) THEN
162     cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
163 jmc 1.31 cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
164 adcroft 1.12 ENDIF
165     C Western boundary
166     IF (OB_Iw(J,bi,bj).NE.0) THEN
167     cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
168 jmc 1.31 cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
169 adcroft 1.12 ENDIF
170     ENDDO
171     ENDIF
172     #endif
173     ENDDO
174     ENDDO
175    
176 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
177 adcroft 1.24 IF (debugMode) THEN
178 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
179     & myThid)
180 adcroft 1.24 ENDIF
181 adcroft 1.23 #endif
182 adcroft 1.12
183 cnh 1.1 C-- Find the surface pressure using a two-dimensional conjugate
184     C-- gradient solver.
185 adcroft 1.22 C see CG2D.h for the interface to this routine.
186     firstResidual=0.
187     lastResidual=0.
188 adcroft 1.19 numIters=cg2dMaxIters
189 cnh 1.1 CALL CG2D(
190 adcroft 1.22 U cg2d_b,
191 cnh 1.6 U cg2d_x,
192 adcroft 1.22 O firstResidual,
193     O lastResidual,
194 adcroft 1.19 U numIters,
195 cnh 1.1 I myThid )
196 adcroft 1.19 _EXCH_XY_R8(cg2d_x, myThid )
197 adcroft 1.23
198 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
199 adcroft 1.24 IF (debugMode) THEN
200 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
201     & myThid)
202 adcroft 1.24 ENDIF
203 adcroft 1.23 #endif
204 cnh 1.1
205 jmc 1.32 C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
206     IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
207     & myTime-deltaTClock) ) THEN
208     _BEGIN_MASTER( myThid )
209     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
210     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
211     WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
212     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
213     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
214     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
215     _END_MASTER( )
216     ENDIF
217 jmc 1.17
218     C-- Transfert the 2D-solution to "etaN" :
219     DO bj=myByLo(myThid),myByHi(myThid)
220     DO bi=myBxLo(myThid),myBxHi(myThid)
221     DO j=1-OLy,sNy+OLy
222     DO i=1-OLx,sNx+OLx
223 jmc 1.18 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
224 jmc 1.17 ENDDO
225     ENDDO
226     ENDDO
227     ENDDO
228 adcroft 1.10
229 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
230     IF ( nonHydrostatic ) THEN
231    
232     C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
233     C see CG3D.h for the interface to this routine.
234     DO bj=myByLo(myThid),myByHi(myThid)
235     DO bi=myBxLo(myThid),myBxHi(myThid)
236     DO j=1,sNy+1
237     DO i=1,sNx+1
238 jmc 1.18 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
239 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
240 jmc 1.18 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
241 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
242     ENDDO
243     ENDDO
244    
245 adcroft 1.12 #ifdef ALLOW_OBCS
246 adcroft 1.14 IF (useOBCS) THEN
247 adcroft 1.9 DO i=1,sNx+1
248     C Northern boundary
249     IF (OB_Jn(I,bi,bj).NE.0) THEN
250     vf(I,OB_Jn(I,bi,bj))=0.
251     ENDIF
252     C Southern boundary
253     IF (OB_Js(I,bi,bj).NE.0) THEN
254     vf(I,OB_Js(I,bi,bj)+1)=0.
255     ENDIF
256     ENDDO
257     DO j=1,sNy+1
258     C Eastern boundary
259     IF (OB_Ie(J,bi,bj).NE.0) THEN
260     uf(OB_Ie(J,bi,bj),J)=0.
261     ENDIF
262     C Western boundary
263     IF (OB_Iw(J,bi,bj).NE.0) THEN
264     uf(OB_Iw(J,bi,bj)+1,J)=0.
265     ENDIF
266     ENDDO
267     ENDIF
268 adcroft 1.12 #endif
269 adcroft 1.9
270 adcroft 1.12 K=1
271     DO j=1,sNy
272     DO i=1,sNx
273     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
274     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
275     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
276     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
277     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
278 jmc 1.18 & +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
279     & -wVel(i,j,k+1,bi,bj)
280 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
281     ENDDO
282     ENDDO
283     DO K=2,Nr-1
284 adcroft 1.9 DO j=1,sNy
285     DO i=1,sNx
286     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
287     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
288     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
289     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
290     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
291 adcroft 1.12 & +( wVel(i,j,k ,bi,bj)
292     & -wVel(i,j,k+1,bi,bj)
293     & )*_rA(i,j,bi,bj)/deltaTmom
294    
295 adcroft 1.9 ENDDO
296     ENDDO
297     ENDDO
298 adcroft 1.12 K=Nr
299     DO j=1,sNy
300     DO i=1,sNx
301     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
302     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
303     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
304     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
305     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
306     & +( wVel(i,j,k ,bi,bj)
307     & )*_rA(i,j,bi,bj)/deltaTmom
308    
309     ENDDO
310     ENDDO
311    
312     #ifdef ALLOW_OBCS
313 adcroft 1.14 IF (useOBCS) THEN
314 adcroft 1.12 DO K=1,Nr
315     DO i=1,sNx
316     C Northern boundary
317     IF (OB_Jn(I,bi,bj).NE.0) THEN
318     cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
319     ENDIF
320     C Southern boundary
321     IF (OB_Js(I,bi,bj).NE.0) THEN
322     cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
323     ENDIF
324     ENDDO
325     DO j=1,sNy
326     C Eastern boundary
327     IF (OB_Ie(J,bi,bj).NE.0) THEN
328     cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
329     ENDIF
330     C Western boundary
331     IF (OB_Iw(J,bi,bj).NE.0) THEN
332     cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
333     ENDIF
334     ENDDO
335     ENDDO
336     ENDIF
337     #endif
338 adcroft 1.9
339     ENDDO ! bi
340     ENDDO ! bj
341    
342 adcroft 1.25 firstResidual=0.
343     lastResidual=0.
344     numIters=cg2dMaxIters
345     CALL CG3D(
346     U cg3d_b,
347     U phi_nh,
348     O firstResidual,
349     O lastResidual,
350     U numIters,
351     I myThid )
352     _EXCH_XYZ_R8(phi_nh, myThid )
353    
354     _BEGIN_MASTER( myThid )
355     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
356     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
357     WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
358     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
359     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
360     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
361     _END_MASTER( )
362 adcroft 1.9
363     ENDIF
364     #endif
365 cnh 1.1
366     RETURN
367     END

  ViewVC Help
Powered by ViewVC 1.1.22