/[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.30 - (hide annotations) (download)
Mon Mar 4 17:26:41 2002 UTC (22 years, 2 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint44f_post, checkpoint44h_pre, checkpoint44g_post, checkpoint44h_post, checkpoint45
Changes since 1.29: +3 -3 lines 
Added PTRACERS package

This allows an arbitrary number of passive tracers to be integrated
forward simultaneously with the dynamicaly model.
 + Implemented so far:
    - basic forward algorithm (time-stepping, advection, diffusion, convection)
    - I/O and checkpointing
    - GM/Redi  *but*  using the GM/Redi coefficient of Salt
 + Not implemented so far:
    - KPP
    - OBCS
 + No specific example supplied (yet) but global_ocean.90x40x15 has the
   necessary data.ptracer file. Simply use -enable=ptracers and uncomment
   line in data.pkg. PTRACER01 then reproduces Salt exactly.
 + This package is disabled by default since it increases storage.
1 adcroft 1.30 C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.29 2002/02/10 00:39:22 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 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     ENDIF
148     C Southern boundary
149     IF (OB_Js(I,bi,bj).NE.0) THEN
150     cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
151     ENDIF
152     ENDDO
153     DO j=1,sNy
154     C Eastern boundary
155     IF (OB_Ie(J,bi,bj).NE.0) THEN
156     cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
157     ENDIF
158     C Western boundary
159     IF (OB_Iw(J,bi,bj).NE.0) THEN
160     cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
161     ENDIF
162     ENDDO
163     ENDIF
164     #endif
165     ENDDO
166     ENDDO
167    
168 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
169 adcroft 1.24 IF (debugMode) THEN
170 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
171     & myThid)
172 adcroft 1.24 ENDIF
173 adcroft 1.23 #endif
174 adcroft 1.12
175 cnh 1.1 C-- Find the surface pressure using a two-dimensional conjugate
176     C-- gradient solver.
177 adcroft 1.22 C see CG2D.h for the interface to this routine.
178     firstResidual=0.
179     lastResidual=0.
180 adcroft 1.19 numIters=cg2dMaxIters
181 cnh 1.1 CALL CG2D(
182 adcroft 1.22 U cg2d_b,
183 cnh 1.6 U cg2d_x,
184 adcroft 1.22 O firstResidual,
185     O lastResidual,
186 adcroft 1.19 U numIters,
187 cnh 1.1 I myThid )
188 adcroft 1.19 _EXCH_XY_R8(cg2d_x, myThid )
189 adcroft 1.23
190 adcroft 1.30 #ifndef DISABLE_DEBUGMODE
191 adcroft 1.24 IF (debugMode) THEN
192 adcroft 1.23 CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
193     & myThid)
194 adcroft 1.24 ENDIF
195 adcroft 1.23 #endif
196 cnh 1.1
197 adcroft 1.19 _BEGIN_MASTER( myThid )
198 adcroft 1.25 WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
199     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
200     WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
201     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
202     WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
203     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
204 adcroft 1.19 _END_MASTER( )
205 jmc 1.17
206     C-- Transfert the 2D-solution to "etaN" :
207     DO bj=myByLo(myThid),myByHi(myThid)
208     DO bi=myBxLo(myThid),myBxHi(myThid)
209     DO j=1-OLy,sNy+OLy
210     DO i=1-OLx,sNx+OLx
211 jmc 1.18 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
212 jmc 1.17 ENDDO
213     ENDDO
214     ENDDO
215     ENDDO
216 adcroft 1.10
217 adcroft 1.9 #ifdef ALLOW_NONHYDROSTATIC
218     IF ( nonHydrostatic ) THEN
219    
220     C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
221     C see CG3D.h for the interface to this routine.
222     DO bj=myByLo(myThid),myByHi(myThid)
223     DO bi=myBxLo(myThid),myBxHi(myThid)
224     DO j=1,sNy+1
225     DO i=1,sNx+1
226 jmc 1.18 uf(i,j)=-_recip_dxC(i,j,bi,bj)*
227 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
228 jmc 1.18 vf(i,j)=-_recip_dyC(i,j,bi,bj)*
229 adcroft 1.9 & (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
230     ENDDO
231     ENDDO
232    
233 adcroft 1.12 #ifdef ALLOW_OBCS
234 adcroft 1.14 IF (useOBCS) THEN
235 adcroft 1.9 DO i=1,sNx+1
236     C Northern boundary
237     IF (OB_Jn(I,bi,bj).NE.0) THEN
238     vf(I,OB_Jn(I,bi,bj))=0.
239     ENDIF
240     C Southern boundary
241     IF (OB_Js(I,bi,bj).NE.0) THEN
242     vf(I,OB_Js(I,bi,bj)+1)=0.
243     ENDIF
244     ENDDO
245     DO j=1,sNy+1
246     C Eastern boundary
247     IF (OB_Ie(J,bi,bj).NE.0) THEN
248     uf(OB_Ie(J,bi,bj),J)=0.
249     ENDIF
250     C Western boundary
251     IF (OB_Iw(J,bi,bj).NE.0) THEN
252     uf(OB_Iw(J,bi,bj)+1,J)=0.
253     ENDIF
254     ENDDO
255     ENDIF
256 adcroft 1.12 #endif
257 adcroft 1.9
258 adcroft 1.12 K=1
259     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 jmc 1.18 & +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
267     & -wVel(i,j,k+1,bi,bj)
268 adcroft 1.12 & )*_rA(i,j,bi,bj)/deltaTmom
269     ENDDO
270     ENDDO
271     DO K=2,Nr-1
272 adcroft 1.9 DO j=1,sNy
273     DO i=1,sNx
274     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
275     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
276     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
277     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
278     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
279 adcroft 1.12 & +( wVel(i,j,k ,bi,bj)
280     & -wVel(i,j,k+1,bi,bj)
281     & )*_rA(i,j,bi,bj)/deltaTmom
282    
283 adcroft 1.9 ENDDO
284     ENDDO
285     ENDDO
286 adcroft 1.12 K=Nr
287     DO j=1,sNy
288     DO i=1,sNx
289     cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
290     & +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
291     & -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
292     & +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
293     & -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
294     & +( wVel(i,j,k ,bi,bj)
295     & )*_rA(i,j,bi,bj)/deltaTmom
296    
297     ENDDO
298     ENDDO
299    
300     #ifdef ALLOW_OBCS
301 adcroft 1.14 IF (useOBCS) THEN
302 adcroft 1.12 DO K=1,Nr
303     DO i=1,sNx
304     C Northern boundary
305     IF (OB_Jn(I,bi,bj).NE.0) THEN
306     cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
307     ENDIF
308     C Southern boundary
309     IF (OB_Js(I,bi,bj).NE.0) THEN
310     cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
311     ENDIF
312     ENDDO
313     DO j=1,sNy
314     C Eastern boundary
315     IF (OB_Ie(J,bi,bj).NE.0) THEN
316     cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
317     ENDIF
318     C Western boundary
319     IF (OB_Iw(J,bi,bj).NE.0) THEN
320     cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
321     ENDIF
322     ENDDO
323     ENDDO
324     ENDIF
325     #endif
326 adcroft 1.9
327     ENDDO ! bi
328     ENDDO ! bj
329    
330 adcroft 1.25 firstResidual=0.
331     lastResidual=0.
332     numIters=cg2dMaxIters
333     CALL CG3D(
334     U cg3d_b,
335     U phi_nh,
336     O firstResidual,
337     O lastResidual,
338     U numIters,
339     I myThid )
340     _EXCH_XYZ_R8(phi_nh, myThid )
341    
342     _BEGIN_MASTER( myThid )
343     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
344     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
345     WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
346     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
347     WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
348     CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
349     _END_MASTER( )
350 adcroft 1.9
351     ENDIF
352     #endif
353 cnh 1.1
354     RETURN
355     END
  ViewVC Help
Powered by ViewVC 1.1.22