/[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.16 by jmc, Tue Feb 20 15:08:34 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 "GRID.h"
21  #include "CG2D.h"  #include "CG2D.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        INTEGER i,j,k,bi,bj
37          _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38          _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39    
40  #ifdef ALLOW_CD        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    
59    #ifdef INCLUDE_CD_CODE
60  C--   Save previous solution.  C--   Save previous solution.
61        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
62         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 38  C--   Save previous solution. Line 69  C--   Save previous solution.
69        ENDDO        ENDDO
70  #endif  #endif
71    
72    C--   Add source term arising from w=d/dt (p_s + p_nh)
73          DO bj=myByLo(myThid),myByHi(myThid)
74           DO bi=myBxLo(myThid),myBxHi(myThid)
75    #ifdef ALLOW_NONHYDROSTATIC
76            DO j=1,sNy
77             DO i=1,sNx
78              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
79         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
80         &          -cg2d_x(I,J,bi,bj)
81         &          -cg3d_x(I,J,1,bi,bj)
82         &        )/deltaTMom/deltaTMom
83              cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
84         &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
85         &         -cg2d_x(I,J,bi,bj)
86         &         -cg3d_x(I,J,1,bi,bj)
87         &       )/deltaTMom/deltaTMom
88             ENDDO
89            ENDDO
90    #else
91            DO j=1,sNy
92             DO i=1,sNx
93              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
94         &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
95         &          -cg2d_x(I,J,bi,bj)
96         &        )/deltaTMom/deltaTMom
97             ENDDO
98            ENDDO
99    #endif
100    
101    #ifdef ALLOW_OBCS
102            IF (useOBCS) 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.
132        CALL CG2D(        CALL CG2D(
133         I           cg2d_b,
134         U           cg2d_x,
135       I           myThid )       I           myThid )
136    
137          _EXCH_XY_R8(cg2d_x, myThid )
138    
139    #ifdef ALLOW_NONHYDROSTATIC
140          IF ( nonHydrostatic ) THEN
141    
142  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 ).
143  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
144  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
145  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
146             DO j=1,sNy+1
147              DO i=1,sNx+1
148               uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*
149         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
150               vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*
151         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
152              ENDDO
153             ENDDO
154    
155    #ifdef ALLOW_OBCS
156             IF (useOBCS) THEN
157              DO i=1,sNx+1
158    C Northern boundary
159              IF (OB_Jn(I,bi,bj).NE.0) THEN
160               vf(I,OB_Jn(I,bi,bj))=0.
161              ENDIF
162    C Southern boundary
163              IF (OB_Js(I,bi,bj).NE.0) THEN
164               vf(I,OB_Js(I,bi,bj)+1)=0.
165              ENDIF
166              ENDDO
167              DO j=1,sNy+1
168    C Eastern boundary
169              IF (OB_Ie(J,bi,bj).NE.0) THEN
170               uf(OB_Ie(J,bi,bj),J)=0.
171              ENDIF
172    C Western boundary
173              IF (OB_Iw(J,bi,bj).NE.0) THEN
174               uf(OB_Iw(J,bi,bj)+1,J)=0.
175              ENDIF
176              ENDDO
177             ENDIF
178    #endif
179    
180             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
198               DO i=1,sNx
199                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)
201         &       -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)
203         &       -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
209              ENDDO
210             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 (useOBCS) 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
253           ENDDO ! bj
254    
255           CALL CG3D( myThid )
256           _EXCH_XYZ_R8(cg3d_x, myThid )
257    
258          ENDIF
259    #endif
260    
261        RETURN        RETURN
262        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22