/[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.48 by jmc, Tue Nov 8 02:14:10 2005 UTC revision 1.64 by jmc, Sun Aug 24 21:40:18 2008 UTC
# Line 11  C     !INTERFACE: Line 11  C     !INTERFACE:
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 45  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 55  C     myThid - Thread number for this in Line 55  C     myThid - Thread number for this in
55  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
56  C     == Local variables ==  C     == Local variables ==
57        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
       _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
58        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
59        _RL tmpFac        _RL tmpFac
60        _RL sumEmP, tileEmP        _RL sumEmP, tileEmP
61        LOGICAL putPmEinXvector        LOGICAL putPmEinXvector
62        INTEGER numIters        INTEGER numIters, ks
63          CHARACTER*10 sufx
64        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
65    #ifdef ALLOW_NONHYDROSTATIC
66          INTEGER kp1
67          _RL     wFacKm, wFacKp
68          LOGICAL zeroPsNH
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  CEOP
75    
76  #ifdef TIME_PER_TIMESTEP  #ifdef ALLOW_NONHYDROSTATIC
77  CCE107 common block for per timestep timing  c       zeroPsNH = .FALSE.
78  C     !TIMING VARIABLES          zeroPsNH = exactConserv
79  C     == Timing variables ==  #else
80        REAL*8 utnew, utold, stnew, stold, wtnew, wtold          cg3d_b(1) = 0.
       COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold  
81  #endif  #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  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  C     instead of simply etaN ; This can speed-up the solver convergence in
94  C     the case where |Global_mean_PmE| is large.  C     the case where |Global_mean_PmE| is large.
95        putPmEinXvector = .FALSE.        putPmEinXvector = .FALSE.
96  c     putPmEinXvector = useRealFreshWaterFlux  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.        sumEmP = 0.
# Line 92  C--   Save previous solution & Initialis Line 108  C--   Save previous solution & Initialis
108            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
109           ENDDO           ENDDO
110          ENDDO          ENDDO
111          IF (useRealFreshWaterFlux) THEN          IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
112           tmpFac = freeSurfFac*convertEmP2rUnit           tmpFac = freeSurfFac*mass2rUnit
113           IF (exactConserv)           IF (exactConserv)
114       &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow       &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
115           DO j=1,sNy           DO j=1,sNy
116            DO i=1,sNx            DO i=1,sNx
117             cg2d_b(i,j,bi,bj) =             cg2d_b(i,j,bi,bj) =
118       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom       &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
119            ENDDO            ENDDO
120           ENDDO           ENDDO
# Line 123  C--   Save previous solution & Initialis Line 139  C--   Save previous solution & Initialis
139         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
140          IF ( putPmEinXvector ) THEN          IF ( putPmEinXvector ) THEN
141            tmpFac = 0.            tmpFac = 0.
142            IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf            IF (globalArea.GT.0.) tmpFac =
143       &                          *convertEmP2rUnit*sumEmP/globalArea       &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
144            DO j=1,sNy            DO j=1,sNy
145             DO i=1,sNx             DO i=1,sNx
146              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
# Line 132  C--   Save previous solution & Initialis Line 148  C--   Save previous solution & Initialis
148             ENDDO             ENDDO
149            ENDDO            ENDDO
150          ENDIF          ENDIF
151          DO K=Nr,1,-1  C- RHS: similar to the divergence of the vertically integrated mass transport:
152           DO j=1,sNy+1  C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
153            DO i=1,sNx+1          DO k=Nr,1,-1
            uf(i,j) = _dyG(i,j,bi,bj)  
      &      *drF(k)*_hFacW(i,j,k,bi,bj)  
            vf(i,j) = _dxG(i,j,bi,bj)  
      &      *drF(k)*_hFacS(i,j,k,bi,bj)  
           ENDDO  
          ENDDO  
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 154  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          IF ( nonHydrostatic ) THEN          IF ( use3Dsolver .AND. zeroPsNH ) THEN
167           DO j=1,sNy           DO j=1,sNy
168            DO i=1,sNx            DO i=1,sNx
169             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
170       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf             IF ( ks.LE.Nr ) THEN
171                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
172         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
173         &         /deltaTMom/deltaTfreesurf
174         &         * etaH(i,j,bi,bj)
175                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
176         &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
177         &         /deltaTMom/deltaTfreesurf
178         &         * etaH(i,j,bi,bj)
179               ENDIF
180              ENDDO
181             ENDDO
182            ELSEIF ( use3Dsolver ) 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       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
191       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
192             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)
193       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -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               ENDIF
198            ENDDO            ENDDO
199           ENDDO           ENDDO
200          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
# Line 173  C--   Add source term arising from w=d/d Line 203  C--   Add source term arising from w=d/d
203  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
204           DO j=1,sNy           DO j=1,sNy
205            DO i=1,sNx            DO i=1,sNx
206               ks = ksurfC(i,j,bi,bj)
207             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
208       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
209         &         /deltaTMom/deltaTfreesurf
210       &         * etaH(i,j,bi,bj)       &         * etaH(i,j,bi,bj)
211            ENDDO            ENDDO
212           ENDDO           ENDDO
213          ELSE          ELSE
214           DO j=1,sNy           DO j=1,sNy
215            DO i=1,sNx            DO i=1,sNx
216               ks = ksurfC(i,j,bi,bj)
217             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
218       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
219         &         /deltaTMom/deltaTfreesurf
220       &         * etaN(i,j,bi,bj)       &         * etaN(i,j,bi,bj)
221            ENDDO            ENDDO
222           ENDDO           ENDDO
# Line 192  C--   Add source term arising from w=d/d Line 226  C--   Add source term arising from w=d/d
226          IF (useOBCS) THEN          IF (useOBCS) THEN
227           DO i=1,sNx           DO i=1,sNx
228  C Northern boundary  C Northern boundary
229            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(i,bi,bj).NE.0) THEN
230             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
231             cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
232            ENDIF            ENDIF
233  C Southern boundary  C Southern boundary
234            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(i,bi,bj).NE.0) THEN
235             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
236             cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
237            ENDIF            ENDIF
238           ENDDO           ENDDO
239           DO j=1,sNy           DO j=1,sNy
240  C Eastern boundary  C Eastern boundary
241            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(j,bi,bj).NE.0) THEN
242             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
243             cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
244            ENDIF            ENDIF
245  C Western boundary  C Western boundary
246            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(j,bi,bj).NE.0) THEN
247             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
248             cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
249            ENDIF            ENDIF
250           ENDDO           ENDDO
251          ENDIF          ENDIF
252  #endif  #endif /* ALLOW_OBCS */
253    C-    end bi,bj loops
254         ENDDO         ENDDO
255        ENDDO        ENDDO
256    
# Line 225  C Western boundary Line 260  C Western boundary
260       &                        myThid)       &                        myThid)
261        ENDIF        ENDIF
262  #endif  #endif
263          IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
264           WRITE(sufx,'(I10.10)') myIter
265           CALL WRITE_FLD_XY_RS( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
266          ENDIF
267    
268  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
269  C--   gradient solver.  C--   gradient solver.
# Line 232  C     see CG2D.h for the interface to th Line 271  C     see CG2D.h for the interface to th
271        firstResidual=0.        firstResidual=0.
272        lastResidual=0.        lastResidual=0.
273        numIters=cg2dMaxIters        numIters=cg2dMaxIters
274    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
275    #ifdef ALLOW_CG2D_NSA
276    C--   Call the not-self-adjoint version of cg2d
277          CALL CG2D_NSA(
278         U           cg2d_b,
279         U           cg2d_x,
280         O           firstResidual,
281         O           lastResidual,
282         U           numIters,
283         I           myThid )
284    #else /* not ALLOW_CG2D_NSA = default */
285        CALL CG2D(        CALL CG2D(
286       U           cg2d_b,       U           cg2d_b,
287       U           cg2d_x,       U           cg2d_x,
# Line 239  C     see CG2D.h for the interface to th Line 289  C     see CG2D.h for the interface to th
289       O           lastResidual,       O           lastResidual,
290       U           numIters,       U           numIters,
291       I           myThid )       I           myThid )
292    #endif /* ALLOW_CG2D_NSA */
293        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
294    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
295    
296  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
297        IF ( debugLevel .GE. debLevB ) THEN        IF ( debugLevel .GE. debLevB ) THEN
# Line 275  C--   Transfert the 2D-solution to "etaN Line 327  C--   Transfert the 2D-solution to "etaN
327        ENDDO        ENDDO
328    
329  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
330        IF ( nonHydrostatic ) THEN        IF ( use3Dsolver ) THEN
331    
332  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 ).
333  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
# Line 294  C     see CG3D.h for the interface to th Line 346  C     see CG3D.h for the interface to th
346           IF (useOBCS) THEN           IF (useOBCS) THEN
347            DO i=1,sNx+1            DO i=1,sNx+1
348  C Northern boundary  C Northern boundary
349            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(i,bi,bj).NE.0) THEN
350             vf(I,OB_Jn(I,bi,bj))=0.             vf(i,OB_Jn(i,bi,bj))=0.
351            ENDIF            ENDIF
352  C Southern boundary  C Southern boundary
353            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(i,bi,bj).NE.0) THEN
354             vf(I,OB_Js(I,bi,bj)+1)=0.             vf(i,OB_Js(i,bi,bj)+1)=0.
355            ENDIF            ENDIF
356            ENDDO            ENDDO
357            DO j=1,sNy+1            DO j=1,sNy+1
358  C Eastern boundary  C Eastern boundary
359            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(j,bi,bj).NE.0) THEN
360             uf(OB_Ie(J,bi,bj),J)=0.             uf(OB_Ie(j,bi,bj),j)=0.
361            ENDIF            ENDIF
362  C Western boundary  C Western boundary
363            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(j,bi,bj).NE.0) THEN
364             uf(OB_Iw(J,bi,bj)+1,J)=0.             uf(OB_Iw(j,bi,bj)+1,J)=0.
365            ENDIF            ENDIF
366            ENDDO            ENDDO
367           ENDIF           ENDIF
368  #endif  #endif /* ALLOW_OBCS */
369    
370           K=1           IF ( usingZCoords ) THEN
371    C-       Z coordinate: assume surface @ level k=1
372               tmpFac = freeSurfFac*deepFac2F(1)
373             ELSE
374    C-       Other than Z coordinate: no assumption on surface level index
375               tmpFac = 0.
376               DO j=1,sNy
377                DO i=1,sNx
378                  ks = ksurfC(i,j,bi,bj)
379                  IF ( ks.LE.Nr ) THEN
380                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
381         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
382         &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
383                  ENDIF
384                ENDDO
385               ENDDO
386             ENDIF
387             k=1
388             kp1 = MIN(k+1,Nr)
389             wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
390             IF (k.GE.Nr) wFacKp = 0.
391           DO j=1,sNy           DO j=1,sNy
392            DO i=1,sNx            DO i=1,sNx
393             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)
394       &       +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)
395       &       -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)
396       &       +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)
397       &       -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 )
398       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
399       &          -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
400       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
401            ENDDO            ENDDO
402           ENDDO           ENDDO
403           DO K=2,Nr-1           DO k=2,Nr
404              kp1 = MIN(k+1,Nr)
405    C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
406    C        both appear in wVel term, but at 2 different levels
407              wFacKm = deepFac2F( k )*rhoFacF( k )
408              wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
409              IF (k.GE.Nr) wFacKp = 0.
410            DO j=1,sNy            DO j=1,sNy
411             DO i=1,sNx             DO i=1,sNx
412              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)
413       &       +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)
414       &       -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)
415       &       +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)
416       &       -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 )
417       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
418       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*wFacKp
419       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
420    
421             ENDDO             ENDDO
422            ENDDO            ENDDO
423           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  
424    
425  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
426           IF (useOBCS) THEN           IF (useOBCS) THEN
427            DO K=1,Nr            DO k=1,Nr
428            DO i=1,sNx            DO i=1,sNx
429  C Northern boundary  C Northern boundary
430            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(i,bi,bj).NE.0) THEN
431             cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.             cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.
432            ENDIF            ENDIF
433  C Southern boundary  C Southern boundary
434            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(i,bi,bj).NE.0) THEN
435             cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.             cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.
436            ENDIF            ENDIF
437            ENDDO            ENDDO
438            DO j=1,sNy            DO j=1,sNy
439  C Eastern boundary  C Eastern boundary
440            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(j,bi,bj).NE.0) THEN
441             cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.             cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.
442            ENDIF            ENDIF
443  C Western boundary  C Western boundary
444            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(j,bi,bj).NE.0) THEN
445             cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.             cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.
446            ENDIF            ENDIF
447            ENDDO            ENDDO
448            ENDDO            ENDDO
449           ENDIF           ENDIF
450  #endif  #endif /* ALLOW_OBCS */
451    C-    end bi,bj loops
452          ENDDO ! bi          ENDDO
453         ENDDO ! bj         ENDDO
454    
455        firstResidual=0.        firstResidual=0.
456        lastResidual=0.        lastResidual=0.
457        numIters=cg2dMaxIters        numIters=cg3dMaxIters
458          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
459        CALL CG3D(        CALL CG3D(
460       U           cg3d_b,       U           cg3d_b,
461       U           phi_nh,       U           phi_nh,
# Line 398  C Western boundary Line 464  C Western boundary
464       U           numIters,       U           numIters,
465       I           myThid )       I           myThid )
466        _EXCH_XYZ_R8(phi_nh, myThid )        _EXCH_XYZ_R8(phi_nh, myThid )
467          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
468    
469        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
470       &   ) THEN       &   ) THEN
# Line 413  C Western boundary Line 480  C Western boundary
480         ENDIF         ENDIF
481        ENDIF        ENDIF
482    
483    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
484          IF ( zeroPsNH ) THEN
485           DO bj=myByLo(myThid),myByHi(myThid)
486            DO bi=myBxLo(myThid),myBxHi(myThid)
487    
488             IF ( usingZCoords ) THEN
489    C-       Z coordinate: assume surface @ level k=1
490              DO k=2,Nr
491               DO j=1-OLy,sNy+OLy
492                DO i=1-OLx,sNx+OLx
493                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
494         &                           - phi_nh(i,j,1,bi,bj)
495                ENDDO
496               ENDDO
497              ENDDO
498              DO j=1-OLy,sNy+OLy
499               DO i=1-OLx,sNx+OLx
500                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
501         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
502                 phi_nh(i,j,1,bi,bj) = 0.
503               ENDDO
504              ENDDO
505             ELSE
506    C-       Other than Z coordinate: no assumption on surface level index
507              DO j=1-OLy,sNy+OLy
508               DO i=1-OLx,sNx+OLx
509                ks = ksurfC(i,j,bi,bj)
510                IF ( ks.LE.Nr ) THEN
511                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
512         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
513                 DO k=Nr,1,-1
514                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
515         &                            - phi_nh(i,j,ks,bi,bj)
516                 ENDDO
517                ENDIF
518               ENDDO
519              ENDDO
520             ENDIF
521    
522            ENDDO
523           ENDDO
524        ENDIF        ENDIF
 #endif  
525    
 #ifdef TIME_PER_TIMESTEP  
 CCE107 Time per timestep information  
       _BEGIN_MASTER( myThid )  
       CALL TIMER_GET_TIME( utnew, stnew, wtnew )  
 C Only output timing information after the 1st timestep  
       IF ( wtold .NE. 0.0D0 ) THEN  
         WRITE(msgBuf,'(A34,3F10.6)')  
      $        'User, system and wallclock time:', utnew - utold,  
      $        stnew - stold, wtnew - wtold  
          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
526        ENDIF        ENDIF
527        utold = utnew  #endif /* ALLOW_NONHYDROSTATIC */
528        stold = stnew  
529        wtold = wtnew  #ifdef ALLOW_SHOWFLOPS
530        _END_MASTER( myThid )        CALL SHOWFLOPS_INSOLVE( myThid)
531  #endif  #endif
532    
533        RETURN        RETURN
534        END        END
   
 #ifdef TIME_PER_TIMESTEP  
 CCE107 Initialization of common block for per timestep timing  
       BLOCK DATA settimers  
 C     !TIMING VARIABLES  
 C     == Timing variables ==  
       REAL*8 utnew, utold, stnew, stold, wtnew, wtold  
       COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold  
       DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/  
       END  
 #endif  
Legend:
Removed from v.1.48  
changed lines
  Added in v.1.64
  ViewVC Help
Powered by ViewVC 1.1.22