/[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.38 by heimbach, Tue Jul 8 15:00:26 2003 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*convertEmP2rUnit
76             IF (exactConserv)
77         &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
78             DO j=1,sNy
79              DO i=1,sNx
80               cg2d_b(i,j,bi,bj) =
81         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
82              ENDDO
83             ENDDO
84            ENDIF
85         ENDDO         ENDDO
86        ENDDO        ENDDO
87    
88          DO bj=myByLo(myThid),myByHi(myThid)
89           DO bi=myBxLo(myThid),myBxHi(myThid)
90            DO K=Nr,1,-1
91             DO j=1,sNy+1
92              DO i=1,sNx+1
93               uf(i,j) = _dyG(i,j,bi,bj)
94         &      *drF(k)*_hFacW(i,j,k,bi,bj)
95               vf(i,j) = _dxG(i,j,bi,bj)
96         &      *drF(k)*_hFacS(i,j,k,bi,bj)
97              ENDDO
98             ENDDO
99             CALL CALC_DIV_GHAT(
100         I       bi,bj,1,sNx,1,sNy,K,
101         I       uf,vf,
102         U       cg2d_b,
103         I       myThid)
104            ENDDO
105           ENDDO
106          ENDDO
107    
108    C--   Add source term arising from w=d/dt (p_s + p_nh)
109          DO bj=myByLo(myThid),myByHi(myThid)
110           DO bi=myBxLo(myThid),myBxHi(myThid)
111    #ifdef ALLOW_NONHYDROSTATIC
112            IF ( nonHydrostatic ) THEN
113             DO j=1,sNy
114              DO i=1,sNx
115               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
116         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
117         &         *( etaN(i,j,bi,bj)
118         &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
119               cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
120         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
121         &         *( etaN(i,j,bi,bj)
122         &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
123              ENDDO
124             ENDDO
125            ELSEIF ( exactConserv ) THEN
126    #else
127            IF ( exactConserv ) THEN
128    #endif
129             DO j=1,sNy
130              DO i=1,sNx
131               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
132         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
133         &         * etaH(i,j,bi,bj)
134              ENDDO
135             ENDDO
136            ELSE
137             DO j=1,sNy
138              DO i=1,sNx
139               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
140         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
141         &         * etaN(i,j,bi,bj)
142              ENDDO
143             ENDDO
144            ENDIF
145    
146    #ifdef ALLOW_OBCS
147            IF (useOBCS) THEN
148             DO i=1,sNx
149    C Northern boundary
150              IF (OB_Jn(I,bi,bj).NE.0) THEN
151               cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
152               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
153              ENDIF
154    C Southern boundary
155              IF (OB_Js(I,bi,bj).NE.0) THEN
156               cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
157               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
158              ENDIF
159             ENDDO
160             DO j=1,sNy
161    C Eastern boundary
162              IF (OB_Ie(J,bi,bj).NE.0) THEN
163               cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
164               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
165              ENDIF
166    C Western boundary
167              IF (OB_Iw(J,bi,bj).NE.0) THEN
168               cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
169               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
170              ENDIF
171             ENDDO
172            ENDIF
173    #endif
174           ENDDO
175          ENDDO
176    
177    #ifndef DISABLE_DEBUGMODE
178          IF ( debugLevel .GE. debLevB ) THEN
179           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
180         &                        myThid)
181          ENDIF
182  #endif  #endif
183    
184  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
185  C--   gradient solver.  C--   gradient solver.
186  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
187          firstResidual=0.
188          lastResidual=0.
189          numIters=cg2dMaxIters
190        CALL CG2D(        CALL CG2D(
191         U           cg2d_b,
192         U           cg2d_x,
193         O           firstResidual,
194         O           lastResidual,
195         U           numIters,
196       I           myThid )       I           myThid )
197          _EXCH_XY_R8(cg2d_x, myThid )
198    
199    #ifndef DISABLE_DEBUGMODE
200          IF ( debugLevel .GE. debLevB ) THEN
201           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
202         &                        myThid)
203          ENDIF
204    #endif
205    
206    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
207          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
208         &                                    myTime-deltaTClock) ) THEN
209           IF ( debugLevel .GE. debLevA ) THEN
210            _BEGIN_MASTER( myThid )
211            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
212            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
213            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
214            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
215            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
216            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
217            _END_MASTER( )
218           ENDIF
219          ENDIF
220    
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              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
227             ENDDO
228            ENDDO
229           ENDDO
230          ENDDO
231    
232    #ifdef ALLOW_NONHYDROSTATIC
233          IF ( nonHydrostatic ) THEN
234    
235  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 ).
236  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
237  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
238  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
239             DO j=1,sNy+1
240              DO i=1,sNx+1
241               uf(i,j)=-_recip_dxC(i,j,bi,bj)*
242         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
243               vf(i,j)=-_recip_dyC(i,j,bi,bj)*
244         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
245              ENDDO
246             ENDDO
247    
248    #ifdef ALLOW_OBCS
249             IF (useOBCS) THEN
250              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    #endif
272    
273             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         &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
282         &          -wVel(i,j,k+1,bi,bj)
283         &        )*_rA(i,j,bi,bj)/deltaTmom
284              ENDDO
285             ENDDO
286             DO K=2,Nr-1
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         &         -wVel(i,j,k+1,bi,bj)
296         &        )*_rA(i,j,bi,bj)/deltaTmom
297    
298               ENDDO
299              ENDDO
300             ENDDO
301             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             IF (useOBCS) THEN
317              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    
342            ENDDO ! bi
343           ENDDO ! bj
344    
345          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          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
358         &                                    myTime-deltaTClock) ) THEN
359           IF ( debugLevel .GE. debLevA ) THEN
360            _BEGIN_MASTER( myThid )
361            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
362            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
363            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
364            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
365            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
366            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
367            _END_MASTER( )
368           ENDIF
369          ENDIF
370    
371          ENDIF
372    #endif
373    
374        RETURN        RETURN
375        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.38
  ViewVC Help
Powered by ViewVC 1.1.22