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