/[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.4 by cnh, Mon Jun 8 21:43:02 1998 UTC revision 1.77 by jmc, Wed Jun 8 01:21:14 2011 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #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 )
 C     | o Controls inversion of two and/or three-dimensional     |  
 C     |   elliptic problems for the pressure field.              |  
 C     \==========================================================/  
11    
12    C     !DESCRIPTION: \bv
13    C     *==========================================================*
14    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 "GRID.h"
27    #include "SURFACE.h"
28    #include "FFIELDS.h"
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "CG2D.h"  #include "SOLVE_FOR_PRESSURE.h"
31    #ifdef ALLOW_NONHYDROSTATIC
32    #include "SOLVE_FOR_PRESSURE3D.h"
33    #include "NH_VARS.h"
34    #endif
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37    #endif
38    
39    C     === Functions ====
40          LOGICAL  DIFFERENT_MULTIPLE
41          EXTERNAL DIFFERENT_MULTIPLE
42    
43    C     !INPUT/OUTPUT PARAMETERS:
44  C     == Routine arguments ==  C     == Routine arguments ==
45  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime :: Current time in simulation
46    C     myIter :: Current iteration number in simulation
47    C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
48          _RL myTime
49          INTEGER myIter
50        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
51    
52  C     Local variables  C     !LOCAL VARIABLES:
53        INTEGER i,j,bi,bj  C     == Local variables ==
54          INTEGER i,j,k,bi,bj
55          INTEGER ks
56          INTEGER numIters
57          _RL firstResidual,lastResidual
58          _RL tmpFac
59          _RL sumEmP, tileEmP(nSx,nSy)
60          LOGICAL putPmEinXvector
61          INTEGER ioUnit
62          CHARACTER*10 sufx
63          CHARACTER*(MAX_LEN_MBUF) msgBuf
64    #ifdef ALLOW_NONHYDROSTATIC
65          LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
66    #else
67          _RL     cg3d_b(1)
68    #endif
69    CEOP
70    
71    #ifdef ALLOW_NONHYDROSTATIC
72            zeroPsNH = .FALSE.
73    c       zeroPsNH = use3Dsolver .AND. exactConserv
74    c    &                         .AND. select_rStar.EQ.0
75            zeroMeanPnh = .FALSE.
76    c       zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
77    c       oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
78    c    &                                .AND. .NOT.zeroPsNH
79            oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
80    #else
81            cg3d_b(1) = 0.
82    #endif
83    
84    C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
85    C anelastic (always Z-coordinate):
86    C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
87    C        (this reduces the number of lines of code to modify)
88    C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
89    C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
90    C       => 2 factors cancel in elliptic eq. for Phi_s ,
91    C       but 1rst factor(a) remains in RHS cg2d_b.
92    
93    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
94    C     instead of simply etaN ; This can speed-up the solver convergence in
95    C     the case where |Global_mean_PmE| is large.
96          putPmEinXvector = .FALSE.
97    c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
98    
99          IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
100            _BEGIN_MASTER( myThid )
101            ioUnit = standardMessageUnit
102            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
103         &       ' putPmEinXvector =', putPmEinXvector
104            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
105    #ifdef ALLOW_NONHYDROSTATIC
106            WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
107         &       ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
108            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
109            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
110         &       ' oldFreeSurfTerm =', oldFreeSurfTerm
111            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
112    #endif
113            _END_MASTER( myThid )
114          ENDIF
115    
116  #ifdef ALLOW_CD  C--   Save previous solution & Initialise Vector solution and source term :
117  C--   Save previous solution.        sumEmP = 0.
118        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
119         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
120          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
121           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
122            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  #ifdef ALLOW_CD_CODE
123              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
124    #endif
125              cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
126              cg2d_b(i,j,bi,bj) = 0.
127             ENDDO
128            ENDDO
129            IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
130             tmpFac = freeSurfFac*mass2rUnit
131             IF (exactConserv)
132         &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
133             DO j=1,sNy
134              DO i=1,sNx
135               cg2d_b(i,j,bi,bj) =
136         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
137              ENDDO
138           ENDDO           ENDDO
139            ENDIF
140            IF ( putPmEinXvector ) THEN
141             tileEmP(bi,bj) = 0.
142             DO j=1,sNy
143              DO i=1,sNx
144                tileEmP(bi,bj) = tileEmP(bi,bj)
145         &                     + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
146         &                                    *maskInC(i,j,bi,bj)
147              ENDDO
148             ENDDO
149            ENDIF
150           ENDDO
151          ENDDO
152          IF ( putPmEinXvector ) THEN
153            CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
154          ENDIF
155    
156          DO bj=myByLo(myThid),myByHi(myThid)
157           DO bi=myBxLo(myThid),myBxHi(myThid)
158            IF ( putPmEinXvector ) THEN
159              tmpFac = 0.
160              IF (globalArea.GT.0.) tmpFac =
161         &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
162              DO j=1,sNy
163               DO i=1,sNx
164                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
165         &                        - tmpFac*Bo_surf(i,j,bi,bj)
166               ENDDO
167              ENDDO
168            ENDIF
169    C- RHS: similar to the divergence of the vertically integrated mass transport:
170    C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
171            DO k=Nr,1,-1
172             CALL CALC_DIV_GHAT(
173         I                       bi,bj,k,
174         U                       cg2d_b, cg3d_b,
175         I                       myThid )
176          ENDDO          ENDDO
177         ENDDO         ENDDO
178        ENDDO        ENDDO
179    
180          DO bj=myByLo(myThid),myByHi(myThid)
181           DO bi=myBxLo(myThid),myBxHi(myThid)
182    #ifdef ALLOW_NONHYDROSTATIC
183            IF ( oldFreeSurfTerm ) THEN
184    C--   Add source term arising from w=d/dt (p_s + p_nh)
185             DO j=1,sNy
186              DO i=1,sNx
187               ks = kSurfC(i,j,bi,bj)
188               IF ( ks.LE.Nr ) THEN
189                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
190         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
191         &         /deltaTMom/deltaTfreesurf
192         &         *( etaN(i,j,bi,bj)
193         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
194                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
195         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
196         &         /deltaTMom/deltaTfreesurf
197         &         *( etaN(i,j,bi,bj)
198         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
199               ENDIF
200              ENDDO
201             ENDDO
202            ELSEIF ( exactConserv ) THEN
203    #else
204    C--   Add source term arising from w=d/dt (p_s)
205            IF ( exactConserv ) THEN
206    #endif /* ALLOW_NONHYDROSTATIC */
207             DO j=1,sNy
208              DO i=1,sNx
209               ks = kSurfC(i,j,bi,bj)
210               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
211         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
212         &         /deltaTMom/deltaTfreesurf
213         &         * etaH(i,j,bi,bj)
214              ENDDO
215             ENDDO
216            ELSE
217             DO j=1,sNy
218              DO i=1,sNx
219               ks = kSurfC(i,j,bi,bj)
220               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
221         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
222         &         /deltaTMom/deltaTfreesurf
223         &         * etaN(i,j,bi,bj)
224              ENDDO
225             ENDDO
226            ENDIF
227    
228    #ifdef ALLOW_OBCS
229    C- Note: solver matrix is trivial outside OB region (main diagonal only)
230    C     => no real need to reset RHS (=cg2d_b) & cg2d_x, except that:
231    C    a) normalisation is fct of Max(RHS), which can be large ouside OB region
232    C      (would be different if we were solving for increment of eta/g
233    C       instead of directly for eta/g).
234    C       => need to reset RHS to ensure that interior solution does not depend
235    C       on ouside OB region.
236    C    b) provide directly the trivial solution cg2d_x == 0 for outside OB region
237    C      (=> no residual => no effect on solver convergence and interior solution)
238            IF (useOBCS) THEN
239             DO j=1,sNy
240              DO i=1,sNx
241               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*maskInC(i,j,bi,bj)
242               cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*maskInC(i,j,bi,bj)
243             ENDDO
244             ENDDO
245            ENDIF
246    #endif /* ALLOW_OBCS */
247    C-    end bi,bj loops
248           ENDDO
249          ENDDO
250    
251    #ifdef ALLOW_DEBUG
252          IF ( debugLevel .GE. debLevD ) THEN
253           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
254         &                        myThid)
255          ENDIF
256  #endif  #endif
257          IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
258           WRITE(sufx,'(I10.10)') myIter
259           CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
260          ENDIF
261    
262  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
263  C--   gradient solver.  C--   gradient solver.
264  C     see CG2D.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
265        CALL CG2D(        firstResidual=0.
266          lastResidual=0.
267          numIters=cg2dMaxIters
268    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
269    #ifdef ALLOW_CG2D_NSA
270    C--   Call the not-self-adjoint version of cg2d
271          CALL CG2D_NSA(
272         U           cg2d_b,
273         U           cg2d_x,
274         O           firstResidual,
275         O           lastResidual,
276         U           numIters,
277       I           myThid )       I           myThid )
278    #else /* not ALLOW_CG2D_NSA = default */
279    #ifdef ALLOW_SRCG
280          IF ( useSRCGSolver ) THEN
281    C--   Call the single reduce CG solver
282           CALL CG2D_SR(
283         U           cg2d_b,
284         U           cg2d_x,
285         O           firstResidual,
286         O           lastResidual,
287         U           numIters,
288         I           myThid )
289          ELSE
290    #else
291          IF (.TRUE.) THEN
292    C--   Call the default CG solver
293    #endif /* ALLOW_SRCG */
294           CALL CG2D(
295         U           cg2d_b,
296         U           cg2d_x,
297         O           firstResidual,
298         O           lastResidual,
299         U           numIters,
300         I           myThid )
301          ENDIF
302    #endif /* ALLOW_CG2D_NSA */
303          _EXCH_XY_RL( cg2d_x, myThid )
304    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
305    
306    #ifdef ALLOW_DEBUG
307          IF ( debugLevel .GE. debLevD ) THEN
308           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
309         &                        myThid)
310          ENDIF
311    #endif
312    
313    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
314          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
315         &   ) THEN
316           IF ( debugLevel .GE. debLevA ) THEN
317            _BEGIN_MASTER( myThid )
318            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
319            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
321            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
322            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
323            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
324            _END_MASTER( myThid )
325           ENDIF
326          ENDIF
327    
328    C--   Transfert the 2D-solution to "etaN" :
329          DO bj=myByLo(myThid),myByHi(myThid)
330           DO bi=myBxLo(myThid),myBxHi(myThid)
331            DO j=1-OLy,sNy+OLy
332             DO i=1-OLx,sNx+OLx
333              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
334             ENDDO
335            ENDDO
336           ENDDO
337          ENDDO
338    
339    #ifdef ALLOW_NONHYDROSTATIC
340          IF ( use3Dsolver ) THEN
341           IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
342            WRITE(sufx,'(I10.10)') myIter
343            CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
344           ENDIF
345    
346  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).  C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
347  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
348  C*P*  CALL CG3D(  
349  C*P* I           myThid )  C--   Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
350    C                             2) Update or Add free-surface contribution
351    C                             3) increment in horiz velocity due to new cg2d_x
352    C                             4) add vertical velocity contribution.
353           CALL PRE_CG3D(
354         I                oldFreeSurfTerm,
355         I                cg2d_x,
356         U                cg3d_b,
357         I                myTime, myIter, myThid )
358    
359    #ifdef ALLOW_DEBUG
360           IF ( debugLevel .GE. debLevD ) THEN
361            CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
362         &                         myThid)
363           ENDIF
364    #endif
365           IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
366            WRITE(sufx,'(I10.10)') myIter
367            CALL WRITE_FLD_XYZ_RL('cg3d_b.',sufx, cg3d_b, myIter,myThid )
368           ENDIF
369    
370           firstResidual=0.
371           lastResidual=0.
372           numIters=cg3dMaxIters
373           CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
374           CALL CG3D(
375         U            cg3d_b,
376         U            phi_nh,
377         O            firstResidual,
378         O            lastResidual,
379         U            numIters,
380         I            myIter, myThid )
381           _EXCH_XYZ_RL( phi_nh, myThid )
382           CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
383    
384           IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
385         &    ) THEN
386            IF ( debugLevel .GE. debLevA ) THEN
387             _BEGIN_MASTER( myThid )
388             WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
389             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
390             WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
391             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
392             WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
393             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
394             _END_MASTER( myThid )
395            ENDIF
396           ENDIF
397    
398    C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
399    C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
400           IF ( nonHydrostatic .AND. exactConserv ) THEN
401            IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
402             WRITE(sufx,'(I10.10)') myIter
403             CALL WRITE_FLD_XYZ_RL('cg3d_x.',sufx, phi_nh, myIter,myThid )
404            ENDIF
405            CALL POST_CG3D(
406         I                  zeroPsNH, zeroMeanPnh,
407         I                  myTime, myIter, myThid )
408           ENDIF
409    
410          ENDIF
411    #endif /* ALLOW_NONHYDROSTATIC */
412    
413    #ifdef ALLOW_SHOWFLOPS
414          CALL SHOWFLOPS_INSOLVE( myThid)
415    #endif
416    
417        RETURN        RETURN
418        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.77
  ViewVC Help
Powered by ViewVC 1.1.22