/[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.11 by adcroft, Mon May 24 15:42:23 1999 UTC revision 1.17 by jmc, Tue Mar 6 16:57:10 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 "CG3D.h"
24  #include "GW.h"  #include "GW.h"
25    #endif
26  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
27  #include "OBCS.h"  #include "OBCS.h"
28  #endif  #endif
 #include "GRID.h"  
 #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 32  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    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        INTEGER i,j,k,bi,bj
 #ifdef ALLOW_NONHYDROSTATIC  
39        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41  #endif        _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    
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)  #ifdef INCLUDE_CD_CODE
50              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
51    #endif
52              cg2d_x(i,j,bi,bj) = etaN(i,j,bi,bj)
53              cg2d_b(i,j,bi,bj) = 0.
54    #ifdef USE_NATURAL_BCS
55         &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*
56         &       EmPmR(I,J,bi,bj)/deltaTMom
57    #endif
58             ENDDO
59            ENDDO
60           ENDDO
61          ENDDO
62    
63          DO bj=myByLo(myThid),myByHi(myThid)
64           DO bi=myBxLo(myThid),myBxHi(myThid)
65            DO K=Nr,1,-1
66             DO j=1,sNy+1
67              DO i=1,sNx+1
68               uf(i,j) = _dyG(i,j,bi,bj)
69         &      *drF(k)*_hFacW(i,j,k,bi,bj)
70               vf(i,j) = _dxG(i,j,bi,bj)
71         &      *drF(k)*_hFacS(i,j,k,bi,bj)
72              ENDDO
73           ENDDO           ENDDO
74             CALL CALC_DIV_GHAT(
75         I       bi,bj,1,sNx,1,sNy,K,
76         I       uf,vf,
77         U       cg2d_b,
78         I       myThid)
79          ENDDO          ENDDO
80         ENDDO         ENDDO
81        ENDDO        ENDDO
82    
83    C--   Add source term arising from w=d/dt (p_s + p_nh)
84          DO bj=myByLo(myThid),myByHi(myThid)
85           DO bi=myBxLo(myThid),myBxHi(myThid)
86    #ifdef ALLOW_NONHYDROSTATIC
87            DO j=1,sNy
88             DO i=1,sNx
89              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
90         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
91         &          -cg2d_x(I,J,bi,bj)
92         &          -cg3d_x(I,J,1,bi,bj)
93         &        )/deltaTMom/deltaTMom
94              cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
95         &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
96         &         -cg2d_x(I,J,bi,bj)
97         &         -cg3d_x(I,J,1,bi,bj)
98         &       )/deltaTMom/deltaTMom
99             ENDDO
100            ENDDO
101    #else
102            DO j=1,sNy
103             DO i=1,sNx
104              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
105         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
106         &          -cg2d_x(I,J,bi,bj)
107         &        )/deltaTMom/deltaTMom
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  #endif
136           ENDDO
137          ENDDO
138    
139    
140  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
141  C--   gradient solver.  C--   gradient solver.
142  C     see CG2D.h for the interface to this routine.  C     see CG2D_INTERNAL.h for the interface to this routine.
143        CALL CG2D(        CALL CG2D(
144       I           cg2d_b,       I           cg2d_b,
145       U           cg2d_x,       U           cg2d_x,
# Line 61  C     see CG2D.h for the interface to th Line 147  C     see CG2D.h for the interface to th
147    
148        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
149    
150    C--   Transfert the 2D-solution to "etaN" :
151          DO bj=myByLo(myThid),myByHi(myThid)
152           DO bi=myBxLo(myThid),myBxHi(myThid)
153            DO j=1-OLy,sNy+OLy
154             DO i=1-OLx,sNx+OLx
155              etaN(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
156             ENDDO
157            ENDDO
158           ENDDO
159          ENDDO
160    
161  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
162        IF ( nonHydrostatic ) THEN        IF ( nonHydrostatic ) THEN
163    
# Line 77  C     see CG3D.h for the interface to th Line 174  C     see CG3D.h for the interface to th
174            ENDDO            ENDDO
175           ENDDO           ENDDO
176    
177           IF (openBoundaries) THEN  #ifdef ALLOW_OBCS
178             IF (useOBCS) THEN
179            DO i=1,sNx+1            DO i=1,sNx+1
180  C Northern boundary  C Northern boundary
181            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.  
182             vf(I,OB_Jn(I,bi,bj))=0.             vf(I,OB_Jn(I,bi,bj))=0.
183            ENDIF            ENDIF
184  C Southern boundary  C Southern boundary
185            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.  
186             vf(I,OB_Js(I,bi,bj)+1)=0.             vf(I,OB_Js(I,bi,bj)+1)=0.
187            ENDIF            ENDIF
188            ENDDO            ENDDO
# Line 94  C Southern boundary Line 190  C Southern boundary
190  C Eastern boundary  C Eastern boundary
191            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
192             uf(OB_Ie(J,bi,bj),J)=0.             uf(OB_Ie(J,bi,bj),J)=0.
            vf(OB_Ie(J,bi,bj),J)=0.  
193            ENDIF            ENDIF
194  C Western boundary  C Western boundary
195            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
196             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.  
197            ENDIF            ENDIF
198            ENDDO            ENDDO
199           ENDIF           ENDIF
200    #endif
201    
202           DO K=1,Nr           K=1
203             DO j=1,sNy
204              DO i=1,sNx
205               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
206         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
207         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
208         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
209         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
210         &       +(
211         &         -wVel(i,j,k+1,bi,bj)
212         &        )*_rA(i,j,bi,bj)/deltaTmom
213         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
214         &          +cg2d_x(I,J,bi,bj)
215         &        )/deltaTMom/deltaTMom
216              ENDDO
217             ENDDO
218             DO K=2,Nr-1
219            DO j=1,sNy            DO j=1,sNy
220             DO i=1,sNx             DO i=1,sNx
 c           cg3d_x(i,j,k,bi,bj) = 0.  
221              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)
222       &       +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)
223       &       -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)
224       &       +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)
225       &       -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 )
226         &       +( wVel(i,j,k  ,bi,bj)
227         &         -wVel(i,j,k+1,bi,bj)
228         &        )*_rA(i,j,bi,bj)/deltaTmom
229    
230             ENDDO             ENDDO
231            ENDDO            ENDDO
232           ENDDO           ENDDO
233             K=Nr
234             DO j=1,sNy
235              DO i=1,sNx
236                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
237         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
238         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
239         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
240         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
241         &       +( wVel(i,j,k  ,bi,bj)
242         &        )*_rA(i,j,bi,bj)/deltaTmom
243    
244              ENDDO
245             ENDDO
246    
247    #ifdef ALLOW_OBCS
248             IF (useOBCS) THEN
249              DO K=1,Nr
250              DO i=1,sNx
251    C Northern boundary
252              IF (OB_Jn(I,bi,bj).NE.0) THEN
253               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
254              ENDIF
255    C Southern boundary
256              IF (OB_Js(I,bi,bj).NE.0) THEN
257               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
258              ENDIF
259              ENDDO
260              DO j=1,sNy
261    C Eastern boundary
262              IF (OB_Ie(J,bi,bj).NE.0) THEN
263               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
264              ENDIF
265    C Western boundary
266              IF (OB_Iw(J,bi,bj).NE.0) THEN
267               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
268              ENDIF
269              ENDDO
270              ENDDO
271             ENDIF
272    #endif
273    
274          ENDDO ! bi          ENDDO ! bi
275         ENDDO ! bj         ENDDO ! bj
276    
277         CALL CG3D( myThid )         CALL CG3D( myThid )
   
278         _EXCH_XYZ_R8(cg3d_x, myThid )         _EXCH_XYZ_R8(cg3d_x, myThid )
279    
280        ENDIF        ENDIF
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22