/[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.12 by adcroft, Tue Mar 14 17:47:26 2000 UTC
# Line 16  C     == Global variables Line 16  C     == Global variables
16  #include "EEPARAMS.h"  #include "EEPARAMS.h"
17  #include "PARAMS.h"  #include "PARAMS.h"
18  #include "DYNVARS.h"  #include "DYNVARS.h"
19    #include "GRID.h"
20  #include "CG2D.h"  #include "CG2D.h"
21  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
22  #include "CG3D.h"  #include "CG3D.h"
23  #include "GW.h"  #include "GW.h"
24    #endif
25  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
26  #include "OBCS.h"  #include "OBCS.h"
27  #endif  #endif
 #include "GRID.h"  
 #endif  
28    
29  C     == Routine arguments ==  C     == Routine arguments ==
30  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
# Line 33  CEndOfInterface Line 33  CEndOfInterface
33    
34  C     Local variables  C     Local variables
35        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
 #ifdef ALLOW_NONHYDROSTATIC  
36        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38    
39    #ifndef DIVG_IN_DYNAMICS
40          DO bj=myByLo(myThid),myByHi(myThid)
41           DO bi=myBxLo(myThid),myBxHi(myThid)
42            DO K=Nr,1,-1
43             DO j=1,sNy+1
44              DO i=1,sNx+1
45               uf(i,j) = _dyG(i,j,bi,bj)
46         &      *drF(k)*_hFacW(i,j,k,bi,bj)
47               vf(i,j) = _dxG(i,j,bi,bj)
48         &      *drF(k)*_hFacS(i,j,k,bi,bj)
49              ENDDO
50             ENDDO
51             CALL CALC_DIV_GHAT(
52         I       bi,bj,1,sNx,1,sNy,K,
53         I       uf,vf,
54         I       myThid)
55            ENDDO
56           ENDDO
57          ENDDO
58  #endif  #endif
59    
60  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
# Line 51  C--   Save previous solution. Line 70  C--   Save previous solution.
70        ENDDO        ENDDO
71  #endif  #endif
72    
73    C--   Add source term arising from w=d/dt (p_s + p_nh)
74          DO bj=myByLo(myThid),myByHi(myThid)
75           DO bi=myBxLo(myThid),myBxHi(myThid)
76            K=1
77            DO j=1,sNy
78             DO i=1,sNx
79              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
80         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
81         &          -cg2d_x(I,J,bi,bj)
82    #ifdef ALLOW_NONHYDROSTATIC
83         &          -cg3d_x(I,J,K,bi,bj)
84    #endif
85         &        )/deltaTMom/deltaTMom
86             ENDDO
87            ENDDO
88    #ifdef ALLOW_NONHYDROSTATIC
89            K=1
90            DO j=1,sNy
91             DO i=1,sNx
92              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
93         &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
94         &         -cg2d_x(I,J,bi,bj)
95         &         -cg3d_x(I,J,K,bi,bj)
96         &       )/deltaTMom/deltaTMom
97             ENDDO
98            ENDDO
99    #endif
100    
101    #ifdef ALLOW_OBCS
102            IF (openBoundaries) THEN
103             DO i=1,sNx
104    C Northern boundary
105              IF (OB_Jn(I,bi,bj).NE.0) THEN
106               cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
107              ENDIF
108    C Southern boundary
109              IF (OB_Js(I,bi,bj).NE.0) THEN
110               cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
111              ENDIF
112             ENDDO
113             DO j=1,sNy
114    C Eastern boundary
115              IF (OB_Ie(J,bi,bj).NE.0) THEN
116               cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
117              ENDIF
118    C Western boundary
119              IF (OB_Iw(J,bi,bj).NE.0) THEN
120               cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
121              ENDIF
122             ENDDO
123            ENDIF
124    #endif
125           ENDDO
126          ENDDO
127    
128    
129  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
130  C--   gradient solver.  C--   gradient solver.
131  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
# Line 77  C     see CG3D.h for the interface to th Line 152  C     see CG3D.h for the interface to th
152            ENDDO            ENDDO
153           ENDDO           ENDDO
154    
155    #ifdef ALLOW_OBCS
156           IF (openBoundaries) THEN           IF (openBoundaries) THEN
157            DO i=1,sNx+1            DO i=1,sNx+1
158  C Northern boundary  C Northern boundary
159            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.  
160             vf(I,OB_Jn(I,bi,bj))=0.             vf(I,OB_Jn(I,bi,bj))=0.
161            ENDIF            ENDIF
162  C Southern boundary  C Southern boundary
163            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.  
164             vf(I,OB_Js(I,bi,bj)+1)=0.             vf(I,OB_Js(I,bi,bj)+1)=0.
165            ENDIF            ENDIF
166            ENDDO            ENDDO
# Line 94  C Southern boundary Line 168  C Southern boundary
168  C Eastern boundary  C Eastern boundary
169            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
170             uf(OB_Ie(J,bi,bj),J)=0.             uf(OB_Ie(J,bi,bj),J)=0.
            vf(OB_Ie(J,bi,bj),J)=0.  
171            ENDIF            ENDIF
172  C Western boundary  C Western boundary
173            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
174             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.  
175            ENDIF            ENDIF
176            ENDDO            ENDDO
177           ENDIF           ENDIF
178    #endif
179    
180           DO K=1,Nr           K=1
181             DO j=1,sNy
182              DO i=1,sNx
183               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
184         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
185         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
186         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
187         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
188         &       +(
189         &         -wVel(i,j,k+1,bi,bj)
190         &        )*_rA(i,j,bi,bj)/deltaTmom
191         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
192         &          +cg2d_x(I,J,bi,bj)
193         &        )/deltaTMom/deltaTMom
194              ENDDO
195             ENDDO
196             DO K=2,Nr-1
197            DO j=1,sNy            DO j=1,sNy
198             DO i=1,sNx             DO i=1,sNx
 c           cg3d_x(i,j,k,bi,bj) = 0.  
199              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)
200       &       +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)
201       &       -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)
202       &       +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)
203       &       -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 )
204         &       +( wVel(i,j,k  ,bi,bj)
205         &         -wVel(i,j,k+1,bi,bj)
206         &        )*_rA(i,j,bi,bj)/deltaTmom
207    
208             ENDDO             ENDDO
209            ENDDO            ENDDO
210           ENDDO           ENDDO
211             K=Nr
212             DO j=1,sNy
213              DO i=1,sNx
214                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
215         &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
216         &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
217         &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
218         &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
219         &       +( wVel(i,j,k  ,bi,bj)
220         &        )*_rA(i,j,bi,bj)/deltaTmom
221    
222              ENDDO
223             ENDDO
224    
225    #ifdef ALLOW_OBCS
226             IF (openBoundaries) THEN
227              DO K=1,Nr
228              DO i=1,sNx
229    C Northern boundary
230              IF (OB_Jn(I,bi,bj).NE.0) THEN
231               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
232              ENDIF
233    C Southern boundary
234              IF (OB_Js(I,bi,bj).NE.0) THEN
235               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
236              ENDIF
237              ENDDO
238              DO j=1,sNy
239    C Eastern boundary
240              IF (OB_Ie(J,bi,bj).NE.0) THEN
241               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
242              ENDIF
243    C Western boundary
244              IF (OB_Iw(J,bi,bj).NE.0) THEN
245               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
246              ENDIF
247              ENDDO
248              ENDDO
249             ENDIF
250    #endif
251    
252          ENDDO ! bi          ENDDO ! bi
253         ENDDO ! bj         ENDDO ! bj
254    
255         CALL CG3D( myThid )         CALL CG3D( myThid )
   
256         _EXCH_XYZ_R8(cg3d_x, myThid )         _EXCH_XYZ_R8(cg3d_x, myThid )
257    
258        ENDIF        ENDIF
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22