/[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.26 - (hide annotations) (download)
Wed Sep 19 13:58:08 2001 UTC (22 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint40
Changes since 1.25: +20 -12 lines
"Volume exact-Conservation" modified for
non-linear free-surface + Crank-Nickelson

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

  ViewVC Help
Powered by ViewVC 1.1.22