/[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.27 - (hide annotations) (download)
Wed Sep 26 18:09:16 2001 UTC (22 years, 8 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint43a-release1mods, release1_b1, checkpoint43, icebear5, icebear4, icebear3, icebear2, release1-branch_tutorials, chkpt44a_post, chkpt44c_pre, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, release1-branch-end, checkpoint44b_post, ecco_ice2, ecco_ice1, ecco_c44_e22, ecco_c44_e25, chkpt44a_pre, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5, release1_beta1, checkpoint44b_pre, checkpoint42, checkpoint41, checkpoint44, release1-branch_branchpoint
Branch point for: c24_e25_ice, release1-branch, release1, ecco-branch, icebear, release1_coupled
Changes since 1.26: +17 -9 lines
Bringing comments up to data and formatting for document extraction.

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

  ViewVC Help
Powered by ViewVC 1.1.22