/[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.31 - (hide annotations) (download)
Wed Mar 27 23:13:39 2002 UTC (22 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint45a_post, checkpoint45b_post, checkpoint45c_post
Changes since 1.30: +5 -1 lines
set cg2d_x to zero where OBCS are applied : this avoid large initial
 cg2d-residual when using OBCS & exactConserv

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

  ViewVC Help
Powered by ViewVC 1.1.22