/[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.17 by jmc, Tue Mar 6 16:57:10 2001 UTC revision 1.42 by edhill, Tue Nov 4 18:40:58 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )  C     !ROUTINE: SOLVE_FOR_PRESSURE
9  C     /==========================================================\  C     !INTERFACE:
10  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
11  C     | o Controls inversion of two and/or three-dimensional     |  
12  C     |   elliptic problems for the pressure field.              |  C     !DESCRIPTION: \bv
13  C     \==========================================================/  C     *==========================================================*
14        IMPLICIT NONE  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            
15    C     | o Controls inversion of two and/or three-dimensional      
16    C     |   elliptic problems for the pressure field.              
17    C     *==========================================================*
18    C     \ev
19    
20    C     !USES:
21          IMPLICIT NONE
22  C     == Global variables  C     == Global variables
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "DYNVARS.h"  #include "DYNVARS.h"
27    #ifdef ALLOW_CD_CODE
28    #include "CD_CODE_VARS.h"
29    #endif
30  #include "GRID.h"  #include "GRID.h"
31  #include "SURFACE.h"  #include "SURFACE.h"
32    #include "FFIELDS.h"
33  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
34  #include "CG3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
35  #include "GW.h"  #include "GW.h"
36  #endif  #endif
37  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
38  #include "OBCS.h"  #include "OBCS.h"
39  #endif  #endif
40    #include "SOLVE_FOR_PRESSURE.h"
41    
42    C     === Functions ====
43          LOGICAL  DIFFERENT_MULTIPLE
44          EXTERNAL DIFFERENT_MULTIPLE
45    
46    C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime - Current time in simulation
49    C     myIter - Current iteration number in simulation
50    C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
51          _RL myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
54    
55  C     Local variables  C     !LOCAL VARIABLES:
56  C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector  C     == Local variables ==
 C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  
57        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
58        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL firstResidual,lastResidual
61        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL tmpFac
62          INTEGER numIters
63          CHARACTER*(MAX_LEN_MBUF) msgBuf
64    CEOP
65    
66  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
67        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
68         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
69          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
70           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
71  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
72            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
73  #endif  #endif
74            cg2d_x(i,j,bi,bj) = etaN(i,j,bi,bj)            cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
75            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
 #ifdef USE_NATURAL_BCS  
      &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*  
      &       EmPmR(I,J,bi,bj)/deltaTMom  
 #endif  
76           ENDDO           ENDDO
77          ENDDO          ENDDO
78            IF (useRealFreshWaterFlux) THEN
79             tmpFac = freeSurfFac*convertEmP2rUnit
80             IF (exactConserv)
81         &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
82             DO j=1,sNy
83              DO i=1,sNx
84               cg2d_b(i,j,bi,bj) =
85         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
86              ENDDO
87             ENDDO
88            ENDIF
89         ENDDO         ENDDO
90        ENDDO        ENDDO
91    
# Line 84  C--   Add source term arising from w=d/d Line 113  C--   Add source term arising from w=d/d
113        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
114         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
115  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
116          DO j=1,sNy          IF ( nonHydrostatic ) THEN
117           DO i=1,sNx           DO j=1,sNy
118            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
119       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
120       &          -cg2d_x(I,J,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
121       &          -cg3d_x(I,J,1,bi,bj)       &         *( etaN(i,j,bi,bj)
122       &        )/deltaTMom/deltaTMom       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
123            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)             cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
124       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
125       &         -cg2d_x(I,J,bi,bj)       &         *( etaN(i,j,bi,bj)
126       &         -cg3d_x(I,J,1,bi,bj)       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
127       &       )/deltaTMom/deltaTMom            ENDDO
128           ENDDO           ENDDO
129          ENDDO          ELSEIF ( exactConserv ) THEN
130  #else  #else
131          DO j=1,sNy          IF ( exactConserv ) THEN
132           DO i=1,sNx  #endif /* ALLOW_NONHYDROSTATIC */
133            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)           DO j=1,sNy
134       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(            DO i=1,sNx
135       &          -cg2d_x(I,J,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
136       &        )/deltaTMom/deltaTMom       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
137         &         * etaH(i,j,bi,bj)
138              ENDDO
139           ENDDO           ENDDO
140          ENDDO          ELSE
141  #endif           DO j=1,sNy
142              DO i=1,sNx
143               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
144         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
145         &         * etaN(i,j,bi,bj)
146              ENDDO
147             ENDDO
148            ENDIF
149    
150  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
151          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 115  C--   Add source term arising from w=d/d Line 153  C--   Add source term arising from w=d/d
153  C Northern boundary  C Northern boundary
154            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(I,bi,bj).NE.0) THEN
155             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
156               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
157            ENDIF            ENDIF
158  C Southern boundary  C Southern boundary
159            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(I,bi,bj).NE.0) THEN
160             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
161               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
162            ENDIF            ENDIF
163           ENDDO           ENDDO
164           DO j=1,sNy           DO j=1,sNy
165  C Eastern boundary  C Eastern boundary
166            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
167             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
168               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
169            ENDIF            ENDIF
170  C Western boundary  C Western boundary
171            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
172             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
173               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
174            ENDIF            ENDIF
175           ENDDO           ENDDO
176          ENDIF          ENDIF
# Line 136  C Western boundary Line 178  C Western boundary
178         ENDDO         ENDDO
179        ENDDO        ENDDO
180    
181    #ifdef ALLOW_DEBUG
182          IF ( debugLevel .GE. debLevB ) THEN
183           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
184         &                        myThid)
185          ENDIF
186    #endif
187    
188  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
189  C--   gradient solver.  C--   gradient solver.
190  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
191          firstResidual=0.
192          lastResidual=0.
193          numIters=cg2dMaxIters
194        CALL CG2D(        CALL CG2D(
195       I           cg2d_b,       U           cg2d_b,
196       U           cg2d_x,       U           cg2d_x,
197         O           firstResidual,
198         O           lastResidual,
199         U           numIters,
200       I           myThid )       I           myThid )
   
201        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
202    
203    #ifdef ALLOW_DEBUG
204          IF ( debugLevel .GE. debLevB ) THEN
205           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
206         &                        myThid)
207          ENDIF
208    #endif
209    
210    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
211          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
212         &                                    myTime-deltaTClock) ) THEN
213           IF ( debugLevel .GE. debLevA ) THEN
214            _BEGIN_MASTER( myThid )
215            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
216            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
217            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
218            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
219            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
220            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
221            _END_MASTER( )
222           ENDIF
223          ENDIF
224    
225  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
226        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
227         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
228          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
229           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
230            etaN(i,j,bi,bj) = cg2d_x(i,j,bi,bj)            etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
231           ENDDO           ENDDO
232          ENDDO          ENDDO
233         ENDDO         ENDDO
# Line 167  C     see CG3D.h for the interface to th Line 242  C     see CG3D.h for the interface to th
242          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
243           DO j=1,sNy+1           DO j=1,sNy+1
244            DO i=1,sNx+1            DO i=1,sNx+1
245             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
246       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
247             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
248       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
249            ENDDO            ENDDO
250           ENDDO           ENDDO
# Line 207  C Western boundary Line 282  C Western boundary
282       &       -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)
283       &       +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)
284       &       -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 )
285       &       +(       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
286       &         -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
287       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
      &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(  
      &          +cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/deltaTMom  
288            ENDDO            ENDDO
289           ENDDO           ENDDO
290           DO K=2,Nr-1           DO K=2,Nr-1
# Line 274  C Western boundary Line 346  C Western boundary
346          ENDDO ! bi          ENDDO ! bi
347         ENDDO ! bj         ENDDO ! bj
348    
349         CALL CG3D( myThid )        firstResidual=0.
350         _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual=0.
351          numIters=cg2dMaxIters
352          CALL CG3D(
353         U           cg3d_b,
354         U           phi_nh,
355         O           firstResidual,
356         O           lastResidual,
357         U           numIters,
358         I           myThid )
359          _EXCH_XYZ_R8(phi_nh, myThid )
360    
361          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
362         &                                    myTime-deltaTClock) ) THEN
363           IF ( debugLevel .GE. debLevA ) THEN
364            _BEGIN_MASTER( myThid )
365            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
366            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
367            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
368            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
369            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
370            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
371            _END_MASTER( )
372           ENDIF
373          ENDIF
374    
375        ENDIF        ENDIF
376  #endif  #endif
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22