/[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.56 by heimbach, Wed Jun 7 01:55:13 2006 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)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL firstResidual,lastResidual
61        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL tmpFac
62          _RL sumEmP, tileEmP
63          LOGICAL putPmEinXvector
64          INTEGER numIters
65          CHARACTER*(MAX_LEN_MBUF) msgBuf
66    #ifdef ALLOW_NONHYDROSTATIC
67          INTEGER ks, kp1
68          _RL     maskKp1
69          LOGICAL zeroPsNH
70    #endif
71    CEOP
72    
73    #ifdef TIME_PER_TIMESTEP_SFP
74    CCE107 common block for per timestep timing
75    C     !TIMING VARIABLES
76    C     == Timing variables ==
77          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
78          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
79    #endif
80    #ifdef USE_PAPI_FLOPS_SFP
81    CCE107 common block for PAPI summary performance
82    #include <fpapi.h>
83          INTEGER*8 flpops, instr
84          INTEGER check
85          REAL*4 real_time, proc_time, mflops, ipc
86          COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
87    #else
88    #ifdef USE_PCL_FLOPS_SFP
89    CCE107 common block for PCL summary performance
90    #include <pclh.f>
91          INTEGER pcl_counter_list(5), flags, nevents, res, ipcl
92          INTEGER*8 i_result(5), descr
93          REAL*8 fp_result(5)
94          COMMON /pclvars/ i_result, descr, fp_result, pcl_counter_list,
95         $     flags, nevents
96          INTEGER nmaxevents
97          PARAMETER (nmaxevents = 61)
98          CHARACTER*22 pcl_counter_name(0:nmaxevents-1)
99          COMMON /pclnames/ pcl_counter_name
100    #endif
101    #endif
102    
103    #ifdef ALLOW_NONHYDROSTATIC
104    c       zeroPsNH = .FALSE.
105            zeroPsNH = exactConserv
106    #endif
107    
108    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
109    C     instead of simply etaN ; This can speed-up the solver convergence in
110    C     the case where |Global_mean_PmE| is large.
111          putPmEinXvector = .FALSE.
112    c     putPmEinXvector = useRealFreshWaterFlux
113    
114  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
115          sumEmP = 0.
116        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
117         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
118          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
119           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
120  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
121            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
122  #endif  #endif
123            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)
124            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  
125           ENDDO           ENDDO
126          ENDDO          ENDDO
127            IF (useRealFreshWaterFlux) THEN
128             tmpFac = freeSurfFac*convertEmP2rUnit
129             IF (exactConserv)
130         &        tmpFac = freeSurfFac*convertEmP2rUnit*implicDiv2DFlow
131             DO j=1,sNy
132              DO i=1,sNx
133               cg2d_b(i,j,bi,bj) =
134         &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
135              ENDDO
136             ENDDO
137            ENDIF
138            IF ( putPmEinXvector ) THEN
139             tileEmP = 0.
140             DO j=1,sNy
141              DO i=1,sNx
142                tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
143         &                                       *maskH(i,j,bi,bj)
144              ENDDO
145             ENDDO
146             sumEmP = sumEmP + tileEmP
147            ENDIF
148         ENDDO         ENDDO
149        ENDDO        ENDDO
150          IF ( putPmEinXvector ) THEN
151            _GLOBAL_SUM_R8( sumEmP, myThid )
152          ENDIF
153    
154        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
155         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
156            IF ( putPmEinXvector ) THEN
157              tmpFac = 0.
158              IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
159         &                          *convertEmP2rUnit*sumEmP/globalArea
160              DO j=1,sNy
161               DO i=1,sNx
162                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
163         &                        - tmpFac*Bo_surf(i,j,bi,bj)
164               ENDDO
165              ENDDO
166            ENDIF
167          DO K=Nr,1,-1          DO K=Nr,1,-1
168           DO j=1,sNy+1           DO j=1,sNy+1
169            DO i=1,sNx+1            DO i=1,sNx+1
# Line 84  C--   Add source term arising from w=d/d Line 186  C--   Add source term arising from w=d/d
186        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
187         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
188  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
189          DO j=1,sNy          IF ( use3Dsolver .AND. zeroPsNH ) THEN
190           DO i=1,sNx           DO j=1,sNy
191            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            DO i=1,sNx
192       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(             ks = ksurfC(i,j,bi,bj)
193       &          -cg2d_x(I,J,bi,bj)             IF ( ks.LE.Nr ) THEN
194       &          -cg3d_x(I,J,1,bi,bj)              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
195       &        )/deltaTMom/deltaTMom       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
196            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)       &         * etaH(i,j,bi,bj)
197       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(              cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
198       &         -cg2d_x(I,J,bi,bj)       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
199       &         -cg3d_x(I,J,1,bi,bj)       &         * etaH(i,j,bi,bj)
200       &       )/deltaTMom/deltaTMom             ENDIF
201              ENDDO
202           ENDDO           ENDDO
203          ENDDO          ELSEIF ( use3Dsolver ) THEN
204             DO j=1,sNy
205              DO i=1,sNx
206               ks = ksurfC(i,j,bi,bj)
207               IF ( ks.LE.Nr ) THEN
208                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
209         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
210         &         *( etaN(i,j,bi,bj)
211         &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
212                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
213         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
214         &         *( etaN(i,j,bi,bj)
215         &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
216               ENDIF
217              ENDDO
218             ENDDO
219            ELSEIF ( exactConserv ) THEN
220  #else  #else
221          DO j=1,sNy          IF ( exactConserv ) THEN
222           DO i=1,sNx  #endif /* ALLOW_NONHYDROSTATIC */
223            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)           DO j=1,sNy
224       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(            DO i=1,sNx
225       &          -cg2d_x(I,J,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
226       &        )/deltaTMom/deltaTMom       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
227         &         * etaH(i,j,bi,bj)
228              ENDDO
229           ENDDO           ENDDO
230          ENDDO          ELSE
231  #endif           DO j=1,sNy
232              DO i=1,sNx
233               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
234         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
235         &         * etaN(i,j,bi,bj)
236              ENDDO
237             ENDDO
238            ENDIF
239    
240  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
241          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 115  C--   Add source term arising from w=d/d Line 243  C--   Add source term arising from w=d/d
243  C Northern boundary  C Northern boundary
244            IF (OB_Jn(I,bi,bj).NE.0) THEN            IF (OB_Jn(I,bi,bj).NE.0) THEN
245             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
246               cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
247            ENDIF            ENDIF
248  C Southern boundary  C Southern boundary
249            IF (OB_Js(I,bi,bj).NE.0) THEN            IF (OB_Js(I,bi,bj).NE.0) THEN
250             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.             cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
251               cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
252            ENDIF            ENDIF
253           ENDDO           ENDDO
254           DO j=1,sNy           DO j=1,sNy
255  C Eastern boundary  C Eastern boundary
256            IF (OB_Ie(J,bi,bj).NE.0) THEN            IF (OB_Ie(J,bi,bj).NE.0) THEN
257             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
258               cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
259            ENDIF            ENDIF
260  C Western boundary  C Western boundary
261            IF (OB_Iw(J,bi,bj).NE.0) THEN            IF (OB_Iw(J,bi,bj).NE.0) THEN
262             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.             cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
263               cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
264            ENDIF            ENDIF
265           ENDDO           ENDDO
266          ENDIF          ENDIF
267  #endif  #endif /* ALLOW_OBCS */
268    C-    end bi,bj loops
269         ENDDO         ENDDO
270        ENDDO        ENDDO
271    
272    #ifdef ALLOW_DEBUG
273          IF ( debugLevel .GE. debLevB ) THEN
274           CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
275         &                        myThid)
276          ENDIF
277    #endif
278    
279  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
280  C--   gradient solver.  C--   gradient solver.
281  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
282          firstResidual=0.
283          lastResidual=0.
284          numIters=cg2dMaxIters
285    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
286    #ifdef ALLOW_CG2D_NSA
287    C--   Call the not-self-adjoint version of cg2d
288          CALL CG2D_NSA(
289         U           cg2d_b,
290         U           cg2d_x,
291         O           firstResidual,
292         O           lastResidual,
293         U           numIters,
294         I           myThid )
295    #else /* not ALLOW_CG2D_NSA = default */
296        CALL CG2D(        CALL CG2D(
297       I           cg2d_b,       U           cg2d_b,
298       U           cg2d_x,       U           cg2d_x,
299         O           firstResidual,
300         O           lastResidual,
301         U           numIters,
302       I           myThid )       I           myThid )
303    #endif /* ALLOW_CG2D_NSA */
304        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
305    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
306    
307    #ifdef ALLOW_DEBUG
308          IF ( debugLevel .GE. debLevB ) THEN
309           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
310         &                        myThid)
311          ENDIF
312    #endif
313    
314    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
315          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
316         &   ) THEN
317           IF ( debugLevel .GE. debLevA ) THEN
318            _BEGIN_MASTER( myThid )
319            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
320            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
321            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
322            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
323            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
324            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
325            _END_MASTER( myThid )
326           ENDIF
327          ENDIF
328    
329  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
330        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
331         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
332          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
333           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
334            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)
335           ENDDO           ENDDO
336          ENDDO          ENDDO
337         ENDDO         ENDDO
338        ENDDO        ENDDO
339    
340  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
341        IF ( nonHydrostatic ) THEN        IF ( use3Dsolver ) THEN
342    
343  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 ).
344  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 346  C     see CG3D.h for the interface to th
346          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
347           DO j=1,sNy+1           DO j=1,sNy+1
348            DO i=1,sNx+1            DO i=1,sNx+1
349             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
350       &         (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))
351             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
352       &         (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))
353            ENDDO            ENDDO
354           ENDDO           ENDDO
# Line 197  C Western boundary Line 376  C Western boundary
376            ENDIF            ENDIF
377            ENDDO            ENDDO
378           ENDIF           ENDIF
379  #endif  #endif /* ALLOW_OBCS */
380    
381             IF ( usingZCoords ) THEN
382    C-       Z coordinate: assume surface @ level k=1
383               tmpFac = freeSurfFac
384             ELSE
385    C-       Other than Z coordinate: no assumption on surface level index
386               tmpFac = 0.
387               DO j=1,sNy
388                DO i=1,sNx
389                  ks = ksurfC(i,j,bi,bj)
390                  IF ( ks.LE.Nr ) THEN
391                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
392         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
393         &                          *_rA(i,j,bi,bj)/deltaTmom
394                  ENDIF
395                ENDDO
396               ENDDO
397             ENDIF
398           K=1           K=1
399             kp1 = MIN(k+1,Nr)
400             maskKp1 = 1.
401             IF (k.GE.Nr) maskKp1 = 0.
402           DO j=1,sNy           DO j=1,sNy
403            DO i=1,sNx            DO i=1,sNx
404             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)
405       &       +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)
406       &       -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)
407       &       +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)
408       &       -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 )
409       &       +(       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
410       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
411       &        )*_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  
412            ENDDO            ENDDO
413           ENDDO           ENDDO
414           DO K=2,Nr-1           DO K=2,Nr
415              kp1 = MIN(k+1,Nr)
416              maskKp1 = 1.
417              IF (k.GE.Nr) maskKp1 = 0.
418            DO j=1,sNy            DO j=1,sNy
419             DO i=1,sNx             DO i=1,sNx
420              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)
421       &       +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)
422       &       -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)
423       &       +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)
424       &       -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 )
425       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j,k  ,bi,bj)*maskC(i,j,k-1,bi,bj)
426       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
427       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
428    
429             ENDDO             ENDDO
430            ENDDO            ENDDO
431           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  
432    
433  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
434           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 269  C Western boundary Line 455  C Western boundary
455            ENDDO            ENDDO
456            ENDDO            ENDDO
457           ENDIF           ENDIF
458  #endif  #endif /* ALLOW_OBCS */
459    C-    end bi,bj loops
460            ENDDO
461           ENDDO
462    
463          ENDDO ! bi        firstResidual=0.
464         ENDDO ! bj        lastResidual=0.
465          numIters=cg3dMaxIters
466          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
467          CALL CG3D(
468         U           cg3d_b,
469         U           phi_nh,
470         O           firstResidual,
471         O           lastResidual,
472         U           numIters,
473         I           myThid )
474          _EXCH_XYZ_R8(phi_nh, myThid )
475          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
476    
477         CALL CG3D( myThid )        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
478         _EXCH_XYZ_R8(cg3d_x, myThid )       &   ) THEN
479           IF ( debugLevel .GE. debLevA ) THEN
480            _BEGIN_MASTER( myThid )
481            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
482            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
483            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
484            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
485            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
486            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
487            _END_MASTER( myThid )
488           ENDIF
489          ENDIF
490    
491    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
492          IF ( zeroPsNH ) THEN
493           DO bj=myByLo(myThid),myByHi(myThid)
494            DO bi=myBxLo(myThid),myBxHi(myThid)
495    
496             IF ( usingZCoords ) THEN
497    C-       Z coordinate: assume surface @ level k=1
498              DO k=2,Nr
499               DO j=1-OLy,sNy+OLy
500                DO i=1-OLx,sNx+OLx
501                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
502         &                           - phi_nh(i,j,1,bi,bj)
503                ENDDO
504               ENDDO
505              ENDDO
506              DO j=1-OLy,sNy+OLy
507               DO i=1-OLx,sNx+OLx
508                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
509         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
510                 phi_nh(i,j,1,bi,bj) = 0.
511               ENDDO
512              ENDDO
513             ELSE
514    C-       Other than Z coordinate: no assumption on surface level index
515              DO j=1-OLy,sNy+OLy
516               DO i=1-OLx,sNx+OLx
517                ks = ksurfC(i,j,bi,bj)
518                IF ( ks.LE.Nr ) THEN
519                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
520         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
521                 DO k=Nr,1,-1
522                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
523         &                            - phi_nh(i,j,ks,bi,bj)
524                 ENDDO
525                ENDIF
526               ENDDO
527              ENDDO
528             ENDIF
529    
530            ENDDO
531           ENDDO
532        ENDIF        ENDIF
 #endif  
533    
534          ENDIF
535    #endif /* ALLOW_NONHYDROSTATIC */
536    
537    #ifdef TIME_PER_TIMESTEP_SFP
538    CCE107 Time per timestep information
539          _BEGIN_MASTER( myThid )
540          CALL TIMER_GET_TIME( utnew, stnew, wtnew )
541    C Only output timing information after the 1st timestep
542          IF ( wtold .NE. 0.0D0 ) THEN
543            WRITE(msgBuf,'(A34,3F10.6)')
544         $        'User, system and wallclock time:', utnew - utold,
545         $        stnew - stold, wtnew - wtold
546             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
547          ENDIF
548          utold = utnew
549          stold = stnew
550          wtold = wtnew
551          _END_MASTER( myThid )
552    #endif
553    #ifdef USE_PAPI_FLOPS_SFP
554    CCE107 PAPI summary performance
555          _BEGIN_MASTER( myThid )
556    #ifdef USE_FLIPS
557          call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
558    #else
559          call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
560    #endif
561          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
562         $     'Mflop/s during this timestep:', mflops, ' ', mflops
563         $     *proc_time/(real_time + 1E-36)
564          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
565    #ifdef PAPI_VERSION
566          call PAPIF_ipc(real_time, proc_time, instr, ipc, check)
567          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
568         $     'IPC during this timestep:', ipc, ' ', ipc*proc_time
569         $     /(real_time + 1E-36),
570          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
571    #endif
572          _END_MASTER( myThid )
573    #else
574    #ifdef USE_PCL_FLOPS_SFP
575    CCE107 PCL summary performance
576          _BEGIN_MASTER( myThid )
577          PCLstop(descr, i_result, fp_result, nevents)
578          do ipcl = 1, nevents
579             WRITE(msgBuf,'(A22,A26,F10.6)'),
580         $        pcl_counter_name(pcl_counter_list(ipcl)),
581         $        'during this timestep:', fp_results(ipcl)
582             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
583          enddo
584          PCLstart(descr, pcl_counter_list, nevents, flags)
585          _END_MASTER( myThid )
586    #endif
587    #endif
588        RETURN        RETURN
589        END        END
590    
591    #ifdef TIME_PER_TIMESTEP_SFP
592    CCE107 Initialization of common block for per timestep timing
593          BLOCK DATA settimers
594    C     !TIMING VARIABLES
595    C     == Timing variables ==
596          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
597          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
598          DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
599          END
600    #endif
601    #ifdef USE_PAPI_FLOPS_SFP
602    CCE107 Initialization of common block for PAPI summary performance
603          BLOCK DATA setpapis
604          INTEGER*8 flpops, instr
605          REAL real_time, proc_time, mflops, ipc
606          COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
607          DATA flpops, instr, real_time, proc_time, mflops, ipc /2*0,4*0.E0/
608          END
609    #endif
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.56
  ViewVC Help
Powered by ViewVC 1.1.22