/[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.17 by jmc, Tue Mar 6 16:57:10 2001 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  CStartOfInterface  CBOP
8        SUBROUTINE SOLVE_FOR_PRESSURE( myThid )  C     !ROUTINE: SOLVE_FOR_PRESSURE
9  C     /==========================================================\  C     !INTERFACE:
10  C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |        SUBROUTINE SOLVE_FOR_PRESSURE( myTime, myIter, myThid )
11  C     | o Controls inversion of two and/or three-dimensional     |  
12  C     |   elliptic problems for the pressure field.              |  C     !DESCRIPTION: \bv
13  C     \==========================================================/  C     *==========================================================*
14        IMPLICIT NONE  C     | SUBROUTINE SOLVE_FOR_PRESSURE
15    C     | o Controls inversion of two and/or three-dimensional
16    C     |   elliptic problems for the pressure field.
17    C     *==========================================================*
18    C     \ev
19    
20    C     !USES:
21          IMPLICIT NONE
22  C     == Global variables  C     == Global variables
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
26  #include "GRID.h"  #include "GRID.h"
27  #include "SURFACE.h"  #include "SURFACE.h"
28    #include "FFIELDS.h"
29    #include "DYNVARS.h"
30    #include "SOLVE_FOR_PRESSURE.h"
31  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
32  #include "CG3D.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
41    
42    C     === Functions ====
43          LOGICAL  DIFFERENT_MULTIPLE
44          EXTERNAL DIFFERENT_MULTIPLE
45    
46    C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myTime :: Current time in simulation
49    C     myIter :: Current iteration number in simulation
50    C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
51          _RL myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
54    
55  C     Local variables  C     !LOCAL VARIABLES:
56  C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector  C     == Local variables ==
 C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  
57        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
58        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL firstResidual,lastResidual
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL tmpFac
60        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL sumEmP, tileEmP(nSx,nSy)
61        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        LOGICAL putPmEinXvector
62          INTEGER numIters, ks, ioUnit
63          CHARACTER*10 sufx
64          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
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) = 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.
 #ifdef USE_NATURAL_BCS  
      &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*  
      &       EmPmR(I,J,bi,bj)/deltaTMom  
 #endif  
133           ENDDO           ENDDO
134          ENDDO          ENDDO
135            IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
136             tmpFac = freeSurfFac*mass2rUnit
137             IF (exactConserv)
138         &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
139             DO j=1,sNy
140              DO i=1,sNx
141               cg2d_b(i,j,bi,bj) =
142         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
143              ENDDO
144             ENDDO
145            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 84  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          DO j=1,sNy  C--   Add EmPmR contribution to top level cg3d_b:
191           DO i=1,sNx  C      (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
192            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)          IF ( use3Dsolver .AND.
193       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &       useRealFreshWaterFlux.AND.fluidIsWater ) THEN
194       &          -cg2d_x(I,J,bi,bj)            tmpFac = freeSurfFac*mass2rUnit
195       &          -cg3d_x(I,J,1,bi,bj)            IF (exactConserv)
196       &        )/deltaTMom/deltaTMom       &     tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
197            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)            ks = 1
198       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(            IF ( usingPCoords ) ks = Nr
199       &         -cg2d_x(I,J,bi,bj)            DO j=1,sNy
200       &         -cg3d_x(I,J,1,bi,bj)             DO i=1,sNx
201       &       )/deltaTMom/deltaTMom              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
209              DO i=1,sNx
210               ks = ksurfC(i,j,bi,bj)
211               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)
216         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
217                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
218         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
219         &         /deltaTMom/deltaTfreesurf
220         &         *( etaN(i,j,bi,bj)
221         &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
222               ENDIF
223              ENDDO
224           ENDDO           ENDDO
225          ENDDO          ELSEIF ( exactConserv ) THEN
226  #else  #else
227          DO j=1,sNy          IF ( exactConserv ) THEN
228           DO i=1,sNx  #endif /* ALLOW_NONHYDROSTATIC */
229            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)           DO j=1,sNy
230       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(            DO i=1,sNx
231       &          -cg2d_x(I,J,bi,bj)             ks = ksurfC(i,j,bi,bj)
232       &        )/deltaTMom/deltaTMom             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
233         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
234         &         /deltaTMom/deltaTfreesurf
235         &         * etaH(i,j,bi,bj)
236              ENDDO
237           ENDDO           ENDDO
238          ENDDO          ELSE
239  #endif           DO j=1,sNy
240              DO i=1,sNx
241               ks = ksurfC(i,j,bi,bj)
242               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
243         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
244         &         /deltaTMom/deltaTfreesurf
245         &         * etaN(i,j,bi,bj)
246              ENDDO
247             ENDDO
248            ENDIF
249    
250  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
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.
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.
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.
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.
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    #ifdef ALLOW_DEBUG
283          IF ( debugLevel .GE. debLevB ) THEN
284           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
285         &                        myThid)
286          ENDIF
287    #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.
295  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
296        CALL CG2D(        firstResidual=0.
297       I           cg2d_b,        lastResidual=0.
298          numIters=cg2dMaxIters
299    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,
315         U           cg2d_x,
316         O           firstResidual,
317         O           lastResidual,
318         U           numIters,
319         I           myThid )
320          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,       U           cg2d_x,
328         O           firstResidual,
329         O           lastResidual,
330         U           numIters,
331       I           myThid )       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    #ifdef ALLOW_DEBUG
338          IF ( debugLevel .GE. debLevB ) THEN
339           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
340         &                        myThid)
341          ENDIF
342    #endif
343    
344        _EXCH_XY_R8(cg2d_x, myThid )  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
345          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
346         &   ) THEN
347           IF ( debugLevel .GE. debLevA ) THEN
348            _BEGIN_MASTER( myThid )
349            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
350            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
351            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
352            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
353            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
354            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
355            _END_MASTER( myThid )
356           ENDIF
357          ENDIF
358    
359  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
360        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
361         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
362          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
363           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
364            etaN(i,j,bi,bj) = cg2d_x(i,j,bi,bj)            etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
365           ENDDO           ENDDO
366          ENDDO          ENDDO
367         ENDDO         ENDDO
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)=-gBaro*_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)=-gBaro*_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 178  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       &       +(       &       +( 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
      &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(  
      &          +cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/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 /* ALLOW_OBCS */
536    C-    end bi,bj loops
537            ENDDO
538           ENDDO
539    
540    #ifdef ALLOW_DEBUG
541          IF ( debugLevel .GE. debLevB ) THEN
542           CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
543         &                        myThid)
544          ENDIF
545  #endif  #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          ENDDO ! bi        firstResidual=0.
552         ENDDO ! bj        lastResidual=0.
553          numIters=cg3dMaxIters
554          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
555          CALL CG3D(
556         U           cg3d_b,
557         U           phi_nh,
558         O           firstResidual,
559         O           lastResidual,
560         U           numIters,
561         I           myIter, myThid )
562          _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         CALL CG3D( myThid )  C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
580         _EXCH_XYZ_R8(cg3d_x, myThid )        IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
581           IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
582            WRITE(sufx,'(I10.10)') myIter
583            CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx, phi_nh, myIter, myThid )
584           ENDIF
585           DO bj=myByLo(myThid),myByHi(myThid)
586            DO bi=myBxLo(myThid),myBxHi(myThid)
587    
588             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
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.17  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22