/[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.19 by adcroft, Fri Mar 9 20:45:09 2001 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  CStartOfInterface
7        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )
# Line 9  C     | SUBROUTINE SOLVE_FOR_PRESSURE Line 10  C     | SUBROUTINE SOLVE_FOR_PRESSURE
10  C     | o Controls inversion of two and/or three-dimensional     |  C     | o Controls inversion of two and/or three-dimensional     |
11  C     |   elliptic problems for the pressure field.              |  C     |   elliptic problems for the pressure field.              |
12  C     \==========================================================/  C     \==========================================================/
13          IMPLICIT NONE
14    
15  C     == Global variables  C     == Global variables
16  #include "SIZE.h"  #include "SIZE.h"
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
23    #include "CG3D.h"
24    #include "GW.h"
25    #endif
26    #ifdef ALLOW_OBCS
27    #include "OBCS.h"
28    #endif
29    
30  C     == Routine arguments ==  C     == Routine arguments ==
31  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
# Line 23  C     myThid - Number of this instance o Line 33  C     myThid - Number of this instance o
33  CEndOfInterface  CEndOfInterface
34    
35  C     Local variables  C     Local variables
36        INTEGER i,j,bi,bj  C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector
37    C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
38          INTEGER i,j,k,bi,bj
39          _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40          _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41          _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42          _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43          _RL tolerance,residual
44          INTEGER numIters
45    
46  #ifdef ALLOW_CD  C--   Save previous solution & Initialise Vector solution and source term :
 C--   Save previous solution.  
47        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
48         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
49          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
50           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
51            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  #ifdef INCLUDE_CD_CODE
52              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
53    #endif
54              cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
55              cg2d_b(i,j,bi,bj) = 0.
56    #ifdef USE_NATURAL_BCS
57         &     + freeSurfFac*_rA(i,j,bi,bj)*
58         &       EmPmR(I,J,bi,bj)/deltaTMom
59    #endif
60           ENDDO           ENDDO
61          ENDDO          ENDDO
62         ENDDO         ENDDO
63        ENDDO        ENDDO
64    
65          DO bj=myByLo(myThid),myByHi(myThid)
66           DO bi=myBxLo(myThid),myBxHi(myThid)
67            DO K=Nr,1,-1
68             DO j=1,sNy+1
69              DO i=1,sNx+1
70               uf(i,j) = _dyG(i,j,bi,bj)
71         &      *drF(k)*_hFacW(i,j,k,bi,bj)
72               vf(i,j) = _dxG(i,j,bi,bj)
73         &      *drF(k)*_hFacS(i,j,k,bi,bj)
74              ENDDO
75             ENDDO
76             CALL CALC_DIV_GHAT(
77         I       bi,bj,1,sNx,1,sNy,K,
78         I       uf,vf,
79         U       cg2d_b,
80         I       myThid)
81            ENDDO
82           ENDDO
83          ENDDO
84    
85    C--   Add source term arising from w=d/dt (p_s + p_nh)
86          DO bj=myByLo(myThid),myByHi(myThid)
87           DO bi=myBxLo(myThid),myBxHi(myThid)
88    #ifdef ALLOW_NONHYDROSTATIC
89            DO j=1,sNy
90             DO i=1,sNx
91              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
92         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
93         &        *( etaN(i,j,bi,bj)
94         &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
95              cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
96         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
97         &        *( etaN(i,j,bi,bj)
98         &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
99    C-jmc
100    c    &      -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
101    c    &        *( cg2d_x(i,j,bi,bj) + cg3d_x(i,j,1,bi,bj) )
102    c    &        /deltaTMom/deltaTMom
103    C-jmc
104             ENDDO
105            ENDDO
106    #else
107            DO j=1,sNy
108             DO i=1,sNx
109              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
110         &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
111         &        * etaN(i,j,bi,bj)
112             ENDDO
113            ENDDO
114  #endif  #endif
115    
116    #ifdef ALLOW_OBCS
117            IF (useOBCS) THEN
118             DO i=1,sNx
119    C Northern boundary
120              IF (OB_Jn(I,bi,bj).NE.0) THEN
121               cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
122              ENDIF
123    C Southern boundary
124              IF (OB_Js(I,bi,bj).NE.0) THEN
125               cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
126              ENDIF
127             ENDDO
128             DO j=1,sNy
129    C Eastern boundary
130              IF (OB_Ie(J,bi,bj).NE.0) THEN
131               cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
132              ENDIF
133    C Western boundary
134              IF (OB_Iw(J,bi,bj).NE.0) THEN
135               cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
136              ENDIF
137             ENDDO
138            ENDIF
139    #endif
140           ENDDO
141          ENDDO
142    
143    
144  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
145  C--   gradient solver.  C--   gradient solver.
146  C     see CG2D.h for the interface to this routine.  C     see CG2D_INTERNAL.h for the interface to this routine.
147          tolerance=cg2dTargetResidual
148          residual=0.
149          numIters=cg2dMaxIters
150        CALL CG2D(        CALL CG2D(
151         I           cg2d_b,
152         U           cg2d_x,
153         U           tolerance,
154         U           residual,
155         U           numIters,
156       I           myThid )       I           myThid )
157          _EXCH_XY_R8(cg2d_x, myThid )
158    
159          _BEGIN_MASTER( myThid )
160           WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
161         &                         0, tolerance
162           WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
163         &                         numIters, residual
164          _END_MASTER( )
165    
166    C--   Transfert the 2D-solution to "etaN" :
167          DO bj=myByLo(myThid),myByHi(myThid)
168           DO bi=myBxLo(myThid),myBxHi(myThid)
169            DO j=1-OLy,sNy+OLy
170             DO i=1-OLx,sNx+OLx
171              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
172             ENDDO
173            ENDDO
174           ENDDO
175          ENDDO
176    
177    #ifdef ALLOW_NONHYDROSTATIC
178          IF ( nonHydrostatic ) THEN
179    
180  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 ).
181  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
182  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
183  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
184             DO j=1,sNy+1
185              DO i=1,sNx+1
186               uf(i,j)=-_recip_dxC(i,j,bi,bj)*
187         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
188               vf(i,j)=-_recip_dyC(i,j,bi,bj)*
189         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
190              ENDDO
191             ENDDO
192    
193    #ifdef ALLOW_OBCS
194             IF (useOBCS) THEN
195              DO i=1,sNx+1
196    C Northern boundary
197              IF (OB_Jn(I,bi,bj).NE.0) THEN
198               vf(I,OB_Jn(I,bi,bj))=0.
199              ENDIF
200    C Southern boundary
201              IF (OB_Js(I,bi,bj).NE.0) THEN
202               vf(I,OB_Js(I,bi,bj)+1)=0.
203              ENDIF
204              ENDDO
205              DO j=1,sNy+1
206    C Eastern boundary
207              IF (OB_Ie(J,bi,bj).NE.0) THEN
208               uf(OB_Ie(J,bi,bj),J)=0.
209              ENDIF
210    C Western boundary
211              IF (OB_Iw(J,bi,bj).NE.0) THEN
212               uf(OB_Iw(J,bi,bj)+1,J)=0.
213              ENDIF
214              ENDDO
215             ENDIF
216    #endif
217    
218             K=1
219             DO j=1,sNy
220              DO i=1,sNx
221               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
222         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
223         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
224         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
225         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
226         &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
227         &          -wVel(i,j,k+1,bi,bj)
228         &        )*_rA(i,j,bi,bj)/deltaTmom
229              ENDDO
230             ENDDO
231             DO K=2,Nr-1
232              DO j=1,sNy
233               DO i=1,sNx
234                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
235         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
236         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
237         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
238         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
239         &       +( wVel(i,j,k  ,bi,bj)
240         &         -wVel(i,j,k+1,bi,bj)
241         &        )*_rA(i,j,bi,bj)/deltaTmom
242    
243               ENDDO
244              ENDDO
245             ENDDO
246             K=Nr
247             DO j=1,sNy
248              DO i=1,sNx
249                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
250         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
251         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
252         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
253         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
254         &       +( wVel(i,j,k  ,bi,bj)
255         &        )*_rA(i,j,bi,bj)/deltaTmom
256    
257              ENDDO
258             ENDDO
259    
260    #ifdef ALLOW_OBCS
261             IF (useOBCS) THEN
262              DO K=1,Nr
263              DO i=1,sNx
264    C Northern boundary
265              IF (OB_Jn(I,bi,bj).NE.0) THEN
266               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
267              ENDIF
268    C Southern boundary
269              IF (OB_Js(I,bi,bj).NE.0) THEN
270               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
271              ENDIF
272              ENDDO
273              DO j=1,sNy
274    C Eastern boundary
275              IF (OB_Ie(J,bi,bj).NE.0) THEN
276               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
277              ENDIF
278    C Western boundary
279              IF (OB_Iw(J,bi,bj).NE.0) THEN
280               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
281              ENDIF
282              ENDDO
283              ENDDO
284             ENDIF
285    #endif
286    
287            ENDDO ! bi
288           ENDDO ! bj
289    
290           CALL CG3D( myThid )
291           _EXCH_XYZ_R8(cg3d_x, myThid )
292    
293          ENDIF
294    #endif
295    
296        RETURN        RETURN
297        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22