/[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.23 by adcroft, Wed Jun 6 14:55:45 2001 UTC revision 1.37 by mlosch, Tue May 13 13:30:05 2003 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
# Line 28  C     == Global variables Line 35  C     == Global variables
35  #endif  #endif
36  #include "SOLVE_FOR_PRESSURE.h"  #include "SOLVE_FOR_PRESSURE.h"
37    
38    C     === Functions ====
39          LOGICAL  DIFFERENT_MULTIPLE
40          EXTERNAL DIFFERENT_MULTIPLE
41    
42    C     !INPUT/OUTPUT PARAMETERS:
43  C     == Routine arguments ==  C     == Routine arguments ==
44  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime - Current time in simulation
45    C     myIter - Current iteration number in simulation
46    C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE
47          _RL myTime
48          INTEGER myIter
49        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
50    
51    C     !LOCAL VARIABLES:
52  C     == Local variables ==  C     == Local variables ==
53        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
54        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
57          _RL tmpFac
58        INTEGER numIters        INTEGER numIters
59          CHARACTER*(MAX_LEN_MBUF) msgBuf
60    CEOP
61    
62  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
63        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
64         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
65          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
66           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
67    #ifdef INCLUDE_CD_CODE
68            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
69    #endif
70            cg2d_x(i,j,bi,bj) = Bo_surf(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)
71            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
 #ifdef USE_NATURAL_BCS  
      &     + freeSurfFac*_rA(i,j,bi,bj)*  
      &       EmPmR(I,J,bi,bj)/deltaTMom  
 #endif  
72           ENDDO           ENDDO
73          ENDDO          ENDDO
74            IF (useRealFreshWaterFlux) THEN
75             tmpFac = freeSurfFac*convertEmP2rUnit
76             IF (exactConserv)
77         &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
78             DO j=1,sNy
79              DO i=1,sNx
80               cg2d_b(i,j,bi,bj) =
81         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
82              ENDDO
83             ENDDO
84            ENDIF
85         ENDDO         ENDDO
86        ENDDO        ENDDO
87    
# Line 81  C--   Add source term arising from w=d/d Line 109  C--   Add source term arising from w=d/d
109        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
110         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
111  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
112          DO j=1,sNy          IF ( nonHydrostatic ) THEN
113           DO i=1,sNx           DO j=1,sNy
114            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
115       &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
116       &        *( etaN(i,j,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
117       &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )       &         *( etaN(i,j,bi,bj)
118            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
119       &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom             cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
120       &        *( etaN(i,j,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
121       &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )       &         *( etaN(i,j,bi,bj)
122  C-jmc       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
123  c    &      -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)            ENDDO
 c    &        *( cg2d_x(i,j,bi,bj) + cg3d_x(i,j,1,bi,bj) )  
 c    &        /deltaTMom/deltaTMom  
 C-jmc  
124           ENDDO           ENDDO
125          ENDDO          ELSEIF ( exactConserv ) THEN
126  #else  #else
127          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)/deltaTMom/deltaTMom  
      &        * etaN(i,j,bi,bj)  
          ENDDO  
         ENDDO  
128  #endif  #endif
129             DO j=1,sNy
130              DO i=1,sNx
131               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
132         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
133         &         * etaH(i,j,bi,bj)
134              ENDDO
135             ENDDO
136            ELSE
137             DO j=1,sNy
138              DO i=1,sNx
139               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
140         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
141         &         * etaN(i,j,bi,bj)
142              ENDDO
143             ENDDO
144            ENDIF
145    
146  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
147          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 114  C-jmc Line 149  C-jmc
149  C Northern boundary  C Northern boundary
150            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(I,bi,bj).NE.0) THEN
151             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
152               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
153            ENDIF            ENDIF
154  C Southern boundary  C Southern boundary
155            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(I,bi,bj).NE.0) THEN
156             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
157               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
158            ENDIF            ENDIF
159           ENDDO           ENDDO
160           DO j=1,sNy           DO j=1,sNy
161  C Eastern boundary  C Eastern boundary
162            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
163             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
164               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
165            ENDIF            ENDIF
166  C Western boundary  C Western boundary
167            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
168             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
169               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
170            ENDIF            ENDIF
171           ENDDO           ENDDO
172          ENDIF          ENDIF
# Line 135  C Western boundary Line 174  C Western boundary
174         ENDDO         ENDDO
175        ENDDO        ENDDO
176    
177  #ifndef EXCLUDE_DEBUGMODE  #ifndef DISABLE_DEBUGMODE
178          IF (debugMode) THEN
179         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
180       &                        myThid)       &                        myThid)
181          ENDIF
182  #endif  #endif
183    
184  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
# Line 155  C     see CG2D.h for the interface to th Line 196  C     see CG2D.h for the interface to th
196       I           myThid )       I           myThid )
197        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
198    
199  #ifndef EXCLUDE_DEBUGMODE  #ifndef DISABLE_DEBUGMODE
200          IF (debugMode) THEN
201         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
202       &                        myThid)       &                        myThid)
203          ENDIF
204  #endif  #endif
205    
206        _BEGIN_MASTER( myThid )  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
207         WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
208       &                               0, firstResidual       &                                    myTime-deltaTClock) ) THEN
209         WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',         _BEGIN_MASTER( myThid )
210       &                         numIters, lastResidual         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
211        _END_MASTER( )         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
212           WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
213           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
214           WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
215           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
216           _END_MASTER( )
217          ENDIF
218    
219  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
220        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 291  C Western boundary Line 340  C Western boundary
340          ENDDO ! bi          ENDDO ! bi
341         ENDDO ! bj         ENDDO ! bj
342    
343         CALL CG3D( myThid )        firstResidual=0.
344         _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual=0.
345          numIters=cg2dMaxIters
346          CALL CG3D(
347         U           cg3d_b,
348         U           phi_nh,
349         O           firstResidual,
350         O           lastResidual,
351         U           numIters,
352         I           myThid )
353          _EXCH_XYZ_R8(phi_nh, myThid )
354    
355          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
356         &                                    myTime-deltaTClock) ) THEN
357           _BEGIN_MASTER( myThid )
358           WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
359           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
360           WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
361           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
362           WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
363           CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
364           _END_MASTER( )
365          ENDIF
366    
367        ENDIF        ENDIF
368  #endif  #endif
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22