/[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.9 by adcroft, Mon Mar 22 15:54:05 1999 UTC revision 1.25 by adcroft, Fri Jun 29 17:14:49 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 16  C     == Global variables Line 17  C     == Global variables
17  #include "EEPARAMS.h"  #include "EEPARAMS.h"
18  #include "PARAMS.h"  #include "PARAMS.h"
19  #include "DYNVARS.h"  #include "DYNVARS.h"
20  #include "CG2D.h"  #include "GRID.h"
21    #include "SURFACE.h"
22  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
23  #include "CG3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
24  #include "GW.h"  #include "GW.h"
25    #endif
26    #ifdef ALLOW_OBCS
27  #include "OBCS.h"  #include "OBCS.h"
 #include "GRID.h"  
28  #endif  #endif
29    #include "SOLVE_FOR_PRESSURE.h"
30    
31  C     == Routine arguments ==  C     == Routine arguments ==
32  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
33        INTEGER myThid        INTEGER myThid
34  CEndOfInterface  CEndOfInterface
35    
36  C     Local variables  C     == Local variables ==
37        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
 #ifdef ALLOW_NONHYDROSTATIC  
38        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40  #endif        _RL firstResidual,lastResidual
41          INTEGER numIters
42          CHARACTER*(MAX_LEN_MBUF) msgBuf
43    
44  #ifdef INCLUDE_CD_CODE  C--   Save previous solution & Initialise Vector solution and source term :
 C--   Save previous solution.  
45        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
46         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
47          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
48           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
49            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
50              cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
51              cg2d_b(i,j,bi,bj) = 0.
52    #ifdef USE_NATURAL_BCS
53         &     + freeSurfFac*_rA(i,j,bi,bj)*
54         &       EmPmR(I,J,bi,bj)/deltaTMom
55    #endif
56             ENDDO
57            ENDDO
58           ENDDO
59          ENDDO
60    
61          DO bj=myByLo(myThid),myByHi(myThid)
62           DO bi=myBxLo(myThid),myBxHi(myThid)
63            DO K=Nr,1,-1
64             DO j=1,sNy+1
65              DO i=1,sNx+1
66               uf(i,j) = _dyG(i,j,bi,bj)
67         &      *drF(k)*_hFacW(i,j,k,bi,bj)
68               vf(i,j) = _dxG(i,j,bi,bj)
69         &      *drF(k)*_hFacS(i,j,k,bi,bj)
70              ENDDO
71             ENDDO
72             CALL CALC_DIV_GHAT(
73         I       bi,bj,1,sNx,1,sNy,K,
74         I       uf,vf,
75         U       cg2d_b,
76         I       myThid)
77            ENDDO
78           ENDDO
79          ENDDO
80    
81    C--   Add source term arising from w=d/dt (p_s + p_nh)
82          DO bj=myByLo(myThid),myByHi(myThid)
83           DO bi=myBxLo(myThid),myBxHi(myThid)
84    #ifdef ALLOW_NONHYDROSTATIC
85            DO j=1,sNy
86             DO i=1,sNx
87              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
88         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
89         &        *( etaN(i,j,bi,bj)
90         &          +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
91              cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
92         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
93         &        *( etaN(i,j,bi,bj)
94         &          +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
95    C-jmc
96    c    &      -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
97    c    &        *( cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj) )
98    c    &        /deltaTMom/deltaTMom
99    C-jmc
100           ENDDO           ENDDO
101          ENDDO          ENDDO
102    #else
103            DO j=1,sNy
104             DO i=1,sNx
105              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
106         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
107         &        * etaN(i,j,bi,bj)
108             ENDDO
109            ENDDO
110    #endif
111    
112    #ifdef ALLOW_OBCS
113            IF (useOBCS) THEN
114             DO i=1,sNx
115    C Northern boundary
116              IF (OB_Jn(I,bi,bj).NE.0) THEN
117               cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
118              ENDIF
119    C Southern boundary
120              IF (OB_Js(I,bi,bj).NE.0) THEN
121               cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
122              ENDIF
123             ENDDO
124             DO j=1,sNy
125    C Eastern boundary
126              IF (OB_Ie(J,bi,bj).NE.0) THEN
127               cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
128              ENDIF
129    C Western boundary
130              IF (OB_Iw(J,bi,bj).NE.0) THEN
131               cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
132              ENDIF
133             ENDDO
134            ENDIF
135    #endif
136         ENDDO         ENDDO
137        ENDDO        ENDDO
138    
139    #ifndef EXCLUDE_DEBUGMODE
140          IF (debugMode) THEN
141           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
142         &                        myThid)
143          ENDIF
144  #endif  #endif
145    
146  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
147  C--   gradient solver.  C--   gradient solver.
148  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
149          firstResidual=0.
150          lastResidual=0.
151          numIters=cg2dMaxIters
152        CALL CG2D(        CALL CG2D(
153       I           cg2d_b,       U           cg2d_b,
154       U           cg2d_x,       U           cg2d_x,
155         O           firstResidual,
156         O           lastResidual,
157         U           numIters,
158       I           myThid )       I           myThid )
159          _EXCH_XY_R8(cg2d_x, myThid )
160    
161    #ifndef EXCLUDE_DEBUGMODE
162          IF (debugMode) THEN
163           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
164         &                        myThid)
165          ENDIF
166    #endif
167    
168          _BEGIN_MASTER( myThid )
169          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
170          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
171          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
172          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
173          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
174          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
175          _END_MASTER( )
176    
177    C--   Transfert the 2D-solution to "etaN" :
178          DO bj=myByLo(myThid),myByHi(myThid)
179           DO bi=myBxLo(myThid),myBxHi(myThid)
180            DO j=1-OLy,sNy+OLy
181             DO i=1-OLx,sNx+OLx
182              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
183             ENDDO
184            ENDDO
185           ENDDO
186          ENDDO
187    
188  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
189        IF ( nonHydrostatic ) THEN        IF ( nonHydrostatic ) THEN
# Line 66  C     see CG3D.h for the interface to th Line 194  C     see CG3D.h for the interface to th
194          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
195           DO j=1,sNy+1           DO j=1,sNy+1
196            DO i=1,sNx+1            DO i=1,sNx+1
197             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
198       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
199             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
200       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
201            ENDDO            ENDDO
202           ENDDO           ENDDO
203    
204           IF (openBoundaries) THEN  #ifdef ALLOW_OBCS
205             IF (useOBCS) THEN
206            DO i=1,sNx+1            DO i=1,sNx+1
207  C Northern boundary  C Northern boundary
208            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(I,bi,bj).NE.0) THEN
            uf(I,OB_Jn(I,bi,bj))=0.  
209             vf(I,OB_Jn(I,bi,bj))=0.             vf(I,OB_Jn(I,bi,bj))=0.
210            ENDIF            ENDIF
211  C Southern boundary  C Southern boundary
212            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(I,bi,bj).NE.0) THEN
            uf(I,OB_Js(I,bi,bj))=0.  
213             vf(I,OB_Js(I,bi,bj)+1)=0.             vf(I,OB_Js(I,bi,bj)+1)=0.
214            ENDIF            ENDIF
215            ENDDO            ENDDO
# Line 90  C Southern boundary Line 217  C Southern boundary
217  C Eastern boundary  C Eastern boundary
218            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
219             uf(OB_Ie(J,bi,bj),J)=0.             uf(OB_Ie(J,bi,bj),J)=0.
            vf(OB_Ie(J,bi,bj),J)=0.  
220            ENDIF            ENDIF
221  C Western boundary  C Western boundary
222            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
223             uf(OB_Iw(J,bi,bj)+1,J)=0.             uf(OB_Iw(J,bi,bj)+1,J)=0.
            vf(OB_Iw(J,bi,bj),J)=0.  
224            ENDIF            ENDIF
225            ENDDO            ENDDO
226           ENDIF           ENDIF
227    #endif
228    
229           DO K=1,Nr           K=1
230             DO j=1,sNy
231              DO i=1,sNx
232               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
233         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
234         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
235         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
236         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
237         &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
238         &          -wVel(i,j,k+1,bi,bj)
239         &        )*_rA(i,j,bi,bj)/deltaTmom
240              ENDDO
241             ENDDO
242             DO K=2,Nr-1
243            DO j=1,sNy            DO j=1,sNy
244             DO i=1,sNx             DO i=1,sNx
 c           cg3d_x(i,j,k,bi,bj) = 0.  
245              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
246       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
247       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
248       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
249       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
250         &       +( wVel(i,j,k  ,bi,bj)
251         &         -wVel(i,j,k+1,bi,bj)
252         &        )*_rA(i,j,bi,bj)/deltaTmom
253    
254             ENDDO             ENDDO
255            ENDDO            ENDDO
256           ENDDO           ENDDO
257             K=Nr
258             DO j=1,sNy
259              DO i=1,sNx
260                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
261         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
262         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
263         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
264         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
265         &       +( wVel(i,j,k  ,bi,bj)
266         &        )*_rA(i,j,bi,bj)/deltaTmom
267    
268              ENDDO
269             ENDDO
270    
271    #ifdef ALLOW_OBCS
272             IF (useOBCS) THEN
273              DO K=1,Nr
274              DO i=1,sNx
275    C Northern boundary
276              IF (OB_Jn(I,bi,bj).NE.0) THEN
277               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
278              ENDIF
279    C Southern boundary
280              IF (OB_Js(I,bi,bj).NE.0) THEN
281               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
282              ENDIF
283              ENDDO
284              DO j=1,sNy
285    C Eastern boundary
286              IF (OB_Ie(J,bi,bj).NE.0) THEN
287               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
288              ENDIF
289    C Western boundary
290              IF (OB_Iw(J,bi,bj).NE.0) THEN
291               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
292              ENDIF
293              ENDDO
294              ENDDO
295             ENDIF
296    #endif
297    
298          ENDDO ! bi          ENDDO ! bi
299         ENDDO ! bj         ENDDO ! bj
300    
301         CALL CG3D( myThid )        firstResidual=0.
302          lastResidual=0.
303          numIters=cg2dMaxIters
304          CALL CG3D(
305         U           cg3d_b,
306         U           phi_nh,
307         O           firstResidual,
308         O           lastResidual,
309         U           numIters,
310         I           myThid )
311          _EXCH_XYZ_R8(phi_nh, myThid )
312    
313          _BEGIN_MASTER( myThid )
314          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
315          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
316          WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
317          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
318          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
319          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320          _END_MASTER( )
321    
322        ENDIF        ENDIF
323  #endif  #endif
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22