/[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.28 by jmc, Fri Feb 8 22:13:39 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        INTEGER myThid  C     myIter - Current iteration number in simulation
42  CEndOfInterface  C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
43          _RL myTime
44  C     Local variables        INTEGER myIter
45  C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector        INTEGER myThid
46  C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  
47    C     !LOCAL VARIABLES:
48    C     == Local variables ==
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)        INTEGER numIters
54          CHARACTER*(MAX_LEN_MBUF) msgBuf
55    CEOP
56    
57  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
58        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 49  C--   Save previous solution & Initialis Line 62  C--   Save previous solution & Initialis
62  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
63            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
64  #endif  #endif
65            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)
66            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
67  #ifdef USE_NATURAL_BCS  #ifdef USE_NATURAL_BCS
68       &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*       &     + freeSurfFac*_rA(i,j,bi,bj)*
69       &       EmPmR(I,J,bi,bj)/deltaTMom       &       EmPmR(I,J,bi,bj)/deltaTMom
70  #endif  #endif
71           ENDDO           ENDDO
# Line 84  C--   Add source term arising from w=d/d Line 97  C--   Add source term arising from w=d/d
97        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
98         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
99  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
100          DO j=1,sNy          IF ( nonHydrostatic ) THEN
101           DO i=1,sNx           DO j=1,sNy
102            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
103       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
104       &          -cg2d_x(I,J,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
105       &          -cg3d_x(I,J,1,bi,bj)       &         *( etaN(i,j,bi,bj)
106       &        )/deltaTMom/deltaTMom       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
107            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)
108       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
109       &         -cg2d_x(I,J,bi,bj)       &         *( etaN(i,j,bi,bj)
110       &         -cg3d_x(I,J,1,bi,bj)       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
111       &       )/deltaTMom/deltaTMom            ENDDO
112           ENDDO           ENDDO
113          ENDDO          ELSEIF ( exactConserv ) THEN
114  #else  #else
115          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  
116  #endif  #endif
117             DO j=1,sNy
118              DO i=1,sNx
119               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
120         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
121         &         * etaH(i,j,bi,bj)
122              ENDDO
123             ENDDO
124            ELSE
125             DO j=1,sNy
126              DO i=1,sNx
127               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
128         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
129         &         * etaN(i,j,bi,bj)
130              ENDDO
131             ENDDO
132            ENDIF
133    
134  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
135          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 136  C Western boundary Line 158  C Western boundary
158         ENDDO         ENDDO
159        ENDDO        ENDDO
160    
161    #ifndef EXCLUDE_DEBUGMODE
162          IF (debugMode) THEN
163           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
164         &                        myThid)
165          ENDIF
166    #endif
167    
168  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
169  C--   gradient solver.  C--   gradient solver.
170  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
171          firstResidual=0.
172          lastResidual=0.
173          numIters=cg2dMaxIters
174        CALL CG2D(        CALL CG2D(
175       I           cg2d_b,       U           cg2d_b,
176       U           cg2d_x,       U           cg2d_x,
177         O           firstResidual,
178         O           lastResidual,
179         U           numIters,
180       I           myThid )       I           myThid )
   
181        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
182    
183    #ifndef EXCLUDE_DEBUGMODE
184          IF (debugMode) THEN
185           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
186         &                        myThid)
187          ENDIF
188    #endif
189    
190          _BEGIN_MASTER( myThid )
191          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
192          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
193          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
194          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
195          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
196          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
197          _END_MASTER( )
198    
199  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
200        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
201         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
202          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
203           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
204            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)
205           ENDDO           ENDDO
206          ENDDO          ENDDO
207         ENDDO         ENDDO
# Line 167  C     see CG3D.h for the interface to th Line 216  C     see CG3D.h for the interface to th
216          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
217           DO j=1,sNy+1           DO j=1,sNy+1
218            DO i=1,sNx+1            DO i=1,sNx+1
219             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
220       &         (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))
221             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
222       &         (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))
223            ENDDO            ENDDO
224           ENDDO           ENDDO
# Line 207  C Western boundary Line 256  C Western boundary
256       &       -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)
257       &       +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)
258       &       -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 )
259       &       +(       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
260       &         -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
261       &        )*_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  
262            ENDDO            ENDDO
263           ENDDO           ENDDO
264           DO K=2,Nr-1           DO K=2,Nr-1
# Line 274  C Western boundary Line 320  C Western boundary
320          ENDDO ! bi          ENDDO ! bi
321         ENDDO ! bj         ENDDO ! bj
322    
323         CALL CG3D( myThid )        firstResidual=0.
324         _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual=0.
325          numIters=cg2dMaxIters
326          CALL CG3D(
327         U           cg3d_b,
328         U           phi_nh,
329         O           firstResidual,
330         O           lastResidual,
331         U           numIters,
332         I           myThid )
333          _EXCH_XYZ_R8(phi_nh, myThid )
334    
335          _BEGIN_MASTER( myThid )
336          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
337          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
338          WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
339          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
340          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
341          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
342          _END_MASTER( )
343    
344        ENDIF        ENDIF
345  #endif  #endif
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22