/[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.74 by jmc, Mon Dec 21 00:24:58 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
58        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        INTEGER ks
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        INTEGER numIters
60        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
61        _RL tmpFac        _RL tmpFac
62        INTEGER numIters        _RL sumEmP, tileEmP(nSx,nSy)
63          LOGICAL putPmEinXvector
64          INTEGER ioUnit
65          CHARACTER*10 sufx
66        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
67    #ifdef ALLOW_NONHYDROSTATIC
68          LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
69    #else
70          _RL     cg3d_b(1)
71    #endif
72  CEOP  CEOP
73    
74    #ifdef ALLOW_NONHYDROSTATIC
75            zeroPsNH = .FALSE.
76    c       zeroPsNH = use3Dsolver .AND. exactConserv
77    c    &                         .AND. select_rStar.EQ.0
78            zeroMeanPnh = .FALSE.
79    c       zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
80    c       oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
81    c    &                                .AND. .NOT.zeroPsNH
82            oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
83    #else
84            cg3d_b(1) = 0.
85    #endif
86    
87    C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
88    C anelastic (always Z-coordinate):
89    C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
90    C        (this reduces the number of lines of code to modify)
91    C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
92    C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
93    C       => 2 factors cancel in elliptic eq. for Phi_s ,
94    C       but 1rst factor(a) remains in RHS cg2d_b.
95    
96    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
97    C     instead of simply etaN ; This can speed-up the solver convergence in
98    C     the case where |Global_mean_PmE| is large.
99          putPmEinXvector = .FALSE.
100    c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
101    
102          IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
103            _BEGIN_MASTER( myThid )
104            ioUnit = standardMessageUnit
105            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
106         &       ' putPmEinXvector =', putPmEinXvector
107            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
108    #ifdef ALLOW_NONHYDROSTATIC
109            WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
110         &       ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
111            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
112            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
113         &       ' oldFreeSurfTerm =', oldFreeSurfTerm
114            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
115    #endif
116            _END_MASTER( myThid )
117          ENDIF
118    
119  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
120          sumEmP = 0.
121        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
122         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
123          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
124           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
125  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
126            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
127  #endif  #endif
128            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)
129            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
130           ENDDO           ENDDO
131          ENDDO          ENDDO
132          IF (useRealFreshWaterFlux) THEN          IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
133           tmpFac = freeSurfFac*convertEmP2rUnit           tmpFac = freeSurfFac*mass2rUnit
134           IF (exactConserv)           IF (exactConserv)
135       &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow       &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
136           DO j=1,sNy           DO j=1,sNy
137            DO i=1,sNx            DO i=1,sNx
138             cg2d_b(i,j,bi,bj) =             cg2d_b(i,j,bi,bj) =
139       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
140            ENDDO            ENDDO
141           ENDDO           ENDDO
142          ENDIF          ENDIF
143            IF ( putPmEinXvector ) THEN
144             tileEmP(bi,bj) = 0.
145             DO j=1,sNy
146              DO i=1,sNx
147                tileEmP(bi,bj) = tileEmP(bi,bj)
148         &                     + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
149         &                                    *maskInC(i,j,bi,bj)
150              ENDDO
151             ENDDO
152            ENDIF
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155          IF ( putPmEinXvector ) THEN
156            CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
157          ENDIF
158    
159        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
160         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
161          DO K=Nr,1,-1          IF ( putPmEinXvector ) THEN
162           DO j=1,sNy+1            tmpFac = 0.
163            DO i=1,sNx+1            IF (globalArea.GT.0.) tmpFac =
164             uf(i,j) = _dyG(i,j,bi,bj)       &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
165       &      *drF(k)*_hFacW(i,j,k,bi,bj)            DO j=1,sNy
166             vf(i,j) = _dxG(i,j,bi,bj)             DO i=1,sNx
167       &      *drF(k)*_hFacS(i,j,k,bi,bj)              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
168         &                        - tmpFac*Bo_surf(i,j,bi,bj)
169               ENDDO
170            ENDDO            ENDDO
171           ENDDO          ENDIF
172    C- RHS: similar to the divergence of the vertically integrated mass transport:
173    C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
174            DO k=Nr,1,-1
175           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
176       I       bi,bj,1,sNx,1,sNy,K,       I                       bi,bj,k,
177       I       uf,vf,       U                       cg2d_b, cg3d_b,
178       U       cg2d_b,       I                       myThid )
      I       myThid)  
179          ENDDO          ENDDO
180         ENDDO         ENDDO
181        ENDDO        ENDDO
182    
 C--   Add source term arising from w=d/dt (p_s + p_nh)  
183        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
184         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
185  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
186          IF ( nonHydrostatic ) THEN          IF ( oldFreeSurfTerm ) THEN
187    C--   Add source term arising from w=d/dt (p_s + p_nh)
188           DO j=1,sNy           DO j=1,sNy
189            DO i=1,sNx            DO i=1,sNx
190             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
191       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf             IF ( ks.LE.Nr ) THEN
192                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
193         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
194         &         /deltaTMom/deltaTfreesurf
195       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
196       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
197             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)
198       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
199         &         /deltaTMom/deltaTfreesurf
200       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
201       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
202               ENDIF
203            ENDDO            ENDDO
204           ENDDO           ENDDO
205          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
206  #else  #else
207    C--   Add source term arising from w=d/dt (p_s)
208          IF ( exactConserv ) THEN          IF ( exactConserv ) THEN
209  #endif  #endif /* ALLOW_NONHYDROSTATIC */
210           DO j=1,sNy           DO j=1,sNy
211            DO i=1,sNx            DO i=1,sNx
212               ks = ksurfC(i,j,bi,bj)
213             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
214       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
215         &         /deltaTMom/deltaTfreesurf
216       &         * etaH(i,j,bi,bj)       &         * etaH(i,j,bi,bj)
217            ENDDO            ENDDO
218           ENDDO           ENDDO
219          ELSE          ELSE
220           DO j=1,sNy           DO j=1,sNy
221            DO i=1,sNx            DO i=1,sNx
222               ks = ksurfC(i,j,bi,bj)
223             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
224       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
225         &         /deltaTMom/deltaTfreesurf
226       &         * etaN(i,j,bi,bj)       &         * etaN(i,j,bi,bj)
227            ENDDO            ENDDO
228           ENDDO           ENDDO
# Line 147  C--   Add source term arising from w=d/d Line 232  C--   Add source term arising from w=d/d
232          IF (useOBCS) THEN          IF (useOBCS) THEN
233           DO i=1,sNx           DO i=1,sNx
234  C Northern boundary  C Northern boundary
235            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(i,bi,bj).NE.0) THEN
236             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
237             cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
238            ENDIF            ENDIF
239  C Southern boundary  C Southern boundary
240            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(i,bi,bj).NE.0) THEN
241             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
242             cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
243            ENDIF            ENDIF
244           ENDDO           ENDDO
245           DO j=1,sNy           DO j=1,sNy
246  C Eastern boundary  C Eastern boundary
247            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(j,bi,bj).NE.0) THEN
248             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
249             cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
250            ENDIF            ENDIF
251  C Western boundary  C Western boundary
252            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(j,bi,bj).NE.0) THEN
253             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
254             cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
255            ENDIF            ENDIF
256           ENDDO           ENDDO
257          ENDIF          ENDIF
258  #endif  #endif /* ALLOW_OBCS */
259    C-    end bi,bj loops
260         ENDDO         ENDDO
261        ENDDO        ENDDO
262    
263  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
264        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
265         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
266       &                        myThid)       &                        myThid)
267        ENDIF        ENDIF
268  #endif  #endif
269          IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
270           WRITE(sufx,'(I10.10)') myIter
271           CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
272          ENDIF
273    
274  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
275  C--   gradient solver.  C--   gradient solver.
# Line 187  C     see CG2D.h for the interface to th Line 277  C     see CG2D.h for the interface to th
277        firstResidual=0.        firstResidual=0.
278        lastResidual=0.        lastResidual=0.
279        numIters=cg2dMaxIters        numIters=cg2dMaxIters
280        CALL CG2D(  c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
281    #ifdef ALLOW_CG2D_NSA
282    C--   Call the not-self-adjoint version of cg2d
283          CALL CG2D_NSA(
284         U           cg2d_b,
285         U           cg2d_x,
286         O           firstResidual,
287         O           lastResidual,
288         U           numIters,
289         I           myThid )
290    #else /* not ALLOW_CG2D_NSA = default */
291    #ifdef ALLOW_SRCG
292          IF ( useSRCGSolver ) THEN
293    C--   Call the single reduce CG solver
294           CALL CG2D_SR(
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 )        ELSE
302    #else
303          IF (.TRUE.) THEN
304    C--   Call the default CG solver
305    #endif /* ALLOW_SRCG */
306           CALL CG2D(
307         U           cg2d_b,
308         U           cg2d_x,
309         O           firstResidual,
310         O           lastResidual,
311         U           numIters,
312         I           myThid )
313          ENDIF
314    #endif /* ALLOW_CG2D_NSA */
315          _EXCH_XY_RL( cg2d_x, myThid )
316    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
317    
318  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
319        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
320         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
321       &                        myThid)       &                        myThid)
322        ENDIF        ENDIF
323  #endif  #endif
324    
325  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) :
326        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
327       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
328         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
329         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
330         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
331         WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
332         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
333         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
334         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
335         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
336            _END_MASTER( myThid )
337           ENDIF
338        ENDIF        ENDIF
339    
340  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
# Line 228  C--   Transfert the 2D-solution to "etaN Line 349  C--   Transfert the 2D-solution to "etaN
349        ENDDO        ENDDO
350    
351  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
352        IF ( nonHydrostatic ) THEN        IF ( use3Dsolver ) THEN
353           IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
354            WRITE(sufx,'(I10.10)') myIter
355            CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
356           ENDIF
357    
358  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 ).
359  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  
360    
361  #ifdef ALLOW_OBCS  C--   Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
362           IF (useOBCS) THEN  C                             2) Update or Add free-surface contribution
363            DO i=1,sNx+1  C                             3) increment in horiz velocity due to new cg2d_x
364  C Northern boundary  C                             4) add vertical velocity contribution.
365            IF (OB_Jn(I,bi,bj).NE.0) THEN         CALL PRE_CG3D(
366             vf(I,OB_Jn(I,bi,bj))=0.       I                oldFreeSurfTerm,
367            ENDIF       I                cg2d_x,
368  C Southern boundary       U                cg3d_b,
369            IF (OB_Js(I,bi,bj).NE.0) THEN       I                myTime, myIter, myThid )
370             vf(I,OB_Js(I,bi,bj)+1)=0.  
371            ENDIF  #ifdef ALLOW_DEBUG
372            ENDDO         IF ( debugLevel .GE. debLevB ) THEN
373            DO j=1,sNy+1          CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
374  C Eastern boundary       &                         myThid)
375            IF (OB_Ie(J,bi,bj).NE.0) THEN         ENDIF
            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  
   
 #ifdef ALLOW_OBCS  
          IF (useOBCS) THEN  
           DO K=1,Nr  
           DO i=1,sNx  
 C Northern boundary  
           IF (OB_Jn(I,bi,bj).NE.0) THEN  
            cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.  
           ENDIF  
 C Southern boundary  
           IF (OB_Js(I,bi,bj).NE.0) THEN  
            cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.  
           ENDIF  
           ENDDO  
           DO j=1,sNy  
 C Eastern boundary  
           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  
376  #endif  #endif
377           IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
378            WRITE(sufx,'(I10.10)') myIter
379            CALL WRITE_FLD_XYZ_RL('cg3d_b.',sufx, cg3d_b, myIter,myThid )
380           ENDIF
381    
382           firstResidual=0.
383           lastResidual=0.
384           numIters=cg3dMaxIters
385           CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
386           CALL CG3D(
387         U            cg3d_b,
388         U            phi_nh,
389         O            firstResidual,
390         O            lastResidual,
391         U            numIters,
392         I            myIter, myThid )
393           _EXCH_XYZ_RL( phi_nh, myThid )
394           CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
395    
396           IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
397         &    ) THEN
398            IF ( debugLevel .GE. debLevA ) THEN
399             _BEGIN_MASTER( myThid )
400             WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
401             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
402             WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
403             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
404             WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
405             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
406             _END_MASTER( myThid )
407            ENDIF
408           ENDIF
409    
410          ENDDO ! bi  C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
411         ENDDO ! bj  C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
412           IF ( nonHydrostatic .AND. exactConserv ) THEN
413        firstResidual=0.          IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
414        lastResidual=0.           WRITE(sufx,'(I10.10)') myIter
415        numIters=cg2dMaxIters           CALL WRITE_FLD_XYZ_RL('cg3d_x.',sufx, phi_nh, myIter,myThid )
416        CALL CG3D(          ENDIF
417       U           cg3d_b,          CALL POST_CG3D(
418       U           phi_nh,       I                  zeroPsNH, zeroMeanPnh,
419       O           firstResidual,       I                  myTime, myIter, myThid )
420       O           lastResidual,         ENDIF
      U           numIters,  
      I           myThid )  
       _EXCH_XYZ_R8(phi_nh, myThid )  
421    
       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( )  
422        ENDIF        ENDIF
423    #endif /* ALLOW_NONHYDROSTATIC */
424    
425        ENDIF  #ifdef ALLOW_SHOWFLOPS
426          CALL SHOWFLOPS_INSOLVE( myThid)
427  #endif  #endif
428    
429        RETURN        RETURN
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22