/[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.37 by mlosch, Tue May 13 13:30:05 2003 UTC revision 1.76 by jmc, Wed May 18 23:41:26 2011 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: SOLVE_FOR_PRESSURE  C     !ROUTINE: SOLVE_FOR_PRESSURE
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)        SUBROUTINE SOLVE_FOR_PRESSURE( myTime, myIter, myThid )
11    
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
14  C     | SUBROUTINE SOLVE_FOR_PRESSURE                              C     | SUBROUTINE SOLVE_FOR_PRESSURE
15  C     | o Controls inversion of two and/or three-dimensional        C     | o Controls inversion of two and/or three-dimensional
16  C     |   elliptic problems for the pressure field.                C     |   elliptic problems for the pressure field.
17  C     *==========================================================*  C     *==========================================================*
18  C     \ev  C     \ev
19    
# Line 22  C     == Global variables Line 23  C     == Global variables
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
26  #include "GRID.h"  #include "GRID.h"
27  #include "SURFACE.h"  #include "SURFACE.h"
28  #include "FFIELDS.h"  #include "FFIELDS.h"
29    #include "DYNVARS.h"
30    #include "SOLVE_FOR_PRESSURE.h"
31  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
32  #include "SOLVE_FOR_PRESSURE3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
33  #include "GW.h"  #include "NH_VARS.h"
34  #endif  #endif
35  #ifdef ALLOW_OBCS  #ifdef ALLOW_CD_CODE
36  #include "OBCS.h"  #include "CD_CODE_VARS.h"
37  #endif  #endif
 #include "SOLVE_FOR_PRESSURE.h"  
38    
39  C     === Functions ====  C     === Functions ====
40        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
# Line 41  C     === Functions ==== Line 42  C     === Functions ====
42    
43  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
44  C     == Routine arguments ==  C     == Routine arguments ==
45  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
46  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
47  C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE  C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
48        _RL myTime        _RL myTime
49        INTEGER myIter        INTEGER myIter
50        INTEGER myThid        INTEGER myThid
# Line 51  C     myThid - Thread number for this in Line 52  C     myThid - Thread number for this in
52  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
53  C     == Local variables ==  C     == Local variables ==
54        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
55        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        INTEGER ks
56        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        INTEGER numIters
57        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
58        _RL tmpFac        _RL tmpFac
59        INTEGER numIters        _RL sumEmP, tileEmP(nSx,nSy)
60          LOGICAL putPmEinXvector
61          INTEGER ioUnit
62          CHARACTER*10 sufx
63        CHARACTER*(MAX_LEN_MBUF) msgBuf        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  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  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
117          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  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
123            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
124  #endif  #endif
125            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)
126            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
127           ENDDO           ENDDO
128          ENDDO          ENDDO
129          IF (useRealFreshWaterFlux) THEN          IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
130           tmpFac = freeSurfFac*convertEmP2rUnit           tmpFac = freeSurfFac*mass2rUnit
131           IF (exactConserv)           IF (exactConserv)
132       &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow       &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
133           DO j=1,sNy           DO j=1,sNy
134            DO i=1,sNx            DO i=1,sNx
135             cg2d_b(i,j,bi,bj) =             cg2d_b(i,j,bi,bj) =
136       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
137            ENDDO            ENDDO
138           ENDDO           ENDDO
139          ENDIF          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         ENDDO
151        ENDDO        ENDDO
152          IF ( putPmEinXvector ) THEN
153            CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
154          ENDIF
155    
156        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
157         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
158          DO K=Nr,1,-1          IF ( putPmEinXvector ) THEN
159           DO j=1,sNy+1            tmpFac = 0.
160            DO i=1,sNx+1            IF (globalArea.GT.0.) tmpFac =
161             uf(i,j) = _dyG(i,j,bi,bj)       &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
162       &      *drF(k)*_hFacW(i,j,k,bi,bj)            DO j=1,sNy
163             vf(i,j) = _dxG(i,j,bi,bj)             DO i=1,sNx
164       &      *drF(k)*_hFacS(i,j,k,bi,bj)              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
165         &                        - tmpFac*Bo_surf(i,j,bi,bj)
166               ENDDO
167            ENDDO            ENDDO
168           ENDDO          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(           CALL CALC_DIV_GHAT(
173       I       bi,bj,1,sNx,1,sNy,K,       I                       bi,bj,k,
174       I       uf,vf,       U                       cg2d_b, cg3d_b,
175       U       cg2d_b,       I                       myThid )
      I       myThid)  
176          ENDDO          ENDDO
177         ENDDO         ENDDO
178        ENDDO        ENDDO
179    
 C--   Add source term arising from w=d/dt (p_s + p_nh)  
180        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
181         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
182  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
183          IF ( nonHydrostatic ) THEN          IF ( oldFreeSurfTerm ) THEN
184    C--   Add source term arising from w=d/dt (p_s + p_nh)
185           DO j=1,sNy           DO j=1,sNy
186            DO i=1,sNx            DO i=1,sNx
187             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
188       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf             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)       &         *( etaN(i,j,bi,bj)
193       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
194             cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)              cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
195       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
196         &         /deltaTMom/deltaTfreesurf
197       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
198       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
199               ENDIF
200            ENDDO            ENDDO
201           ENDDO           ENDDO
202          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
203  #else  #else
204    C--   Add source term arising from w=d/dt (p_s)
205          IF ( exactConserv ) THEN          IF ( exactConserv ) THEN
206  #endif  #endif /* ALLOW_NONHYDROSTATIC */
207           DO j=1,sNy           DO j=1,sNy
208            DO i=1,sNx            DO i=1,sNx
209               ks = ksurfC(i,j,bi,bj)
210             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
211       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
212         &         /deltaTMom/deltaTfreesurf
213       &         * etaH(i,j,bi,bj)       &         * etaH(i,j,bi,bj)
214            ENDDO            ENDDO
215           ENDDO           ENDDO
216          ELSE          ELSE
217           DO j=1,sNy           DO j=1,sNy
218            DO i=1,sNx            DO i=1,sNx
219               ks = ksurfC(i,j,bi,bj)
220             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
221       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
222         &         /deltaTMom/deltaTfreesurf
223       &         * etaN(i,j,bi,bj)       &         * etaN(i,j,bi,bj)
224            ENDDO            ENDDO
225           ENDDO           ENDDO
226          ENDIF          ENDIF
227    
228  #ifdef ALLOW_OBCS  #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          IF (useOBCS) THEN
          DO i=1,sNx  
 C Northern boundary  
           IF (OB_Jn(I,bi,bj).NE.0) THEN  
            cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.  
            cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.  
           ENDIF  
 C Southern boundary  
           IF (OB_Js(I,bi,bj).NE.0) THEN  
            cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.  
            cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.  
           ENDIF  
          ENDDO  
239           DO j=1,sNy           DO j=1,sNy
240  C Eastern boundary            DO i=1,sNx
241            IF (OB_Ie(J,bi,bj).NE.0) THEN             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*maskInC(i,j,bi,bj)
242             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*maskInC(i,j,bi,bj)
243             cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.           ENDDO
           ENDIF  
 C Western boundary  
           IF (OB_Iw(J,bi,bj).NE.0) THEN  
            cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.  
            cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.  
           ENDIF  
244           ENDDO           ENDDO
245          ENDIF          ENDIF
246  #endif  #endif /* ALLOW_OBCS */
247    C-    end bi,bj loops
248         ENDDO         ENDDO
249        ENDDO        ENDDO
250    
251  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
252        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
253         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
254       &                        myThid)       &                        myThid)
255        ENDIF        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.
# Line 187  C     see CG2D.h for the interface to th Line 265  C     see CG2D.h for the interface to th
265        firstResidual=0.        firstResidual=0.
266        lastResidual=0.        lastResidual=0.
267        numIters=cg2dMaxIters        numIters=cg2dMaxIters
268        CALL CG2D(  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 )
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,       U           cg2d_b,
296       U           cg2d_x,       U           cg2d_x,
297       O           firstResidual,       O           firstResidual,
298       O           lastResidual,       O           lastResidual,
299       U           numIters,       U           numIters,
300       I           myThid )       I           myThid )
301        _EXCH_XY_R8(cg2d_x, myThid )        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  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
307        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
308         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
309       &                        myThid)       &                        myThid)
310        ENDIF        ENDIF
311  #endif  #endif
312    
313  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
314        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
315       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
316         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
317         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
318         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
319         WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
321         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
322         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
323         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
324            _END_MASTER( myThid )
325           ENDIF
326        ENDIF        ENDIF
327    
328  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
# Line 228  C--   Transfert the 2D-solution to "etaN Line 337  C--   Transfert the 2D-solution to "etaN
337        ENDDO        ENDDO
338    
339  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
340        IF ( nonHydrostatic ) THEN        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.
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO j=1,sNy+1  
           DO i=1,sNx+1  
            uf(i,j)=-_recip_dxC(i,j,bi,bj)*  
      &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))  
            vf(i,j)=-_recip_dyC(i,j,bi,bj)*  
      &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))  
           ENDDO  
          ENDDO  
   
 #ifdef ALLOW_OBCS  
          IF (useOBCS) THEN  
           DO i=1,sNx+1  
 C Northern boundary  
           IF (OB_Jn(I,bi,bj).NE.0) THEN  
            vf(I,OB_Jn(I,bi,bj))=0.  
           ENDIF  
 C Southern boundary  
           IF (OB_Js(I,bi,bj).NE.0) THEN  
            vf(I,OB_Js(I,bi,bj)+1)=0.  
           ENDIF  
           ENDDO  
           DO j=1,sNy+1  
 C Eastern boundary  
           IF (OB_Ie(J,bi,bj).NE.0) THEN  
            uf(OB_Ie(J,bi,bj),J)=0.  
           ENDIF  
 C Western boundary  
           IF (OB_Iw(J,bi,bj).NE.0) THEN  
            uf(OB_Iw(J,bi,bj)+1,J)=0.  
           ENDIF  
           ENDDO  
          ENDIF  
 #endif  
   
          K=1  
          DO j=1,sNy  
           DO i=1,sNx  
            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)  
      &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)  
      &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)  
      &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)  
      &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )  
      &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom  
      &          -wVel(i,j,k+1,bi,bj)  
      &        )*_rA(i,j,bi,bj)/deltaTmom  
           ENDDO  
          ENDDO  
          DO K=2,Nr-1  
           DO j=1,sNy  
            DO i=1,sNx  
             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)  
      &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)  
      &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)  
      &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)  
      &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )  
      &       +( wVel(i,j,k  ,bi,bj)  
      &         -wVel(i,j,k+1,bi,bj)  
      &        )*_rA(i,j,bi,bj)/deltaTmom  
   
            ENDDO  
           ENDDO  
          ENDDO  
          K=Nr  
          DO j=1,sNy  
           DO i=1,sNx  
             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)  
      &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)  
      &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)  
      &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)  
      &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )  
      &       +( wVel(i,j,k  ,bi,bj)  
      &        )*_rA(i,j,bi,bj)/deltaTmom  
   
           ENDDO  
          ENDDO  
348    
349  #ifdef ALLOW_OBCS  C--   Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
350           IF (useOBCS) THEN  C                             2) Update or Add free-surface contribution
351            DO K=1,Nr  C                             3) increment in horiz velocity due to new cg2d_x
352            DO i=1,sNx  C                             4) add vertical velocity contribution.
353  C Northern boundary         CALL PRE_CG3D(
354            IF (OB_Jn(I,bi,bj).NE.0) THEN       I                oldFreeSurfTerm,
355             cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.       I                cg2d_x,
356            ENDIF       U                cg3d_b,
357  C Southern boundary       I                myTime, myIter, myThid )
358            IF (OB_Js(I,bi,bj).NE.0) THEN  
359             cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.  #ifdef ALLOW_DEBUG
360            ENDIF         IF ( debugLevel .GE. debLevB ) THEN
361            ENDDO          CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
362            DO j=1,sNy       &                         myThid)
363  C Eastern boundary         ENDIF
           IF (OB_Ie(J,bi,bj).NE.0) THEN  
            cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.  
           ENDIF  
 C Western boundary  
           IF (OB_Iw(J,bi,bj).NE.0) THEN  
            cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.  
           ENDIF  
           ENDDO  
           ENDDO  
          ENDIF  
364  #endif  #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          ENDDO ! bi  C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
399         ENDDO ! bj  C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
400           IF ( nonHydrostatic .AND. exactConserv ) THEN
401        firstResidual=0.          IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
402        lastResidual=0.           WRITE(sufx,'(I10.10)') myIter
403        numIters=cg2dMaxIters           CALL WRITE_FLD_XYZ_RL('cg3d_x.',sufx, phi_nh, myIter,myThid )
404        CALL CG3D(          ENDIF
405       U           cg3d_b,          CALL POST_CG3D(
406       U           phi_nh,       I                  zeroPsNH, zeroMeanPnh,
407       O           firstResidual,       I                  myTime, myIter, myThid )
408       O           lastResidual,         ENDIF
      U           numIters,  
      I           myThid )  
       _EXCH_XYZ_R8(phi_nh, myThid )  
409    
       IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,  
      &                                    myTime-deltaTClock) ) THEN  
        _BEGIN_MASTER( myThid )  
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual  
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
        WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters  
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual  
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
        _END_MASTER( )  
410        ENDIF        ENDIF
411    #endif /* ALLOW_NONHYDROSTATIC */
412    
413        ENDIF  #ifdef ALLOW_SHOWFLOPS
414          CALL SHOWFLOPS_INSOLVE( myThid)
415  #endif  #endif
416    
417        RETURN        RETURN
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.76
  ViewVC Help
Powered by ViewVC 1.1.22