/[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.29 by jmc, Sun Feb 10 00:39:22 2002 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )  C     !ROUTINE: SOLVE_FOR_PRESSURE
8  C     /==========================================================\  C     !INTERFACE:
9  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
10  C     | o Controls inversion of two and/or three-dimensional     |  
11  C     |   elliptic problems for the pressure field.              |  C     !DESCRIPTION: \bv
12  C     \==========================================================/  C     *==========================================================*
13        IMPLICIT NONE  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            
14    C     | o Controls inversion of two and/or three-dimensional      
15    C     |   elliptic problems for the pressure field.              
16    C     *==========================================================*
17    C     \ev
18    
19    C     !USES:
20          IMPLICIT NONE
21  C     == Global variables  C     == Global variables
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 19  C     == Global variables Line 25  C     == Global variables
25  #include "DYNVARS.h"  #include "DYNVARS.h"
26  #include "GRID.h"  #include "GRID.h"
27  #include "SURFACE.h"  #include "SURFACE.h"
28    #include "FFIELDS.h"
29  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
30  #include "CG3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
31  #include "GW.h"  #include "GW.h"
32  #endif  #endif
33  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
34  #include "OBCS.h"  #include "OBCS.h"
35  #endif  #endif
36    #include "SOLVE_FOR_PRESSURE.h"
37    
38    C     !INPUT/OUTPUT PARAMETERS:
39  C     == Routine arguments ==  C     == Routine arguments ==
40  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime - Current time in simulation
41    C     myIter - Current iteration number in simulation
42    C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
43          _RL myTime
44          INTEGER myIter
45        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
46    
47  C     Local variables  C     !LOCAL VARIABLES:
48  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  
49        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
50        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL firstResidual,lastResidual
53        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL tmpFac
54          INTEGER numIters
55          CHARACTER*(MAX_LEN_MBUF) msgBuf
56    CEOP
57    
58  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
59        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 49  C--   Save previous solution & Initialis Line 63  C--   Save previous solution & Initialis
63  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
64            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
65  #endif  #endif
66            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)
67            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  
68           ENDDO           ENDDO
69          ENDDO          ENDDO
70            IF (useRealFreshWaterFlux) THEN
71             tmpFac = freeSurfFac
72             IF (exactConserv) tmpFac = freeSurfFac*implicDiv2DFlow
73             DO j=1,sNy
74              DO i=1,sNx
75               cg2d_b(i,j,bi,bj) =
76         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
77              ENDDO
78             ENDDO
79            ENDIF
80         ENDDO         ENDDO
81        ENDDO        ENDDO
82    
# Line 84  C--   Add source term arising from w=d/d Line 104  C--   Add source term arising from w=d/d
104        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
105         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
106  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
107          DO j=1,sNy          IF ( nonHydrostatic ) THEN
108           DO i=1,sNx           DO j=1,sNy
109            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
110       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
111       &          -cg2d_x(I,J,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
112       &          -cg3d_x(I,J,1,bi,bj)       &         *( etaN(i,j,bi,bj)
113       &        )/deltaTMom/deltaTMom       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
114            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)
115       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
116       &         -cg2d_x(I,J,bi,bj)       &         *( etaN(i,j,bi,bj)
117       &         -cg3d_x(I,J,1,bi,bj)       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
118       &       )/deltaTMom/deltaTMom            ENDDO
119           ENDDO           ENDDO
120          ENDDO          ELSEIF ( exactConserv ) THEN
121  #else  #else
122          DO j=1,sNy          IF ( exactConserv ) THEN
          DO i=1,sNx  
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)  
      &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(  
      &          -cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/deltaTMom  
          ENDDO  
         ENDDO  
123  #endif  #endif
124             DO j=1,sNy
125              DO i=1,sNx
126               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
127         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
128         &         * etaH(i,j,bi,bj)
129              ENDDO
130             ENDDO
131            ELSE
132             DO j=1,sNy
133              DO i=1,sNx
134               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
135         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
136         &         * etaN(i,j,bi,bj)
137              ENDDO
138             ENDDO
139            ENDIF
140    
141  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
142          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 136  C Western boundary Line 165  C Western boundary
165         ENDDO         ENDDO
166        ENDDO        ENDDO
167    
168    #ifndef EXCLUDE_DEBUGMODE
169          IF (debugMode) THEN
170           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
171         &                        myThid)
172          ENDIF
173    #endif
174    
175  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
176  C--   gradient solver.  C--   gradient solver.
177  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
178          firstResidual=0.
179          lastResidual=0.
180          numIters=cg2dMaxIters
181        CALL CG2D(        CALL CG2D(
182       I           cg2d_b,       U           cg2d_b,
183       U           cg2d_x,       U           cg2d_x,
184         O           firstResidual,
185         O           lastResidual,
186         U           numIters,
187       I           myThid )       I           myThid )
   
188        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
189    
190    #ifndef EXCLUDE_DEBUGMODE
191          IF (debugMode) THEN
192           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
193         &                        myThid)
194          ENDIF
195    #endif
196    
197          _BEGIN_MASTER( myThid )
198          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
199          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
200          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
201          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
202          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
203          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
204          _END_MASTER( )
205    
206  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
207        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
208         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
209          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
210           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
211            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)
212           ENDDO           ENDDO
213          ENDDO          ENDDO
214         ENDDO         ENDDO
# Line 167  C     see CG3D.h for the interface to th Line 223  C     see CG3D.h for the interface to th
223          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
224           DO j=1,sNy+1           DO j=1,sNy+1
225            DO i=1,sNx+1            DO i=1,sNx+1
226             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
227       &         (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))
228             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
229       &         (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))
230            ENDDO            ENDDO
231           ENDDO           ENDDO
# Line 207  C Western boundary Line 263  C Western boundary
263       &       -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)
264       &       +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)
265       &       -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 )
266       &       +(       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
267       &         -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
268       &        )*_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  
269            ENDDO            ENDDO
270           ENDDO           ENDDO
271           DO K=2,Nr-1           DO K=2,Nr-1
# Line 274  C Western boundary Line 327  C Western boundary
327          ENDDO ! bi          ENDDO ! bi
328         ENDDO ! bj         ENDDO ! bj
329    
330         CALL CG3D( myThid )        firstResidual=0.
331         _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual=0.
332          numIters=cg2dMaxIters
333          CALL CG3D(
334         U           cg3d_b,
335         U           phi_nh,
336         O           firstResidual,
337         O           lastResidual,
338         U           numIters,
339         I           myThid )
340          _EXCH_XYZ_R8(phi_nh, myThid )
341    
342          _BEGIN_MASTER( myThid )
343          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
344          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
345          WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
346          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
347          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
348          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
349          _END_MASTER( )
350    
351        ENDIF        ENDIF
352  #endif  #endif
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22