/[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.51 by jmc, Tue Dec 20 20:31:28 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
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    
81    #ifdef ALLOW_NONHYDROSTATIC
82    c       zeroPsNH = .FALSE.
83            zeroPsNH = exactConserv
84    #endif
85    
86    C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
87    C     instead of simply etaN ; This can speed-up the solver convergence in
88    C     the case where |Global_mean_PmE| is large.
89          putPmEinXvector = .FALSE.
90    c     putPmEinXvector = useRealFreshWaterFlux
91    
92  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
93          sumEmP = 0.
94        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
95         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
96          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
97           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
98  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
99            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
100  #endif  #endif
101            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 113  C--   Save previous solution & Initialis
113            ENDDO            ENDDO
114           ENDDO           ENDDO
115          ENDIF          ENDIF
116            IF ( putPmEinXvector ) THEN
117             tileEmP = 0.
118             DO j=1,sNy
119              DO i=1,sNx
120                tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
121         &                                       *maskH(i,j,bi,bj)
122              ENDDO
123             ENDDO
124             sumEmP = sumEmP + tileEmP
125            ENDIF
126         ENDDO         ENDDO
127        ENDDO        ENDDO
128          IF ( putPmEinXvector ) THEN
129            _GLOBAL_SUM_R8( sumEmP, myThid )
130          ENDIF
131    
132        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
133         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
134            IF ( putPmEinXvector ) THEN
135              tmpFac = 0.
136              IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
137         &                          *convertEmP2rUnit*sumEmP/globalArea
138              DO j=1,sNy
139               DO i=1,sNx
140                cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
141         &                        - tmpFac*Bo_surf(i,j,bi,bj)
142               ENDDO
143              ENDDO
144            ENDIF
145          DO K=Nr,1,-1          DO K=Nr,1,-1
146           DO j=1,sNy+1           DO j=1,sNy+1
147            DO i=1,sNx+1            DO i=1,sNx+1
# Line 109  C--   Add source term arising from w=d/d Line 164  C--   Add source term arising from w=d/d
164        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
165         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
166  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
167          IF ( nonHydrostatic ) THEN          IF ( nonHydrostatic .AND. zeroPsNH ) THEN
168           DO j=1,sNy           DO j=1,sNy
169            DO i=1,sNx            DO i=1,sNx
170             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
171               IF ( ks.LE.Nr ) THEN
172                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
173         &       -freeSurfFac*_rA(i,j,bi,bj)/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)/deltaTMom/deltaTfreesurf
177         &         * etaH(i,j,bi,bj)
178               ENDIF
179              ENDDO
180             ENDDO
181            ELSEIF ( nonHydrostatic ) THEN
182             DO j=1,sNy
183              DO i=1,sNx
184               ks = ksurfC(i,j,bi,bj)
185               IF ( ks.LE.Nr ) THEN
186                cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
187       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
188       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
189       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
190             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)
191       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf       &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
192       &         *( etaN(i,j,bi,bj)       &         *( etaN(i,j,bi,bj)
193       &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )       &           +phi_nh(i,j,ks,bi,bj)*horiVertRatio/gravity )
194               ENDIF
195            ENDDO            ENDDO
196           ENDDO           ENDDO
197          ELSEIF ( exactConserv ) THEN          ELSEIF ( exactConserv ) THEN
198  #else  #else
199          IF ( exactConserv ) THEN          IF ( exactConserv ) THEN
200  #endif  #endif /* ALLOW_NONHYDROSTATIC */
201           DO j=1,sNy           DO j=1,sNy
202            DO i=1,sNx            DO i=1,sNx
203             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 242  C Western boundary
242            ENDIF            ENDIF
243           ENDDO           ENDDO
244          ENDIF          ENDIF
245  #endif  #endif /* ALLOW_OBCS */
246    C-    end bi,bj loops
247         ENDDO         ENDDO
248        ENDDO        ENDDO
249    
250  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
251        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
252         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
253       &                        myThid)       &                        myThid)
254        ENDIF        ENDIF
# Line 187  C     see CG2D.h for the interface to th Line 260  C     see CG2D.h for the interface to th
260        firstResidual=0.        firstResidual=0.
261        lastResidual=0.        lastResidual=0.
262        numIters=cg2dMaxIters        numIters=cg2dMaxIters
263    c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
264        CALL CG2D(        CALL CG2D(
265       U           cg2d_b,       U           cg2d_b,
266       U           cg2d_x,       U           cg2d_x,
# Line 195  C     see CG2D.h for the interface to th Line 269  C     see CG2D.h for the interface to th
269       U           numIters,       U           numIters,
270       I           myThid )       I           myThid )
271        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
272    c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
273    
274  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
275        IF (debugMode) THEN        IF ( debugLevel .GE. debLevB ) THEN
276         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',         CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
277       &                        myThid)       &                        myThid)
278        ENDIF        ENDIF
279  #endif  #endif
280    
281  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) :
282        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
283       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
284         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
285         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
286         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
287         WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
288         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
289         WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
290         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
291         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
292            _END_MASTER( myThid )
293           ENDIF
294        ENDIF        ENDIF
295    
296  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
# Line 266  C Western boundary Line 343  C Western boundary
343            ENDIF            ENDIF
344            ENDDO            ENDDO
345           ENDIF           ENDIF
346  #endif  #endif /* ALLOW_OBCS */
347    
348             IF ( usingZCoords ) THEN
349    C-       Z coordinate: assume surface @ level k=1
350               tmpFac = freeSurfFac
351             ELSE
352    C-       Other than Z coordinate: no assumption on surface level index
353               tmpFac = 0.
354               DO j=1,sNy
355                DO i=1,sNx
356                  ks = ksurfC(i,j,bi,bj)
357                  IF ( ks.LE.Nr ) THEN
358                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
359         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
360         &                          *_rA(i,j,bi,bj)/deltaTmom
361                  ENDIF
362                ENDDO
363               ENDDO
364             ENDIF
365           K=1           K=1
366             kp1 = MIN(k+1,Nr)
367             maskKp1 = 1.
368             IF (k.GE.Nr) maskKp1 = 0.
369           DO j=1,sNy           DO j=1,sNy
370            DO i=1,sNx            DO i=1,sNx
371             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)
372       &       +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)
373       &       -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)
374       &       +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)
375       &       -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 )
376       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom       &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
377       &          -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
378       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
379            ENDDO            ENDDO
380           ENDDO           ENDDO
381           DO K=2,Nr-1           DO K=2,Nr
382              kp1 = MIN(k+1,Nr)
383              maskKp1 = 1.
384              IF (k.GE.Nr) maskKp1 = 0.
385            DO j=1,sNy            DO j=1,sNy
386             DO i=1,sNx             DO i=1,sNx
387              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)
388       &       +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)
389       &       -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)
390       &       +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)
391       &       -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 )
392       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j,k  ,bi,bj)*maskC(i,j,k-1,bi,bj)
393       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,kp1,bi,bj)*maskKp1
394       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
395    
396             ENDDO             ENDDO
397            ENDDO            ENDDO
398           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  
399    
400  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
401           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 335  C Western boundary Line 422  C Western boundary
422            ENDDO            ENDDO
423            ENDDO            ENDDO
424           ENDIF           ENDIF
425  #endif  #endif /* ALLOW_OBCS */
426    C-    end bi,bj loops
427          ENDDO ! bi          ENDDO
428         ENDDO ! bj         ENDDO
429    
430        firstResidual=0.        firstResidual=0.
431        lastResidual=0.        lastResidual=0.
432        numIters=cg2dMaxIters        numIters=cg3dMaxIters
433          CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
434        CALL CG3D(        CALL CG3D(
435       U           cg3d_b,       U           cg3d_b,
436       U           phi_nh,       U           phi_nh,
# Line 351  C Western boundary Line 439  C Western boundary
439       U           numIters,       U           numIters,
440       I           myThid )       I           myThid )
441        _EXCH_XYZ_R8(phi_nh, myThid )        _EXCH_XYZ_R8(phi_nh, myThid )
442          CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
443    
444        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,        IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
445       &                                    myTime-deltaTClock) ) THEN       &   ) THEN
446         _BEGIN_MASTER( myThid )         IF ( debugLevel .GE. debLevA ) THEN
447         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual          _BEGIN_MASTER( myThid )
448         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
449         WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
450         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
451         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
452         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)          WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
453         _END_MASTER( )          CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
454            _END_MASTER( myThid )
455           ENDIF
456        ENDIF        ENDIF
457    
458    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
459          IF ( zeroPsNH ) THEN
460           DO bj=myByLo(myThid),myByHi(myThid)
461            DO bi=myBxLo(myThid),myBxHi(myThid)
462    
463             IF ( usingZCoords ) THEN
464    C-       Z coordinate: assume surface @ level k=1
465              DO k=2,Nr
466               DO j=1-OLy,sNy+OLy
467                DO i=1-OLx,sNx+OLx
468                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
469         &                           - phi_nh(i,j,1,bi,bj)
470                ENDDO
471               ENDDO
472              ENDDO
473              DO j=1-OLy,sNy+OLy
474               DO i=1-OLx,sNx+OLx
475                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
476         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
477                 phi_nh(i,j,1,bi,bj) = 0.
478               ENDDO
479              ENDDO
480             ELSE
481    C-       Other than Z coordinate: no assumption on surface level index
482              DO j=1-OLy,sNy+OLy
483               DO i=1-OLx,sNx+OLx
484                ks = ksurfC(i,j,bi,bj)
485                IF ( ks.LE.Nr ) THEN
486                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
487         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
488                 DO k=Nr,1,-1
489                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
490         &                            - phi_nh(i,j,ks,bi,bj)
491                 ENDDO
492                ENDIF
493               ENDDO
494              ENDDO
495             ENDIF
496    
497            ENDDO
498           ENDDO
499          ENDIF
500    
501          ENDIF
502    #endif /* ALLOW_NONHYDROSTATIC */
503    
504    #ifdef TIME_PER_TIMESTEP
505    CCE107 Time per timestep information
506          _BEGIN_MASTER( myThid )
507          CALL TIMER_GET_TIME( utnew, stnew, wtnew )
508    C Only output timing information after the 1st timestep
509          IF ( wtold .NE. 0.0D0 ) THEN
510            WRITE(msgBuf,'(A34,3F10.6)')
511         $        'User, system and wallclock time:', utnew - utold,
512         $        stnew - stold, wtnew - wtold
513             CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
514        ENDIF        ENDIF
515          utold = utnew
516          stold = stnew
517          wtold = wtnew
518          _END_MASTER( myThid )
519  #endif  #endif
520    
521        RETURN        RETURN
522        END        END
523    
524    #ifdef TIME_PER_TIMESTEP
525    CCE107 Initialization of common block for per timestep timing
526          BLOCK DATA settimers
527    C     !TIMING VARIABLES
528    C     == Timing variables ==
529          REAL*8 utnew, utold, stnew, stold, wtnew, wtold
530          COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
531          DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
532          END
533    #endif
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22