/[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.72 by jmc, Mon Nov 30 19:20:07 2009 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
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37  #endif  #endif
38  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
39  #include "OBCS.h"  #include "OBCS.h"
40  #endif  #endif
 #include "SOLVE_FOR_PRESSURE.h"  
41    
42  C     === Functions ====  C     === Functions ====
43        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
# Line 41  C     === Functions ==== Line 45  C     === Functions ====
45    
46  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
49  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
50  C     myThid - Thread number for this instance of SOLVE_FOR_PRESSURE  C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
51        _RL myTime        _RL myTime
52        INTEGER myIter        INTEGER myIter
53        INTEGER myThid        INTEGER myThid
# Line 51  C     myThid - Thread number for this in Line 55  C     myThid - Thread number for this in
55  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
56  C     == Local variables ==  C     == Local variables ==
57        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
       _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
58        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
59        _RL tmpFac        _RL tmpFac
60        INTEGER numIters        _RL sumEmP, tileEmP(nSx,nSy)
61          LOGICAL putPmEinXvector
62          INTEGER numIters, ks, ioUnit
63          CHARACTER*10 sufx
64        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
65    #ifdef ALLOW_NONHYDROSTATIC
66          INTEGER kp1
67          _RL     wFacKm, wFacKp
68          LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
69          _RL     tmpVar(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70          _RL     uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71          _RL     vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72    #else
73          _RL     cg3d_b(1)
74    #endif
75  CEOP  CEOP
76    
77    #ifdef ALLOW_NONHYDROSTATIC
78            zeroPsNH = .FALSE.
79    c       zeroPsNH = use3Dsolver .AND. exactConserv
80    c    &                         .AND. select_rStar.EQ.0
81            zeroMeanPnh = .FALSE.
82    c       zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
83    c       oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
84    c    &                                .AND. .NOT.zeroPsNH
85            oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
86    #else
87            cg3d_b(1) = 0.
88    #endif
89    
90    C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
91    C anelastic (always Z-coordinate):
92    C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
93    C        (this reduces the number of lines of code to modify)
94    C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
95    C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
96    C       => 2 factors cancel in elliptic eq. for Phi_s ,
97    C       but 1rst factor(a) remains in RHS cg2d_b.
98    
99    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
100    C     instead of simply etaN ; This can speed-up the solver convergence in
101    C     the case where |Global_mean_PmE| is large.
102          putPmEinXvector = .FALSE.
103    c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
104    
105          IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
106            _BEGIN_MASTER( myThid )
107            ioUnit = standardMessageUnit
108            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
109         &       ' putPmEinXvector =', putPmEinXvector
110            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
111    #ifdef ALLOW_NONHYDROSTATIC
112            WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
113         &       ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
114            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
115            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
116         &       ' oldFreeSurfTerm =', oldFreeSurfTerm
117            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
118    #endif
119            _END_MASTER( myThid )
120          ENDIF
121    
122  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
123          sumEmP = 0.
124        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
125         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
126          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
127           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
128  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
129            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
130  #endif  #endif
131            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)
132            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
133           ENDDO           ENDDO
134          ENDDO          ENDDO
135          IF (useRealFreshWaterFlux) THEN          IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
136           tmpFac = freeSurfFac*convertEmP2rUnit           tmpFac = freeSurfFac*mass2rUnit
137           IF (exactConserv)           IF (exactConserv)
138       &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow       &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
139           DO j=1,sNy           DO j=1,sNy
140            DO i=1,sNx            DO i=1,sNx
141             cg2d_b(i,j,bi,bj) =             cg2d_b(i,j,bi,bj) =
142       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
143            ENDDO            ENDDO
144           ENDDO           ENDDO
145          ENDIF          ENDIF
146            IF ( putPmEinXvector ) THEN
147             tileEmP(bi,bj) = 0.
148             DO j=1,sNy
149              DO i=1,sNx
150                tileEmP(bi,bj) = tileEmP(bi,bj)
151         &                     + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
152         &                                    *maskH(i,j,bi,bj)
153              ENDDO
154             ENDDO
155            ENDIF
156         ENDDO         ENDDO
157        ENDDO        ENDDO
158          IF ( putPmEinXvector ) THEN
159            CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
160          ENDIF
161    
162        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
163         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
164          DO K=Nr,1,-1          IF ( putPmEinXvector ) THEN
165           DO j=1,sNy+1            tmpFac = 0.
166            DO i=1,sNx+1            IF (globalArea.GT.0.) tmpFac =
167             uf(i,j) = _dyG(i,j,bi,bj)       &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
168       &      *drF(k)*_hFacW(i,j,k,bi,bj)            DO j=1,sNy
169             vf(i,j) = _dxG(i,j,bi,bj)             DO i=1,sNx
170       &      *drF(k)*_hFacS(i,j,k,bi,bj)              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
171         &                        - tmpFac*Bo_surf(i,j,bi,bj)
172               ENDDO
173            ENDDO            ENDDO
174           ENDDO          ENDIF
175    C- RHS: similar to the divergence of the vertically integrated mass transport:
176    C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
177            DO k=Nr,1,-1
178           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
179       I       bi,bj,1,sNx,1,sNy,K,       I                       bi,bj,k,
180       I       uf,vf,       U                       cg2d_b, cg3d_b,
181       U       cg2d_b,       I                       myThid )
      I       myThid)  
182          ENDDO          ENDDO
183         ENDDO         ENDDO
184        ENDDO        ENDDO
# Line 109  C--   Add source term arising from w=d/d Line 187  C--   Add source term arising from w=d/d
187        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
188         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
189  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
190          IF ( nonHydrostatic ) THEN  C--   Add EmPmR contribution to top level cg3d_b:
191    C      (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
192            IF ( use3Dsolver .AND.
193         &       useRealFreshWaterFlux.AND.fluidIsWater ) THEN
194              tmpFac = freeSurfFac*mass2rUnit
195              IF (exactConserv)
196         &     tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
197              ks = 1
198              IF ( usingPCoords ) ks = Nr
199              DO j=1,sNy
200               DO i=1,sNx
201                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
202         &        + tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
203               ENDDO
204              ENDDO
205            ENDIF
206    
207            IF ( oldFreeSurfTerm ) THEN
208           DO j=1,sNy           DO j=1,sNy
209            DO i=1,sNx            DO i=1,sNx
210             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
211       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf             IF ( ks.LE.Nr ) THEN
212                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
213         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
214         &         /deltaTMom/deltaTfreesurf
215       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
216       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
217             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)
218       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
219         &         /deltaTMom/deltaTfreesurf
220       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
221       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
222               ENDIF
223            ENDDO            ENDDO
224           ENDDO           ENDDO
225          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
226  #else  #else
227          IF ( exactConserv ) THEN          IF ( exactConserv ) THEN
228  #endif  #endif /* ALLOW_NONHYDROSTATIC */
229           DO j=1,sNy           DO j=1,sNy
230            DO i=1,sNx            DO i=1,sNx
231               ks = ksurfC(i,j,bi,bj)
232             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
233       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
234         &         /deltaTMom/deltaTfreesurf
235       &         * etaH(i,j,bi,bj)       &         * etaH(i,j,bi,bj)
236            ENDDO            ENDDO
237           ENDDO           ENDDO
238          ELSE          ELSE
239           DO j=1,sNy           DO j=1,sNy
240            DO i=1,sNx            DO i=1,sNx
241               ks = ksurfC(i,j,bi,bj)
242             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
243       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
244         &         /deltaTMom/deltaTfreesurf
245       &         * etaN(i,j,bi,bj)       &         * etaN(i,j,bi,bj)
246            ENDDO            ENDDO
247           ENDDO           ENDDO
# Line 147  C--   Add source term arising from w=d/d Line 251  C--   Add source term arising from w=d/d
251          IF (useOBCS) THEN          IF (useOBCS) THEN
252           DO i=1,sNx           DO i=1,sNx
253  C Northern boundary  C Northern boundary
254            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(i,bi,bj).NE.0) THEN
255             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
256             cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
257            ENDIF            ENDIF
258  C Southern boundary  C Southern boundary
259            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(i,bi,bj).NE.0) THEN
260             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
261             cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
262            ENDIF            ENDIF
263           ENDDO           ENDDO
264           DO j=1,sNy           DO j=1,sNy
265  C Eastern boundary  C Eastern boundary
266            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(j,bi,bj).NE.0) THEN
267             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
268             cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
269            ENDIF            ENDIF
270  C Western boundary  C Western boundary
271            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(j,bi,bj).NE.0) THEN
272             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
273             cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
274            ENDIF            ENDIF
275           ENDDO           ENDDO
276          ENDIF          ENDIF
277  #endif  #endif /* ALLOW_OBCS */
278    C-    end bi,bj loops
279         ENDDO         ENDDO
280        ENDDO        ENDDO
281    
282  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
283        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
284         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
285       &                        myThid)       &                        myThid)
286        ENDIF        ENDIF
287  #endif  #endif
288          IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
289           WRITE(sufx,'(I10.10)') myIter
290           CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
291          ENDIF
292    
293  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
294  C--   gradient solver.  C--   gradient solver.
# Line 187  C     see CG2D.h for the interface to th Line 296  C     see CG2D.h for the interface to th
296        firstResidual=0.        firstResidual=0.
297        lastResidual=0.        lastResidual=0.
298        numIters=cg2dMaxIters        numIters=cg2dMaxIters
299        CALL CG2D(  c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
300    #ifdef ALLOW_CG2D_NSA
301    C--   Call the not-self-adjoint version of cg2d
302          CALL CG2D_NSA(
303         U           cg2d_b,
304         U           cg2d_x,
305         O           firstResidual,
306         O           lastResidual,
307         U           numIters,
308         I           myThid )
309    #else /* not ALLOW_CG2D_NSA = default */
310    #ifdef ALLOW_SRCG
311          IF ( useSRCGSolver ) THEN
312    C--   Call the single reduce CG solver
313           CALL CG2D_SR(
314       U           cg2d_b,       U           cg2d_b,
315       U           cg2d_x,       U           cg2d_x,
316       O           firstResidual,       O           firstResidual,
317       O           lastResidual,       O           lastResidual,
318       U           numIters,       U           numIters,
319       I           myThid )       I           myThid )
320        _EXCH_XY_R8(cg2d_x, myThid )        ELSE
321    #else
322          IF (.TRUE.) THEN
323    C--   Call the default CG solver
324    #endif /* ALLOW_SRCG */
325           CALL CG2D(
326         U           cg2d_b,
327         U           cg2d_x,
328         O           firstResidual,
329         O           lastResidual,
330         U           numIters,
331         I           myThid )
332          ENDIF
333    #endif /* ALLOW_CG2D_NSA */
334          _EXCH_XY_RL( cg2d_x, myThid )
335    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
336    
337  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
338        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
339         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
340       &                        myThid)       &                        myThid)
341        ENDIF        ENDIF
342  #endif  #endif
343    
344  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) :
345        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
346       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
347         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
348         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
349         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
350         WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
351         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
352         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
353         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
354         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
355            _END_MASTER( myThid )
356           ENDIF
357        ENDIF        ENDIF
358    
359  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
# Line 228  C--   Transfert the 2D-solution to "etaN Line 368  C--   Transfert the 2D-solution to "etaN
368        ENDDO        ENDDO
369    
370  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
371        IF ( nonHydrostatic ) THEN        IF ( use3Dsolver ) THEN
372           IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
373            WRITE(sufx,'(I10.10)') myIter
374            CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
375           ENDIF
376    
377  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 ).
378  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
379         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
380          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
381    
382    C--   Update or Add free-surface contribution to cg3d_b:
383    c        IF ( select_rStar.EQ.0 .AND. exactConserv ) THEN
384             IF ( select_rStar.EQ.0 .AND. .NOT.oldFreeSurfTerm ) THEN
385               tmpFac = 0.
386               DO j=1,sNy
387                DO i=1,sNx
388                  ks = ksurfC(i,j,bi,bj)
389                  IF ( ks.LE.Nr ) THEN
390                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
391         &                  +freeSurfFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
392         &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
393         &                              /deltaTMom/deltaTfreesurf
394                  ENDIF
395                ENDDO
396               ENDDO
397    #ifdef NONLIN_FRSURF
398             ELSEIF ( select_rStar.NE.0 ) THEN
399               tmpFac = 0.
400               DO j=1,sNy
401                DO i=1,sNx
402                  ks = ksurfC(i,j,bi,bj)
403                  tmpVar(i,j) = freeSurfFac
404         &                    *( etaN(i,j,bi,bj) - etaH(i,j,bi,bj) )
405         &                    *_rA(i,j,bi,bj)*deepFac2F(ks)
406         &                    /deltaTMom/deltaTfreesurf
407         &                    *recip_Rcol(i,j,bi,bj)
408                ENDDO
409               ENDDO
410               DO k=1,Nr
411                DO j=1,sNy
412                 DO i=1,sNx
413                  cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
414         &          + tmpVar(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
415                 ENDDO
416                ENDDO
417               ENDDO
418    #endif /* NONLIN_FRSURF */
419             ELSEIF ( usingZCoords ) THEN
420    C-       Z coordinate: assume surface @ level k=1
421               tmpFac = freeSurfFac*deepFac2F(1)
422             ELSE
423    C-       Other than Z coordinate: no assumption on surface level index
424               tmpFac = 0.
425               DO j=1,sNy
426                DO i=1,sNx
427                  ks = ksurfC(i,j,bi,bj)
428                  IF ( ks.LE.Nr ) THEN
429                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
430         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
431         &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
432                  ENDIF
433                ENDDO
434               ENDDO
435             ENDIF
436    
437    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
438    
439    C--   Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
440    C                             2) add vertical velocity contribution.
441           DO j=1,sNy+1           DO j=1,sNy+1
442            DO i=1,sNx+1            DO i=1,sNx+1
443             uf(i,j)=-_recip_dxC(i,j,bi,bj)*             uf(i,j) = -_recip_dxC(i,j,bi,bj)
444       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &             * implicSurfPress*implicDiv2DFlow
445             vf(i,j)=-_recip_dyC(i,j,bi,bj)*       &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
446       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))             vf(i,j) = -_recip_dyC(i,j,bi,bj)
447         &             * implicSurfPress*implicDiv2DFlow
448         &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
449            ENDDO            ENDDO
450           ENDDO           ENDDO
451    
# Line 247  C     see CG3D.h for the interface to th Line 453  C     see CG3D.h for the interface to th
453           IF (useOBCS) THEN           IF (useOBCS) THEN
454            DO i=1,sNx+1            DO i=1,sNx+1
455  C Northern boundary  C Northern boundary
456            IF (OB_Jn(I,bi,bj).NE.0) THEN             IF (OB_Jn(i,bi,bj).NE.0)
457             vf(I,OB_Jn(I,bi,bj))=0.       &      vf(i,OB_Jn(i,bi,bj)) = 0.
           ENDIF  
458  C Southern boundary  C Southern boundary
459            IF (OB_Js(I,bi,bj).NE.0) THEN             IF (OB_Js(i,bi,bj).NE.0)
460             vf(I,OB_Js(I,bi,bj)+1)=0.       &      vf(i,OB_Js(i,bi,bj)+1) = 0.
           ENDIF  
461            ENDDO            ENDDO
462            DO j=1,sNy+1            DO j=1,sNy+1
463  C Eastern boundary  C Eastern boundary
464            IF (OB_Ie(J,bi,bj).NE.0) THEN             IF (OB_Ie(j,bi,bj).NE.0)
465             uf(OB_Ie(J,bi,bj),J)=0.       &      uf(OB_Ie(j,bi,bj),j) = 0.
           ENDIF  
466  C Western boundary  C Western boundary
467            IF (OB_Iw(J,bi,bj).NE.0) THEN             IF (OB_Iw(j,bi,bj).NE.0)
468             uf(OB_Iw(J,bi,bj)+1,J)=0.       &      uf(OB_Iw(j,bi,bj)+1,j) = 0.
           ENDIF  
469            ENDDO            ENDDO
470           ENDIF           ENDIF
471  #endif  #endif /* ALLOW_OBCS */
472    
473           K=1  C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
474    C       uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
475    C       always recomputed from continuity eq (like eta when exactConserv=T)
476             k=1
477             kp1 = MIN(k+1,Nr)
478             wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
479             IF (k.GE.Nr) wFacKp = 0.
480           DO j=1,sNy           DO j=1,sNy
481            DO i=1,sNx            DO i=1,sNx
482             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
483       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
484       &       -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)
485       &       +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)
486       &       -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 )
487       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
488       &          -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
489       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
490            ENDDO            ENDDO
491           ENDDO           ENDDO
492           DO K=2,Nr-1           DO k=2,Nr
493              kp1 = MIN(k+1,Nr)
494    C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
495    C        both appear in wVel term, but at 2 different levels
496              wFacKm = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k )
497              wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
498              IF (k.GE.Nr) wFacKp = 0.
499            DO j=1,sNy            DO j=1,sNy
500             DO i=1,sNx             DO i=1,sNx
501              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
502       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
503       &       -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)
504       &       +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)
505       &       -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 )
506       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
507       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
508       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
509    
510             ENDDO             ENDDO
511            ENDDO            ENDDO
512           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  
513    
514  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
515           IF (useOBCS) THEN           IF (useOBCS) THEN
516            DO K=1,Nr            DO k=1,Nr
517            DO i=1,sNx             DO i=1,sNx
518  C Northern boundary  C Northern boundary
519            IF (OB_Jn(I,bi,bj).NE.0) THEN              IF (OB_Jn(i,bi,bj).NE.0)
520             cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.       &       cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj) = 0.
           ENDIF  
521  C Southern boundary  C Southern boundary
522            IF (OB_Js(I,bi,bj).NE.0) THEN              IF (OB_Js(i,bi,bj).NE.0)
523             cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.       &       cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj) = 0.
524            ENDIF             ENDDO
525            ENDDO             DO j=1,sNy
           DO j=1,sNy  
526  C Eastern boundary  C Eastern boundary
527            IF (OB_Ie(J,bi,bj).NE.0) THEN              IF (OB_Ie(j,bi,bj).NE.0)
528             cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.       &       cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj) = 0.
           ENDIF  
529  C Western boundary  C Western boundary
530            IF (OB_Iw(J,bi,bj).NE.0) THEN              IF (OB_Iw(j,bi,bj).NE.0)
531             cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.       &       cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj) = 0.
532            ENDIF             ENDDO
           ENDDO  
533            ENDDO            ENDDO
534           ENDIF           ENDIF
535  #endif  #endif /* ALLOW_OBCS */
536    C-    end bi,bj loops
537            ENDDO
538           ENDDO
539    
540          ENDDO ! bi  #ifdef ALLOW_DEBUG
541         ENDDO ! bj        IF ( debugLevel .GE. debLevB ) THEN
542           CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
543         &                        myThid)
544          ENDIF
545    #endif
546          IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
547           WRITE(sufx,'(I10.10)') myIter
548           CALL WRITE_FLD_XYZ_RL( 'cg3d_b.',sufx, cg3d_b, myIter, myThid )
549          ENDIF
550    
551        firstResidual=0.        firstResidual=0.
552        lastResidual=0.        lastResidual=0.
553        numIters=cg2dMaxIters        numIters=cg3dMaxIters
554          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
555        CALL CG3D(        CALL CG3D(
556       U           cg3d_b,       U           cg3d_b,
557       U           phi_nh,       U           phi_nh,
558       O           firstResidual,       O           firstResidual,
559       O           lastResidual,       O           lastResidual,
560       U           numIters,       U           numIters,
561       I           myThid )       I           myIter, myThid )
562        _EXCH_XYZ_R8(phi_nh, myThid )        _EXCH_XYZ_RL( phi_nh, myThid )
563          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
564    
565          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
566         &   ) THEN
567           IF ( debugLevel .GE. debLevA ) THEN
568            _BEGIN_MASTER( myThid )
569            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
570            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
571            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
572            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
573            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
574            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
575            _END_MASTER( myThid )
576           ENDIF
577          ENDIF
578    
579        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,  C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
580       &                                    myTime-deltaTClock) ) THEN        IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
581         _BEGIN_MASTER( myThid )         IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
582         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual          WRITE(sufx,'(I10.10)') myIter
583         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx, phi_nh, myIter, myThid )
584         WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters         ENDIF
585         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         DO bj=myByLo(myThid),myByHi(myThid)
586         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual          DO bi=myBxLo(myThid),myBxHi(myThid)
587         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
588         _END_MASTER( )           IF ( zeroPsNH .AND. usingZCoords ) THEN
589    C-       Z coordinate: assume surface @ level k=1
590              DO j=1-OLy,sNy+OLy
591               DO i=1-OLx,sNx+OLx
592                 tmpVar(i,j) = phi_nh(i,j,1,bi,bj)
593               ENDDO
594              ENDDO
595             ELSEIF ( zeroPsNH ) THEN
596    C-       Other than Z coordinate: no assumption on surface level index
597              DO j=1-OLy,sNy+OLy
598               DO i=1-OLx,sNx+OLx
599                ks = ksurfC(i,j,bi,bj)
600                IF ( ks.LE.Nr ) THEN
601                 tmpVar(i,j) = phi_nh(i,j,ks,bi,bj)
602                ELSE
603                 tmpVar(i,j) = 0.
604                ENDIF
605               ENDDO
606              ENDDO
607    #ifdef NONLIN_FRSURF
608             ELSE
609    C        zeroMeanPnh : transfert vertical average of P_NH to EtaN
610              DO j=1-OLy,sNy+OLy
611               DO i=1-OLx,sNx+OLx
612                 tmpVar(i,j) = 0.
613               ENDDO
614              ENDDO
615              DO k=1,Nr
616               DO j=1-OLy,sNy+OLy
617                DO i=1-OLx,sNx+OLx
618                 tmpVar(i,j) = tmpVar(i,j)
619         &         + phi_nh(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
620                ENDDO
621               ENDDO
622              ENDDO
623              DO j=1-OLy,sNy+OLy
624               DO i=1-OLx,sNx+OLx
625                 tmpVar(i,j) = tmpVar(i,j)*recip_Rcol(i,j,bi,bj)
626               ENDDO
627              ENDDO
628    #endif /* NONLIN_FRSURF */
629             ENDIF
630             DO k=1,Nr
631              DO j=1-OLy,sNy+OLy
632               DO i=1-OLx,sNx+OLx
633                phi_nh(i,j,k,bi,bj) = ( phi_nh(i,j,k,bi,bj)
634         &                            - tmpVar(i,j)
635         &                            )*maskC(i,j,k,bi,bj)
636               ENDDO
637              ENDDO
638             ENDDO
639             DO j=1-OLy,sNy+OLy
640              DO i=1-OLx,sNx+OLx
641                etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
642         &                      *( cg2d_x(i,j,bi,bj) + tmpVar(i,j) )
643              ENDDO
644             ENDDO
645    
646            ENDDO
647           ENDDO
648        ENDIF        ENDIF
649    
650        ENDIF        ENDIF
651    #endif /* ALLOW_NONHYDROSTATIC */
652    
653    #ifdef ALLOW_SHOWFLOPS
654          CALL SHOWFLOPS_INSOLVE( myThid)
655  #endif  #endif
656    
657        RETURN        RETURN
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22