/[MITgcm]/MITgcm/model/src/solve_for_pressure.F
ViewVC logotype

Diff of /MITgcm/model/src/solve_for_pressure.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.37 by mlosch, Tue May 13 13:30:05 2003 UTC revision 1.52 by ce107, Thu Dec 22 01:08:57 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 22  C     == Global variables Line 23  C     == Global variables
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
26  #include "GRID.h"  #include "GRID.h"
27  #include "SURFACE.h"  #include "SURFACE.h"
28  #include "FFIELDS.h"  #include "FFIELDS.h"
29    #include "DYNVARS.h"
30    #include "SOLVE_FOR_PRESSURE.h"
31  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
32  #include "SOLVE_FOR_PRESSURE3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
33  #include "GW.h"  #include "NH_VARS.h"
34    #endif
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37  #endif  #endif
38  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
39  #include "OBCS.h"  #include "OBCS.h"
40  #endif  #endif
 #include "SOLVE_FOR_PRESSURE.h"  
41    
42  C     === Functions ====  C     === Functions ====
43        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
# Line 55  C     == Local variables == Line 59  C     == Local variables ==
59        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL firstResidual,lastResidual        _RL firstResidual,lastResidual
61        _RL tmpFac        _RL tmpFac
62          _RL sumEmP, tileEmP
63          LOGICAL putPmEinXvector
64        INTEGER numIters        INTEGER numIters
65        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
66    #ifdef ALLOW_NONHYDROSTATIC
67          INTEGER ks, kp1
68          _RL     maskKp1
69          LOGICAL zeroPsNH
70    #endif
71  CEOP  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
84          INTEGER check
85          REAL real_time, proc_time, mflops
86          COMMON /papivars/ flpops, real_time, proc_time, mflops, check
87    #endif
88    
89    #ifdef ALLOW_NONHYDROSTATIC
90    c       zeroPsNH = .FALSE.
91            zeroPsNH = exactConserv
92    #endif
93    
94    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
95    C     instead of simply etaN ; This can speed-up the solver convergence in
96    C     the case where |Global_mean_PmE| is large.
97          putPmEinXvector = .FALSE.
98    c     putPmEinXvector = useRealFreshWaterFlux
99    
100  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
101          sumEmP = 0.
102        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
103         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
104          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
105           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
106  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
107            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
108  #endif  #endif
109            cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)            cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
# Line 82  C--   Save previous solution & Initialis Line 121  C--   Save previous solution & Initialis
121            ENDDO            ENDDO
122           ENDDO           ENDDO
123          ENDIF          ENDIF
124            IF ( putPmEinXvector ) THEN
125             tileEmP = 0.
126             DO j=1,sNy
127              DO i=1,sNx
128                tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
129         &                                       *maskH(i,j,bi,bj)
130              ENDDO
131             ENDDO
132             sumEmP = sumEmP + tileEmP
133            ENDIF
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136          IF ( putPmEinXvector ) THEN
137            _GLOBAL_SUM_R8( sumEmP, myThid )
138          ENDIF
139    
140        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
141         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
142            IF ( putPmEinXvector ) THEN
143              tmpFac = 0.
144              IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
145         &                          *convertEmP2rUnit*sumEmP/globalArea
146              DO j=1,sNy
147               DO i=1,sNx
148                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
149         &                        - tmpFac*Bo_surf(i,j,bi,bj)
150               ENDDO
151              ENDDO
152            ENDIF
153          DO K=Nr,1,-1          DO K=Nr,1,-1
154           DO j=1,sNy+1           DO j=1,sNy+1
155            DO i=1,sNx+1            DO i=1,sNx+1
# Line 109  C--   Add source term arising from w=d/d Line 172  C--   Add source term arising from w=d/d
172        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
173         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
174  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
175          IF ( nonHydrostatic ) THEN          IF ( nonHydrostatic .AND. zeroPsNH ) THEN
176           DO j=1,sNy           DO j=1,sNy
177            DO i=1,sNx            DO i=1,sNx
178             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
179               IF ( ks.LE.Nr ) THEN
180                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
181         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
182         &         * etaH(i,j,bi,bj)
183                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
184         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
185         &         * etaH(i,j,bi,bj)
186               ENDIF
187              ENDDO
188             ENDDO
189            ELSEIF ( nonHydrostatic ) THEN
190             DO j=1,sNy
191              DO i=1,sNx
192               ks = ksurfC(i,j,bi,bj)
193               IF ( ks.LE.Nr ) THEN
194                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
195       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
196       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
197       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
198             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)
199       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
200       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
201       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
202               ENDIF
203            ENDDO            ENDDO
204           ENDDO           ENDDO
205          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
206  #else  #else
207          IF ( exactConserv ) THEN          IF ( exactConserv ) THEN
208  #endif  #endif /* ALLOW_NONHYDROSTATIC */
209           DO j=1,sNy           DO j=1,sNy
210            DO i=1,sNx            DO i=1,sNx
211             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
# Line 170  C Western boundary Line 250  C Western boundary
250            ENDIF            ENDIF
251           ENDDO           ENDDO
252          ENDIF          ENDIF
253  #endif  #endif /* ALLOW_OBCS */
254    C-    end bi,bj loops
255         ENDDO         ENDDO
256        ENDDO        ENDDO
257    
258  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
259        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
260         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
261       &                        myThid)       &                        myThid)
262        ENDIF        ENDIF
# Line 187  C     see CG2D.h for the interface to th Line 268  C     see CG2D.h for the interface to th
268        firstResidual=0.        firstResidual=0.
269        lastResidual=0.        lastResidual=0.
270        numIters=cg2dMaxIters        numIters=cg2dMaxIters
271    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
272        CALL CG2D(        CALL CG2D(
273       U           cg2d_b,       U           cg2d_b,
274       U           cg2d_x,       U           cg2d_x,
# Line 195  C     see CG2D.h for the interface to th Line 277  C     see CG2D.h for the interface to th
277       U           numIters,       U           numIters,
278       I           myThid )       I           myThid )
279        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
280    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
281    
282  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
283        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
284         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
285       &                        myThid)       &                        myThid)
286        ENDIF        ENDIF
287  #endif  #endif
288    
289  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :  C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
290        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
291       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
292         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
293         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
294         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
295         WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
296         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
297         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
298         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
299         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
300            _END_MASTER( myThid )
301           ENDIF
302        ENDIF        ENDIF
303    
304  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
# Line 266  C Western boundary Line 351  C Western boundary
351            ENDIF            ENDIF
352            ENDDO            ENDDO
353           ENDIF           ENDIF
354  #endif  #endif /* ALLOW_OBCS */
355    
356             IF ( usingZCoords ) THEN
357    C-       Z coordinate: assume surface @ level k=1
358               tmpFac = freeSurfFac
359             ELSE
360    C-       Other than Z coordinate: no assumption on surface level index
361               tmpFac = 0.
362               DO j=1,sNy
363                DO i=1,sNx
364                  ks = ksurfC(i,j,bi,bj)
365                  IF ( ks.LE.Nr ) THEN
366                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
367         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
368         &                          *_rA(i,j,bi,bj)/deltaTmom
369                  ENDIF
370                ENDDO
371               ENDDO
372             ENDIF
373           K=1           K=1
374             kp1 = MIN(k+1,Nr)
375             maskKp1 = 1.
376             IF (k.GE.Nr) maskKp1 = 0.
377           DO j=1,sNy           DO j=1,sNy
378            DO i=1,sNx            DO i=1,sNx
379             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)
380       &       +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)
381       &       -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)
382       &       +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)
383       &       -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 )
384       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
385       &          -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
386       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
387            ENDDO            ENDDO
388           ENDDO           ENDDO
389           DO K=2,Nr-1           DO K=2,Nr
390              kp1 = MIN(k+1,Nr)
391              maskKp1 = 1.
392              IF (k.GE.Nr) maskKp1 = 0.
393            DO j=1,sNy            DO j=1,sNy
394             DO i=1,sNx             DO i=1,sNx
395              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)
396       &       +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)
397       &       -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)
398       &       +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)
399       &       -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 )
400       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j,k  ,bi,bj)*maskC(i,j,k-1,bi,bj)
401       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
402       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
403    
404             ENDDO             ENDDO
405            ENDDO            ENDDO
406           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  
407    
408  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
409           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 335  C Western boundary Line 430  C Western boundary
430            ENDDO            ENDDO
431            ENDDO            ENDDO
432           ENDIF           ENDIF
433  #endif  #endif /* ALLOW_OBCS */
434    C-    end bi,bj loops
435          ENDDO ! bi          ENDDO
436         ENDDO ! bj         ENDDO
437    
438        firstResidual=0.        firstResidual=0.
439        lastResidual=0.        lastResidual=0.
440        numIters=cg2dMaxIters        numIters=cg3dMaxIters
441          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
442        CALL CG3D(        CALL CG3D(
443       U           cg3d_b,       U           cg3d_b,
444       U           phi_nh,       U           phi_nh,
# Line 351  C Western boundary Line 447  C Western boundary
447       U           numIters,       U           numIters,
448       I           myThid )       I           myThid )
449        _EXCH_XYZ_R8(phi_nh, myThid )        _EXCH_XYZ_R8(phi_nh, myThid )
450          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
451    
452          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
453         &   ) THEN
454           IF ( debugLevel .GE. debLevA ) THEN
455            _BEGIN_MASTER( myThid )
456            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
457            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
458            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
459            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
460            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
461            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
462            _END_MASTER( myThid )
463           ENDIF
464          ENDIF
465    
466    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
467          IF ( zeroPsNH ) THEN
468           DO bj=myByLo(myThid),myByHi(myThid)
469            DO bi=myBxLo(myThid),myBxHi(myThid)
470    
471        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,           IF ( usingZCoords ) THEN
472       &                                    myTime-deltaTClock) ) THEN  C-       Z coordinate: assume surface @ level k=1
473         _BEGIN_MASTER( myThid )            DO k=2,Nr
474         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual             DO j=1-OLy,sNy+OLy
475         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)              DO i=1-OLx,sNx+OLx
476         WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters               phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
477         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &                           - phi_nh(i,j,1,bi,bj)
478         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual              ENDDO
479         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)             ENDDO
480         _END_MASTER( )            ENDDO
481              DO j=1-OLy,sNy+OLy
482               DO i=1-OLx,sNx+OLx
483                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
484         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
485                 phi_nh(i,j,1,bi,bj) = 0.
486               ENDDO
487              ENDDO
488             ELSE
489    C-       Other than Z coordinate: no assumption on surface level index
490              DO j=1-OLy,sNy+OLy
491               DO i=1-OLx,sNx+OLx
492                ks = ksurfC(i,j,bi,bj)
493                IF ( ks.LE.Nr ) THEN
494                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
495         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
496                 DO k=Nr,1,-1
497                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
498         &                            - phi_nh(i,j,ks,bi,bj)
499                 ENDDO
500                ENDIF
501               ENDDO
502              ENDDO
503             ENDIF
504    
505            ENDDO
506           ENDDO
507        ENDIF        ENDIF
508    
509        ENDIF        ENDIF
510  #endif  #endif /* ALLOW_NONHYDROSTATIC */
511    
512    #ifdef TIME_PER_TIMESTEP_SFP
513    CCE107 Time per timestep information
514          _BEGIN_MASTER( myThid )
515          CALL TIMER_GET_TIME( utnew, stnew, wtnew )
516    C Only output timing information after the 1st timestep
517          IF ( wtold .NE. 0.0D0 ) THEN
518            WRITE(msgBuf,'(A34,3F10.6)')
519         $        'User, system and wallclock time:', utnew - utold,
520         $        stnew - stold, wtnew - wtold
521             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
522          ENDIF
523          utold = utnew
524          stold = stnew
525          wtold = wtnew
526          _END_MASTER( myThid )
527    #endif
528    #ifdef USE_PAPI_FLOPS_SFP
529    CCE107 PAPI summary performance
530          _BEGIN_MASTER( myThid )
531          call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
532          WRITE(msgBuf,'(A34,F10.6)')
533         $        'Mflop/s during this timestep:', mflops
534          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
535          _END_MASTER( myThid )
536    #endif
537        RETURN        RETURN
538        END        END
539    
540    #ifdef TIME_PER_TIMESTEP_SFP
541    CCE107 Initialization of common block for per timestep timing
542          BLOCK DATA settimers
543    C     !TIMING VARIABLES
544    C     == Timing variables ==
545          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
546          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
547          DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
548          END
549    #endif
550    #ifdef USE_PAPI_FLOPS_SFP
551    CCE107 Initialization of common block for PAPI summary performance
552          BLOCK DATA setpapis
553          INTEGER*8 flpops
554          INTEGER check
555          REAL real_time, proc_time, mflops
556          COMMON /papivars/ flpops, real_time, proc_time, mflops, check
557          DATA flpops, real_time, proc_time, mflops, check /0, 3*0.0E0, 0/
558          END
559    #endif
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22