/[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.15 by cnh, Sun Feb 4 14:38:48 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"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "DYNVARS.h"  #include "DYNVARS.h"
26  #include "GRID.h"  #include "GRID.h"
27  #include "CG2D.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     === 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  C     !LOCAL VARIABLES:
52    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
57          _RL tmpFac
58          INTEGER numIters
59          CHARACTER*(MAX_LEN_MBUF) msgBuf
60    CEOP
61    
62    C--   Save previous solution & Initialise Vector solution and source term :
63          DO bj=myByLo(myThid),myByHi(myThid)
64           DO bi=myBxLo(myThid),myBxHi(myThid)
65            DO j=1-OLy,sNy+OLy
66             DO i=1-OLx,sNx+OLx
67    #ifdef INCLUDE_CD_CODE
68              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)
71              cg2d_b(i,j,bi,bj) = 0.
72             ENDDO
73            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
86          ENDDO
87    
 #ifndef DIVG_IN_DYNAMICS  
88        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
89         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
90          DO K=Nr,1,-1          DO K=Nr,1,-1
# Line 52  C     Local variables Line 99  C     Local variables
99           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
100       I       bi,bj,1,sNx,1,sNy,K,       I       bi,bj,1,sNx,1,sNy,K,
101       I       uf,vf,       I       uf,vf,
102         U       cg2d_b,
103       I       myThid)       I       myThid)
104          ENDDO          ENDDO
105         ENDDO         ENDDO
106        ENDDO        ENDDO
 #endif  
   
 #ifdef INCLUDE_CD_CODE  
 C--   Save previous solution.  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
 #endif  
107    
108  C--   Add source term arising from w=d/dt (p_s + p_nh)  C--   Add source term arising from w=d/dt (p_s + p_nh)
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)*horiVertRatio*(             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
116       &          -cg2d_x(I,J,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
117       &          -cg3d_x(I,J,1,bi,bj)       &         *( etaN(i,j,bi,bj)
118       &        )/deltaTMom/deltaTMom       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
119            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)
120       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
121       &         -cg2d_x(I,J,bi,bj)       &         *( etaN(i,j,bi,bj)
122       &         -cg3d_x(I,J,1,bi,bj)       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
123       &       )/deltaTMom/deltaTMom            ENDDO
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)*horiVertRatio*(  
      &          -cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/deltaTMom  
          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 106  C--   Add source term arising from w=d/d Line 149  C--   Add source term arising from w=d/d
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 127  C Western boundary Line 174  C Western boundary
174         ENDDO         ENDDO
175        ENDDO        ENDDO
176    
177    #ifndef DISABLE_DEBUGMODE
178          IF (debugMode) THEN
179           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
180         &                        myThid)
181          ENDIF
182    #endif
183    
184  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
185  C--   gradient solver.  C--   gradient solver.
186  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
187          firstResidual=0.
188          lastResidual=0.
189          numIters=cg2dMaxIters
190        CALL CG2D(        CALL CG2D(
191       I           cg2d_b,       U           cg2d_b,
192       U           cg2d_x,       U           cg2d_x,
193         O           firstResidual,
194         O           lastResidual,
195         U           numIters,
196       I           myThid )       I           myThid )
   
197        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
198    
199    #ifndef DISABLE_DEBUGMODE
200          IF (debugMode) THEN
201           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
202         &                        myThid)
203          ENDIF
204    #endif
205    
206    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
207          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,
208         &                                    myTime-deltaTClock) ) THEN
209           _BEGIN_MASTER( myThid )
210           WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
211           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" :
220          DO bj=myByLo(myThid),myByHi(myThid)
221           DO bi=myBxLo(myThid),myBxHi(myThid)
222            DO j=1-OLy,sNy+OLy
223             DO i=1-OLx,sNx+OLx
224              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
225             ENDDO
226            ENDDO
227           ENDDO
228          ENDDO
229    
230  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
231        IF ( nonHydrostatic ) THEN        IF ( nonHydrostatic ) THEN
232    
# Line 147  C     see CG3D.h for the interface to th Line 236  C     see CG3D.h for the interface to th
236          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
237           DO j=1,sNy+1           DO j=1,sNy+1
238            DO i=1,sNx+1            DO i=1,sNx+1
239             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
240       &         (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))
241             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
242       &         (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))
243            ENDDO            ENDDO
244           ENDDO           ENDDO
# Line 187  C Western boundary Line 276  C Western boundary
276       &       -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)
277       &       +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)
278       &       -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 )
279       &       +(       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
280       &         -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
281       &        )*_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  
282            ENDDO            ENDDO
283           ENDDO           ENDDO
284           DO K=2,Nr-1           DO K=2,Nr-1
# Line 254  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.15  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22