/[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.2 by cnh, Fri Apr 24 02:05:42 1998 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_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
16    #include "SIZE.h"
17    #include "EEPARAMS.h"
18    #include "PARAMS.h"
19    #include "DYNVARS.h"
20    #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
32        INTEGER myThid        INTEGER myThid
33  CEndOfInterface  CEndOfInterface
34    
35    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
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    
44    C--   Save previous solution & Initialise Vector solution and source term :
45          DO bj=myByLo(myThid),myByHi(myThid)
46           DO bi=myBxLo(myThid),myBxHi(myThid)
47            DO j=1-OLy,sNy+OLy
48             DO i=1-OLx,sNx+OLx
49    #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
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
80           ENDDO
81          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
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,
145         U           cg2d_x,
146       I           myThid )       I           myThid )
147    
148          _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
162          IF ( nonHydrostatic ) THEN
163    
164  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 ).
165  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
166  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
167  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
168             DO j=1,sNy+1
169              DO i=1,sNx+1
170               uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*
171         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
172               vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*
173         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
174              ENDDO
175             ENDDO
176    
177    #ifdef ALLOW_OBCS
178             IF (useOBCS) THEN
179              DO i=1,sNx+1
180    C Northern boundary
181              IF (OB_Jn(I,bi,bj).NE.0) THEN
182               vf(I,OB_Jn(I,bi,bj))=0.
183              ENDIF
184    C Southern boundary
185              IF (OB_Js(I,bi,bj).NE.0) THEN
186               vf(I,OB_Js(I,bi,bj)+1)=0.
187              ENDIF
188              ENDDO
189              DO j=1,sNy+1
190    C Eastern boundary
191              IF (OB_Ie(J,bi,bj).NE.0) THEN
192               uf(OB_Ie(J,bi,bj),J)=0.
193              ENDIF
194    C Western boundary
195              IF (OB_Iw(J,bi,bj).NE.0) THEN
196               uf(OB_Iw(J,bi,bj)+1,J)=0.
197              ENDIF
198              ENDDO
199             ENDIF
200    #endif
201    
202             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
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         &       +( wVel(i,j,k  ,bi,bj)
227         &         -wVel(i,j,k+1,bi,bj)
228         &        )*_rA(i,j,bi,bj)/deltaTmom
229    
230               ENDDO
231              ENDDO
232             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
275           ENDDO ! bj
276    
277           CALL CG3D( myThid )
278           _EXCH_XYZ_R8(cg3d_x, myThid )
279    
280          ENDIF
281    #endif
282    
283        RETURN        RETURN
284        END        END
 C $Id$  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22