/[MITgcm]/MITgcm/model/src/solve_for_pressure.F
ViewVC logotype

Diff of /MITgcm/model/src/solve_for_pressure.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by cnh, Mon Jun 8 21:43:02 1998 UTC revision 1.32 by jmc, Sat Jun 15 03:18:07 2002 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )  C     !ROUTINE: SOLVE_FOR_PRESSURE
8  C     /==========================================================\  C     !INTERFACE:
9  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
 C     | o Controls inversion of two and/or three-dimensional     |  
 C     |   elliptic problems for the pressure field.              |  
 C     \==========================================================/  
10    
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          IMPLICIT NONE
21  C     == Global variables  C     == Global variables
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "DYNVARS.h"  #include "DYNVARS.h"
26  #include "CG2D.h"  #include "GRID.h"
27    #include "SURFACE.h"
28    #include "FFIELDS.h"
29    #ifdef ALLOW_NONHYDROSTATIC
30    #include "SOLVE_FOR_PRESSURE3D.h"
31    #include "GW.h"
32    #endif
33    #ifdef ALLOW_OBCS
34    #include "OBCS.h"
35    #endif
36    #include "SOLVE_FOR_PRESSURE.h"
37    
38    C     === Functions ====
39          LOGICAL  DIFFERENT_MULTIPLE
40          EXTERNAL DIFFERENT_MULTIPLE
41    
42    C     !INPUT/OUTPUT PARAMETERS:
43  C     == Routine arguments ==  C     == Routine arguments ==
44  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  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        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
50    
51  C     Local variables  C     !LOCAL VARIABLES:
52        INTEGER i,j,bi,bj  C     == Local variables ==
53          INTEGER i,j,k,bi,bj
54          _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL firstResidual,lastResidual
57          _RL tmpFac
58          INTEGER numIters
59          CHARACTER*(MAX_LEN_MBUF) msgBuf
60    CEOP
61    
62  #ifdef ALLOW_CD  C--   Save previous solution & Initialise Vector solution and source term :
 C--   Save previous solution.  
63        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
64         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
65          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
66           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
67            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  #ifdef INCLUDE_CD_CODE
68              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
69    #endif
70              cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
71              cg2d_b(i,j,bi,bj) = 0.
72           ENDDO           ENDDO
73          ENDDO          ENDDO
74            IF (useRealFreshWaterFlux) THEN
75             tmpFac = freeSurfFac
76             IF (exactConserv) tmpFac = freeSurfFac*implicDiv2DFlow
77             DO j=1,sNy
78              DO i=1,sNx
79               cg2d_b(i,j,bi,bj) =
80         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
81              ENDDO
82             ENDDO
83            ENDIF
84         ENDDO         ENDDO
85        ENDDO        ENDDO
86    
87          DO bj=myByLo(myThid),myByHi(myThid)
88           DO bi=myBxLo(myThid),myBxHi(myThid)
89            DO K=Nr,1,-1
90             DO j=1,sNy+1
91              DO i=1,sNx+1
92               uf(i,j) = _dyG(i,j,bi,bj)
93         &      *drF(k)*_hFacW(i,j,k,bi,bj)
94               vf(i,j) = _dxG(i,j,bi,bj)
95         &      *drF(k)*_hFacS(i,j,k,bi,bj)
96              ENDDO
97             ENDDO
98             CALL CALC_DIV_GHAT(
99         I       bi,bj,1,sNx,1,sNy,K,
100         I       uf,vf,
101         U       cg2d_b,
102         I       myThid)
103            ENDDO
104           ENDDO
105          ENDDO
106    
107    C--   Add source term arising from w=d/dt (p_s + p_nh)
108          DO bj=myByLo(myThid),myByHi(myThid)
109           DO bi=myBxLo(myThid),myBxHi(myThid)
110    #ifdef ALLOW_NONHYDROSTATIC
111            IF ( nonHydrostatic ) THEN
112             DO j=1,sNy
113              DO i=1,sNx
114               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,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               cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
119         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
120         &         *( etaN(i,j,bi,bj)
121         &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
122              ENDDO
123             ENDDO
124            ELSEIF ( exactConserv ) THEN
125    #else
126            IF ( exactConserv ) THEN
127    #endif
128             DO j=1,sNy
129              DO i=1,sNx
130               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
131         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
132         &         * etaH(i,j,bi,bj)
133              ENDDO
134             ENDDO
135            ELSE
136             DO j=1,sNy
137              DO i=1,sNx
138               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
139         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
140         &         * etaN(i,j,bi,bj)
141              ENDDO
142             ENDDO
143            ENDIF
144    
145    #ifdef ALLOW_OBCS
146            IF (useOBCS) THEN
147             DO i=1,sNx
148    C Northern boundary
149              IF (OB_Jn(I,bi,bj).NE.0) THEN
150               cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
151               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
152              ENDIF
153    C Southern boundary
154              IF (OB_Js(I,bi,bj).NE.0) THEN
155               cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
156               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
157              ENDIF
158             ENDDO
159             DO j=1,sNy
160    C Eastern boundary
161              IF (OB_Ie(J,bi,bj).NE.0) THEN
162               cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
163               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
164              ENDIF
165    C Western boundary
166              IF (OB_Iw(J,bi,bj).NE.0) THEN
167               cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
168               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
169              ENDIF
170             ENDDO
171            ENDIF
172    #endif
173           ENDDO
174          ENDDO
175    
176    #ifndef DISABLE_DEBUGMODE
177          IF (debugMode) THEN
178           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
179         &                        myThid)
180          ENDIF
181  #endif  #endif
182    
183  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
184  C--   gradient solver.  C--   gradient solver.
185  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
186          firstResidual=0.
187          lastResidual=0.
188          numIters=cg2dMaxIters
189        CALL CG2D(        CALL CG2D(
190         U           cg2d_b,
191         U           cg2d_x,
192         O           firstResidual,
193         O           lastResidual,
194         U           numIters,
195       I           myThid )       I           myThid )
196          _EXCH_XY_R8(cg2d_x, myThid )
197    
198    #ifndef DISABLE_DEBUGMODE
199          IF (debugMode) THEN
200           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
201         &                        myThid)
202          ENDIF
203    #endif
204    
205    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
206          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
207         &                                    myTime-deltaTClock) ) THEN
208           _BEGIN_MASTER( myThid )
209           WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
210           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
211           WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
212           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
213           WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
214           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
215           _END_MASTER( )
216          ENDIF
217    
218    C--   Transfert the 2D-solution to "etaN" :
219          DO bj=myByLo(myThid),myByHi(myThid)
220           DO bi=myBxLo(myThid),myBxHi(myThid)
221            DO j=1-OLy,sNy+OLy
222             DO i=1-OLx,sNx+OLx
223              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
224             ENDDO
225            ENDDO
226           ENDDO
227          ENDDO
228    
229    #ifdef ALLOW_NONHYDROSTATIC
230          IF ( nonHydrostatic ) THEN
231    
232  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
233  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
234  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
235  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
236             DO j=1,sNy+1
237              DO i=1,sNx+1
238               uf(i,j)=-_recip_dxC(i,j,bi,bj)*
239         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
240               vf(i,j)=-_recip_dyC(i,j,bi,bj)*
241         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
242              ENDDO
243             ENDDO
244    
245    #ifdef ALLOW_OBCS
246             IF (useOBCS) THEN
247              DO i=1,sNx+1
248    C Northern boundary
249              IF (OB_Jn(I,bi,bj).NE.0) THEN
250               vf(I,OB_Jn(I,bi,bj))=0.
251              ENDIF
252    C Southern boundary
253              IF (OB_Js(I,bi,bj).NE.0) THEN
254               vf(I,OB_Js(I,bi,bj)+1)=0.
255              ENDIF
256              ENDDO
257              DO j=1,sNy+1
258    C Eastern boundary
259              IF (OB_Ie(J,bi,bj).NE.0) THEN
260               uf(OB_Ie(J,bi,bj),J)=0.
261              ENDIF
262    C Western boundary
263              IF (OB_Iw(J,bi,bj).NE.0) THEN
264               uf(OB_Iw(J,bi,bj)+1,J)=0.
265              ENDIF
266              ENDDO
267             ENDIF
268    #endif
269    
270             K=1
271             DO j=1,sNy
272              DO i=1,sNx
273               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
274         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
275         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
276         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
277         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
278         &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
279         &          -wVel(i,j,k+1,bi,bj)
280         &        )*_rA(i,j,bi,bj)/deltaTmom
281              ENDDO
282             ENDDO
283             DO K=2,Nr-1
284              DO j=1,sNy
285               DO i=1,sNx
286                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
287         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
288         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
289         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
290         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
291         &       +( wVel(i,j,k  ,bi,bj)
292         &         -wVel(i,j,k+1,bi,bj)
293         &        )*_rA(i,j,bi,bj)/deltaTmom
294    
295               ENDDO
296              ENDDO
297             ENDDO
298             K=Nr
299             DO j=1,sNy
300              DO i=1,sNx
301                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
302         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
303         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
304         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
305         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
306         &       +( wVel(i,j,k  ,bi,bj)
307         &        )*_rA(i,j,bi,bj)/deltaTmom
308    
309              ENDDO
310             ENDDO
311    
312    #ifdef ALLOW_OBCS
313             IF (useOBCS) THEN
314              DO K=1,Nr
315              DO i=1,sNx
316    C Northern boundary
317              IF (OB_Jn(I,bi,bj).NE.0) THEN
318               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
319              ENDIF
320    C Southern boundary
321              IF (OB_Js(I,bi,bj).NE.0) THEN
322               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
323              ENDIF
324              ENDDO
325              DO j=1,sNy
326    C Eastern boundary
327              IF (OB_Ie(J,bi,bj).NE.0) THEN
328               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
329              ENDIF
330    C Western boundary
331              IF (OB_Iw(J,bi,bj).NE.0) THEN
332               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
333              ENDIF
334              ENDDO
335              ENDDO
336             ENDIF
337    #endif
338    
339            ENDDO ! bi
340           ENDDO ! bj
341    
342          firstResidual=0.
343          lastResidual=0.
344          numIters=cg2dMaxIters
345          CALL CG3D(
346         U           cg3d_b,
347         U           phi_nh,
348         O           firstResidual,
349         O           lastResidual,
350         U           numIters,
351         I           myThid )
352          _EXCH_XYZ_R8(phi_nh, myThid )
353    
354          _BEGIN_MASTER( myThid )
355          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
356          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
357          WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
358          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
359          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
360          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
361          _END_MASTER( )
362    
363          ENDIF
364    #endif
365    
366        RETURN        RETURN
367        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22