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