/[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.4 by cnh, Mon Jun 8 21:43:02 1998 UTC revision 1.55 by ce107, Tue May 9 16:07:52 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #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)
 C     | o Controls inversion of two and/or three-dimensional     |  
 C     |   elliptic problems for the pressure field.              |  
 C     \==========================================================/  
11    
12    C     !DESCRIPTION: \bv
13    C     *==========================================================*
14    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"
26    #include "GRID.h"
27    #include "SURFACE.h"
28    #include "FFIELDS.h"
29  #include "DYNVARS.h"  #include "DYNVARS.h"
30  #include "CG2D.h"  #include "SOLVE_FOR_PRESSURE.h"
31    #ifdef ALLOW_NONHYDROSTATIC
32    #include "SOLVE_FOR_PRESSURE3D.h"
33    #include "NH_VARS.h"
34    #endif
35    #ifdef ALLOW_CD_CODE
36    #include "CD_CODE_VARS.h"
37    #endif
38    #ifdef ALLOW_OBCS
39    #include "OBCS.h"
40    #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        INTEGER i,j,bi,bj  C     == Local variables ==
57          INTEGER i,j,k,bi,bj
58          _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59          _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60          _RL firstResidual,lastResidual
61          _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  #ifdef ALLOW_CD  C--   Save previous solution & Initialise Vector solution and source term :
115  C--   Save previous solution.        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            cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)  #ifdef ALLOW_CD_CODE
121              etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
122    #endif
123              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.
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)
155           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
168             DO j=1,sNy+1
169              DO i=1,sNx+1
170               uf(i,j) = _dyG(i,j,bi,bj)
171         &      *drF(k)*_hFacW(i,j,k,bi,bj)
172               vf(i,j) = _dxG(i,j,bi,bj)
173         &      *drF(k)*_hFacS(i,j,k,bi,bj)
174              ENDDO
175             ENDDO
176             CALL CALC_DIV_GHAT(
177         I       bi,bj,1,sNx,1,sNy,K,
178         I       uf,vf,
179         U       cg2d_b,
180         I       myThid)
181            ENDDO
182           ENDDO
183          ENDDO
184    
185    C--   Add source term arising from w=d/dt (p_s + p_nh)
186          DO bj=myByLo(myThid),myByHi(myThid)
187           DO bi=myBxLo(myThid),myBxHi(myThid)
188    #ifdef ALLOW_NONHYDROSTATIC
189            IF ( use3Dsolver .AND. zeroPsNH ) 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
196         &         * etaH(i,j,bi,bj)
197                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
198         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
199         &         * etaH(i,j,bi,bj)
200               ENDIF
201              ENDDO
202             ENDDO
203            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
221            IF ( exactConserv ) THEN
222    #endif /* ALLOW_NONHYDROSTATIC */
223             DO j=1,sNy
224              DO i=1,sNx
225               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
226         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
227         &         * etaH(i,j,bi,bj)
228              ENDDO
229             ENDDO
230            ELSE
231             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
241            IF (useOBCS) THEN
242             DO i=1,sNx
243    C Northern boundary
244              IF (OB_Jn(I,bi,bj).NE.0) THEN
245               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
248    C Southern boundary
249              IF (OB_Js(I,bi,bj).NE.0) THEN
250               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
253             ENDDO
254             DO j=1,sNy
255    C Eastern boundary
256              IF (OB_Ie(J,bi,bj).NE.0) THEN
257               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
260    C Western boundary
261              IF (OB_Iw(J,bi,bj).NE.0) THEN
262               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
265             ENDDO
266            ENDIF
267    #endif /* ALLOW_OBCS */
268    C-    end bi,bj loops
269           ENDDO
270          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  #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.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        CALL CG2D(        CALL CG2D(
287         U           cg2d_b,
288         U           cg2d_x,
289         O           firstResidual,
290         O           lastResidual,
291         U           numIters,
292       I           myThid )       I           myThid )
293          _EXCH_XY_R8(cg2d_x, myThid )
294    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
295    
296    #ifdef ALLOW_DEBUG
297          IF ( debugLevel .GE. debLevB ) THEN
298           CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
299         &                        myThid)
300          ENDIF
301    #endif
302    
303    C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
304          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
305         &   ) THEN
306           IF ( debugLevel .GE. debLevA ) THEN
307            _BEGIN_MASTER( myThid )
308            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
309            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
310            WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
311            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
312            WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
313            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
314            _END_MASTER( myThid )
315           ENDIF
316          ENDIF
317    
318    C--   Transfert the 2D-solution to "etaN" :
319          DO bj=myByLo(myThid),myByHi(myThid)
320           DO bi=myBxLo(myThid),myBxHi(myThid)
321            DO j=1-OLy,sNy+OLy
322             DO i=1-OLx,sNx+OLx
323              etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
324             ENDDO
325            ENDDO
326           ENDDO
327          ENDDO
328    
329    #ifdef ALLOW_NONHYDROSTATIC
330          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.
334  C*P*  CALL CG3D(         DO bj=myByLo(myThid),myByHi(myThid)
335  C*P* I           myThid )          DO bi=myBxLo(myThid),myBxHi(myThid)
336             DO j=1,sNy+1
337              DO i=1,sNx+1
338               uf(i,j)=-_recip_dxC(i,j,bi,bj)*
339         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
340               vf(i,j)=-_recip_dyC(i,j,bi,bj)*
341         &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
342              ENDDO
343             ENDDO
344    
345    #ifdef ALLOW_OBCS
346             IF (useOBCS) THEN
347              DO i=1,sNx+1
348    C Northern boundary
349              IF (OB_Jn(I,bi,bj).NE.0) THEN
350               vf(I,OB_Jn(I,bi,bj))=0.
351              ENDIF
352    C Southern boundary
353              IF (OB_Js(I,bi,bj).NE.0) THEN
354               vf(I,OB_Js(I,bi,bj)+1)=0.
355              ENDIF
356              ENDDO
357              DO j=1,sNy+1
358    C Eastern boundary
359              IF (OB_Ie(J,bi,bj).NE.0) THEN
360               uf(OB_Ie(J,bi,bj),J)=0.
361              ENDIF
362    C Western boundary
363              IF (OB_Iw(J,bi,bj).NE.0) THEN
364               uf(OB_Iw(J,bi,bj)+1,J)=0.
365              ENDIF
366              ENDDO
367             ENDIF
368    #endif /* ALLOW_OBCS */
369    
370             IF ( usingZCoords ) THEN
371    C-       Z coordinate: assume surface @ level k=1
372               tmpFac = freeSurfFac
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)/deltaTmom
383                  ENDIF
384                ENDDO
385               ENDDO
386             ENDIF
387             K=1
388             kp1 = MIN(k+1,Nr)
389             maskKp1 = 1.
390             IF (k.GE.Nr) maskKp1 = 0.
391             DO j=1,sNy
392              DO i=1,sNx
393                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)
395         &       -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)
397         &       -drF(K)*dxG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
398         &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
399         &         -wVel(i,j,kp1,bi,bj)*maskKp1
400         &        )*_rA(i,j,bi,bj)/deltaTmom
401              ENDDO
402             ENDDO
403             DO K=2,Nr
404              kp1 = MIN(k+1,Nr)
405              maskKp1 = 1.
406              IF (k.GE.Nr) maskKp1 = 0.
407              DO j=1,sNy
408               DO i=1,sNx
409                cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
410         &       +drF(K)*dyG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
411         &       -drF(K)*dyG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
412         &       +drF(K)*dxG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
413         &       -drF(K)*dxG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
414         &       +( wVel(i,j,k  ,bi,bj)*maskC(i,j,k-1,bi,bj)
415         &         -wVel(i,j,kp1,bi,bj)*maskKp1
416         &        )*_rA(i,j,bi,bj)/deltaTmom
417    
418               ENDDO
419              ENDDO
420             ENDDO
421    
422    #ifdef ALLOW_OBCS
423             IF (useOBCS) THEN
424              DO K=1,Nr
425              DO i=1,sNx
426    C Northern boundary
427              IF (OB_Jn(I,bi,bj).NE.0) THEN
428               cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
429              ENDIF
430    C Southern boundary
431              IF (OB_Js(I,bi,bj).NE.0) THEN
432               cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
433              ENDIF
434              ENDDO
435              DO j=1,sNy
436    C Eastern boundary
437              IF (OB_Ie(J,bi,bj).NE.0) THEN
438               cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
439              ENDIF
440    C Western boundary
441              IF (OB_Iw(J,bi,bj).NE.0) THEN
442               cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
443              ENDIF
444              ENDDO
445              ENDDO
446             ENDIF
447    #endif /* ALLOW_OBCS */
448    C-    end bi,bj loops
449            ENDDO
450           ENDDO
451    
452          firstResidual=0.
453          lastResidual=0.
454          numIters=cg3dMaxIters
455          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
456          CALL CG3D(
457         U           cg3d_b,
458         U           phi_nh,
459         O           firstResidual,
460         O           lastResidual,
461         U           numIters,
462         I           myThid )
463          _EXCH_XYZ_R8(phi_nh, myThid )
464          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
465    
466          IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
467         &   ) THEN
468           IF ( debugLevel .GE. debLevA ) THEN
469            _BEGIN_MASTER( myThid )
470            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
471            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
472            WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
473            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
474            WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
475            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
476            _END_MASTER( myThid )
477           ENDIF
478          ENDIF
479    
480    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
481          IF ( zeroPsNH ) THEN
482           DO bj=myByLo(myThid),myByHi(myThid)
483            DO bi=myBxLo(myThid),myBxHi(myThid)
484    
485             IF ( usingZCoords ) THEN
486    C-       Z coordinate: assume surface @ level k=1
487              DO k=2,Nr
488               DO j=1-OLy,sNy+OLy
489                DO i=1-OLx,sNx+OLx
490                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
491         &                           - phi_nh(i,j,1,bi,bj)
492                ENDDO
493               ENDDO
494              ENDDO
495              DO j=1-OLy,sNy+OLy
496               DO i=1-OLx,sNx+OLx
497                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
498         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
499                 phi_nh(i,j,1,bi,bj) = 0.
500               ENDDO
501              ENDDO
502             ELSE
503    C-       Other than Z coordinate: no assumption on surface level index
504              DO j=1-OLy,sNy+OLy
505               DO i=1-OLx,sNx+OLx
506                ks = ksurfC(i,j,bi,bj)
507                IF ( ks.LE.Nr ) THEN
508                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
509         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
510                 DO k=Nr,1,-1
511                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
512         &                            - phi_nh(i,j,ks,bi,bj)
513                 ENDDO
514                ENDIF
515               ENDDO
516              ENDDO
517             ENDIF
518    
519            ENDDO
520           ENDDO
521          ENDIF
522    
523          ENDIF
524    #endif /* ALLOW_NONHYDROSTATIC */
525    
526    #ifdef TIME_PER_TIMESTEP_SFP
527    CCE107 Time per timestep information
528          _BEGIN_MASTER( myThid )
529          CALL TIMER_GET_TIME( utnew, stnew, wtnew )
530    C Only output timing information after the 1st timestep
531          IF ( wtold .NE. 0.0D0 ) THEN
532            WRITE(msgBuf,'(A34,3F10.6)')
533         $        'User, system and wallclock time:', utnew - utold,
534         $        stnew - stold, wtnew - wtold
535             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
536          ENDIF
537          utold = utnew
538          stold = stnew
539          wtold = wtnew
540          _END_MASTER( myThid )
541    #endif
542    #ifdef USE_PAPI_FLOPS_SFP
543    CCE107 PAPI summary performance
544          _BEGIN_MASTER( myThid )
545    #ifdef USE_FLIPS
546          call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
547    #else
548          call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
549    #endif
550          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
551         $     'Mflop/s during this timestep:', mflops, ' ', mflops
552         $     *proc_time/(real_time + 1E-36)
553          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
554    #ifdef PAPI_VERSION
555          call PAPIF_ipc(real_time, proc_time, instr, ipc, check)
556          WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
557         $     'IPC during this timestep:', ipc, ' ', ipc*proc_time
558         $     /(real_time + 1E-36),
559          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
560    #endif
561          _END_MASTER( myThid )
562    #else
563    #ifdef USE_PCL_FLOPS_SFP
564    CCE107 PCL summary performance
565          _BEGIN_MASTER( myThid )
566          PCLstop(descr, i_result, fp_result, nevents)
567          do ipcl = 1, nevents
568             WRITE(msgBuf,'(A22,A26,F10.6)'),
569         $        pcl_counter_name(pcl_counter_list(ipcl)),
570         $        'during this timestep:', fp_results(ipcl)
571             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
572          enddo
573          PCLstart(descr, pcl_counter_list, nevents, flags)
574          _END_MASTER( myThid )
575    #endif
576    #endif
577        RETURN        RETURN
578        END        END
579    
580    #ifdef TIME_PER_TIMESTEP_SFP
581    CCE107 Initialization of common block for per timestep timing
582          BLOCK DATA settimers
583    C     !TIMING VARIABLES
584    C     == Timing variables ==
585          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
586          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
587          DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
588          END
589    #endif
590    #ifdef USE_PAPI_FLOPS_SFP
591    CCE107 Initialization of common block for PAPI summary performance
592          BLOCK DATA setpapis
593          INTEGER*8 flpops, instr
594          REAL real_time, proc_time, mflops, ipc
595          COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
596          DATA flpops, instr, real_time, proc_time, mflops, ipc /2*0,4*0.E0/
597          END
598    #endif
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.55
  ViewVC Help
Powered by ViewVC 1.1.22