/[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.34 - (hide annotations) (download)
Mon Sep 23 16:13:31 2002 UTC (21 years, 8 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint46h_pre
Changes since 1.33: +4 -1 lines
Added code to convert surface volume flux (fresh water) into
a mass flux when using P coordinates in the ocean (OCEANICP).
Note: It assumes you have set rho0=rhoConst=density of fresh water.

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

  ViewVC Help
Powered by ViewVC 1.1.22