/[MITgcm]/MITgcm/pkg/kpp/kpp_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/kpp/kpp_calc.F

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

revision 1.3 by heimbach, Wed Sep 13 17:07:11 2000 UTC revision 1.4 by heimbach, Mon Nov 13 16:37:02 2000 UTC
# Line 121  c     myTime - Current time in simulatio Line 121  c     myTime - Current time in simulatio
121    
122  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
123    
124    c Local constants
125    c     minusone, p0, p5, p25, p125, p0625
126    c     imin, imax, jmin, jmax  - array computation indices
127    
128          _RL        minusone
129          parameter( minusone=-1.0)
130          _KPP_RL    p0    , p5    , p25     , p125      , p0625
131          parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
132          integer    imin      , imax        , jmin      , jmax
133    #ifdef FRUGAL_KPP
134          parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )
135    #else
136          parameter( imin=-2   , imax=sNx+3  , jmin=-2   , jmax=sNy+3   )
137    #endif
138    
139  c Local arrays and variables  c Local arrays and variables
140  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
141  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)
# Line 145  c     zRef   (nx,ny)       - Reference d Line 160  c     zRef   (nx,ny)       - Reference d
160  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
161  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
162    
163        _RS worka (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
164        _RS workb (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        integer work1 ( ibot:itop    , jbot:jtop                    )
165  #ifdef FRUGAL_KPP        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )
166        integer work1(sNx,sNy)        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )
167        _RS work2    (sNx,sNy)        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )
168        _RS ustar    (sNx,sNy)        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )
169        _RS bo       (sNx,sNy)        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )
170        _RS bosol    (sNx,sNy)        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )
171        _RS shsq     (sNx,sNy,Nr)        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )
172        _RS dVsq     (sNx,sNy,Nr)        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )
173        _RS dbloc    (sNx,sNy,Nr)        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )
174        _RS Ritop    (sNx,sNy,Nr)        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )
175        _RS vddiff   (sNx,sNy,0:Nrp1,mdiff)        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )
       _RS ghat     (sNx,sNy,Nr)  
       _RS hbl      (sNx,sNy)  
 #ifdef KPP_ESTIMATE_UREF  
       _RS z0       (sNx,sNy)  
       _RS zRef     (sNx,sNy)  
       _RS uRef     (sNx,sNy)  
       _RS vRef     (sNx,sNy)  
 #endif /* KPP_ESTIMATE_UREF */  
 #else /* FRUGAL_KPP */  
       integer work1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS work2    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS ustar    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS bo       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS bosol    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS shsq     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS dVsq     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS dbloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS Ritop    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS vddiff   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:Nrp1,mdiff)  
       _RS ghat     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS hbl      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
176  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
177        _RS z0       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )
178        _RS zRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )
179        _RS uRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )
180        _RS vRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )
181  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
 #endif /* FRUGAL_KPP */  
   
 c     imin,imax,jmin,jmax  - array indices  
       integer    imin   , imax      , jmin   , jmax  
       parameter( imin=-2, imax=sNx+3, jmin=-2, jmax=sNy+3 )  
   
 c     mixing process switches  
       logical    lri  
       parameter( lri = .true. )  
   
       _RS        m1  
       parameter( m1=-1.0)  
       _RS        p0    , p5    , p25     , p125      , p0625  
       parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )  
182                
183        _RL     tempVar        _KPP_RL tempvar1, tempvar2
184        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, im1, ip1, jm1, jp1
185    
186  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
187        _RS     dBdz1, dBdz2, ustarX, ustarY        _KPP_RL dBdz1, dBdz2, ustarX, ustarY
188  #endif  #endif
189    
190  c     Check to see if new vertical mixing coefficient should be computed now?  c     Check to see if new vertical mixing coefficient should be computed now?
# Line 249  c--------------------------------------- Line 229  c---------------------------------------
229        CALL STATEKPP(        CALL STATEKPP(
230       I       bi, bj, myThid       I       bi, bj, myThid
231       O     , work2, dbloc, Ritop       O     , work2, dbloc, Ritop
232  #ifdef FRUGAL_KPP       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)
      O     , vddiff(1    ,1    ,1,1), vddiff(1    ,1    ,1,2)  
 #else  
      O     , vddiff(1-OLx,1-OLy,1,1), vddiff(1-OLx,1-OLy,1,2)  
 #endif  
233       &     )       &     )
234        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
235    
 #ifdef KPP_SMOOTH_DBLOC  
 c     horizontally smooth dbloc with a 121 filter  
 c     (stored in ghat to save space)  
   
236        DO k = 1, Nr        DO k = 1, Nr
237           CALL SMOOTH_HORIZ_RS (           DO j = jbot, jtop
238       I        k, bi, bj,              DO i = ibot, itop
      I        dbloc(1-OLx,1-OLy,k),  
      O        ghat (1-OLx,1-OLy,k) )  
       ENDDO  
   
 #else /* KPP_SMOOTH_DBLOC */  
   
       DO k = 1, Nr  
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             DO i = 1, sNx  
 #else  
          DO j = 1-OLy, sNy+OLy  
             DO i = imin, imax  
 #endif  
239                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
240              ENDDO              ENDDO
241           ENDDO           ENDDO
242        ENDDO        ENDDO
243    
244    #ifdef KPP_SMOOTH_DBLOC
245    c     horizontally smooth dbloc with a 121 filter
246    c     smooth dbloc stored in ghat to save space
247    c     dbloc(k) is buoyancy gradientnote between k and k+1
248    c     levels therefore k+1 mask must be used
249    
250          DO k = 1, Nr-1
251             CALL KPP_SMOOTH_HORIZ (
252         I        k+1, bi, bj,
253         U        ghat (ibot,jbot,k) )
254          ENDDO
255    
256  #endif /* KPP_SMOOTH_DBLOC */  #endif /* KPP_SMOOTH_DBLOC */
257    
258  #ifdef KPP_SMOOTH_DENS  #ifdef KPP_SMOOTH_DENS
259  c     horizontally smooth density related quantities with 121 filters  c     horizontally smooth density related quantities with 121 filters
260        CALL SMOOTH_HORIZ_RS (        CALL KPP_SMOOTH_HORIZ (
261       I     k, bi, bj,       I     1, bi, bj,
262       I     work2,       U     work2 )
      O     work2 )  
263        DO k = 1, Nr        DO k = 1, Nr
264           CALL SMOOTH_HORIZ_RS (           CALL KPP_SMOOTH_HORIZ (
265       I        k, bi, bj,       I        k+1, bi, bj,
266       I        dbloc (1-OLx,1-OLy,k)  ,       U        dbloc (ibot,jbot,k) )
267       O        dbloc (1-OLx,1-OLy,k)   )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
268       I        k, bi, bj,       I        k, bi, bj,
269       I        Ritop (1-OLx,1-OLy,k)  ,       U        Ritop (ibot,jbot,k)  )
270       O        Ritop (1-OLx,1-OLy,k)   )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
271       I        k, bi, bj,       I        k, bi, bj,
272       I        vddiff(1-OLx,1-OLy,k,1),       U        vddiff(ibot,jbot,k,1) )
273       O        vddiff(1-OLx,1-OLy,k,1) )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
274       I        k, bi, bj,       I        k, bi, bj,
275       I        vddiff(1-OLx,1-OLy,k,2),       U        vddiff(ibot,jbot,k,2) )
      O        vddiff(1-OLx,1-OLy,k,2) )  
276        ENDDO        ENDDO
277  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
278    
279        DO k = 1, Nr        DO k = 1, Nr
280  #ifdef FRUGAL_KPP           DO j = jbot, jtop
281           DO j = 1, sNy              DO i = ibot, itop
             DO i = 1, sNx  
 #else  
          DO j = 1-OLy, sNy+OLy  
             DO i = 1-OLx, sNx+OLx  
 #endif  
282    
283  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
284  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
# Line 351  c                     beta *SurfaceTende Line 311  c                     beta *SurfaceTende
311  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)
312  c------------------------------------------------------------------------  c------------------------------------------------------------------------
313    
314  #ifdef FRUGAL_KPP  c initialize arrays to zero
315        DO j = 1, sNy        DO j = jbot, jtop
316           jp1 = j + 1           DO i = ibot, itop
317           DO i = 1, sNx              ustar(i,j) = p0
318  #else              bo   (I,J) = p0
319                bosol(I,J) = p0
320             END DO
321          END DO
322    
323        DO j = jmin, jmax        DO j = jmin, jmax
324         jp1 = j + 1         jp1 = j + 1
325         DO i = imin, imax         DO i = imin, imax
 #endif  
326          ip1 = i+1          ip1 = i+1
327          tempVar =          tempVar1 =
328       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *
329       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +
330       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *
331       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))
332          if ( tempVar .lt. (epsln*epsln) ) then          if ( tempVar1 .lt. (phepsi*phepsi) ) then
333             ustar(i,j) = SQRT( epsln * p5 * delZ(1) )             ustar(i,j) = SQRT( phepsi * p5 * delZ(1) )
334          else          else
335             ustar(i,j) = SQRT( SQRT( tempVar ) * p5 * delZ(1) )             tempVar2 =  SQRT( tempVar1 ) * p5 * delZ(1)
336               ustar(i,j) = SQRT( tempVar2 )
337          endif          endif
338          bo(I,J) = - gravity *          bo(I,J) = - gravity *
339       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +
# Line 381  c--------------------------------------- Line 345  c---------------------------------------
345         END DO         END DO
346        END DO        END DO
347    
 #ifndef FRUGAL_KPP  
 c     set array edges to zero  
       DO j = jmin, jmax  
          DO i = 1-OLx, imin-1  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
          DO i = imax+1, sNx+OLx  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
       DO i = 1-OLx, sNx+OLx  
          DO j = 1-OLy, jmin-1  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
          DO j = jmax+1, sNy+OLy  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
 #endif  
   
348  c------------------------------------------------------------------------  c------------------------------------------------------------------------
349  c     velocity shear  c     velocity shear
350  c     --------------  c     --------------
# Line 418  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k)) Line 354  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))
354  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces
355  c------------------------------------------------------------------------  c------------------------------------------------------------------------
356    
357    c initialize arrays to zero
358          DO k = 1, Nr
359             DO j = jbot, jtop
360                DO i = ibot, itop
361                   shsq(i,j,k) = p0
362                   dVsq(i,j,k) = p0
363                END DO
364             END DO
365          END DO
366    
367  c     dVsq computation  c     dVsq computation
368    
369  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
# Line 428  c     thickness in the model.  First det Line 374  c     thickness in the model.  First det
374  c     Second zRef = espilon * hMix.  Third determine roughness length  c     Second zRef = espilon * hMix.  Third determine roughness length
375  c     scale z0.  Third estimate reference velocity.  c     scale z0.  Third estimate reference velocity.
376    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          jp1 = j + 1  
          DO i = 1, sNx  
 #else  
377        DO j = jmin, jmax        DO j = jmin, jmax
378           jp1 = j + 1           jp1 = j + 1
379           DO i = imin, imax           DO i = imin, imax
 #endif /* FRUGAL_KPP */  
380              ip1 = i + 1              ip1 = i + 1
381    
382  c     Determine mixed layer depth hMix as the shallowest depth at which  c     Determine mixed layer depth hMix as the shallowest depth at which
# Line 465  c     Linearly interpolate to find hMix. Line 405  c     Linearly interpolate to find hMix.
405              ENDIF              ENDIF
406    
407  c     Compute roughness length scale z0 subject to 0 < z0  c     Compute roughness length scale z0 subject to 0 < z0
408                 tempVar = p5 * (                 tempVar1 = p5 * (
409       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
410       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
411       &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *       &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *
# Line 474  c     Compute roughness length scale z0 Line 414  c     Compute roughness length scale z0
414       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
415       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
416       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
417                 if ( tempVar .lt. (epsln*epsln) ) then                 if ( tempVar1 .lt. (epsln*epsln) ) then
418                    tempVar = epsln                    tempVar2 = epsln
419                 else                 else
420                    tempVar = SQRT ( tempVar )                    tempVar2 = SQRT ( tempVar1 )
421                 endif                 endif
422                 z0(i,j) = rF(2) *                 z0(i,j) = rF(2) *
423       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /
424       &                     ( rF(3) - rF(2) ) -       &                     ( rF(3) - rF(2) ) -
425       &                     tempVar * vonK /       &                     tempVar2 * vonK /
426       &                     MAX ( ustar(i,j), phepsi ) )       &                     MAX ( ustar(i,j), phepsi ) )
427                 z0(i,j) = MAX ( z0(i,j), phepsi )                 z0(i,j) = MAX ( z0(i,j), phepsi )
428    
# Line 500  c     Estimate reference velocity uRef a Line 440  c     Estimate reference velocity uRef a
440       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5
441                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +
442       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5
443                    tempVar = ustarX * ustarX + ustarY * ustarY                    tempVar1 = ustarX * ustarX + ustarY * ustarY
444                    if ( tempVar .lt. (epsln*epsln) ) then                    if ( tempVar1 .lt. (epsln*epsln) ) then
445                       tempVar = epsln                       tempVar2 = epsln
446                    else                    else
447                       tempVar = SQRT ( tempVar )                       tempVar2 = SQRT ( tempVar1 )
448                    endif                    endif
449                    tempVar = ustar(i,j) *                    tempVar2 = ustar(i,j) *
450       &                 ( LOG ( zRef(i,j) / rF(2) ) +       &                 ( LOG ( zRef(i,j) / rF(2) ) +
451       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /
452       &                 vonK / tempVar       &                 vonK / tempVar2
453                    uRef(i,j) = uRef(i,j) + ustarX * tempVar                    uRef(i,j) = uRef(i,j) + ustarX * tempVar2
454                    vRef(i,j) = vRef(i,j) + ustarY * tempVar                    vRef(i,j) = vRef(i,j) + ustarY * tempVar2
455                 ENDIF                 ENDIF
456    
457           END DO           END DO
458        END DO        END DO
459    
       IF (KPPmixingMaps) THEN  
 #ifdef FRUGAL_KPP  
          CALL PRINT_MAPRS(  
      I        zRef, 'zRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        z0, 'z0', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        uRef, 'uRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        vRef, 'vRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #else  
          CALL PRINT_MAPRS(  
      I        zRef, 'zRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        z0, 'z0', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        uRef, 'uRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        vRef, 'vRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #endif  
       ENDIF  
   
460        DO k = 1, Nr        DO k = 1, Nr
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
461           DO j = jmin, jmax           DO j = jmin, jmax
462              jm1 = j - 1              jm1 = j - 1
463              jp1 = j + 1              jp1 = j + 1
464              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
465                 im1 = i - 1                 im1 = i - 1
466                 ip1 = i + 1                 ip1 = i + 1
467                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 604  c     Estimate reference velocity uRef a Line 499  c     Estimate reference velocity uRef a
499  #else /* KPP_ESTIMATE_UREF */  #else /* KPP_ESTIMATE_UREF */
500    
501        DO k = 1, Nr        DO k = 1, Nr
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
502           DO j = jmin, jmax           DO j = jmin, jmax
503              jm1 = j - 1              jm1 = j - 1
504              jp1 = j + 1              jp1 = j + 1
505              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
506                 im1 = i - 1                 im1 = i - 1
507                 ip1 = i + 1                 ip1 = i + 1
508                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 654  c     Estimate reference velocity uRef a Line 542  c     Estimate reference velocity uRef a
542  c     shsq computation  c     shsq computation
543        DO k = 1, Nrm1        DO k = 1, Nrm1
544           kp1 = k + 1           kp1 = k + 1
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
545           DO j = jmin, jmax           DO j = jmin, jmax
546              jm1 = j - 1              jm1 = j - 1
547              jp1 = j + 1              jp1 = j + 1
548              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
549                 im1 = i - 1                 im1 = i - 1
550                 ip1 = i + 1                 ip1 = i + 1
551                 shsq(i,j,k) = p5 * (                 shsq(i,j,k) = p5 * (
# Line 699  c     shsq computation Line 580  c     shsq computation
580           END DO           END DO
581        END DO        END DO
582    
 c     shsq @ Nr computation  
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
       DO j = jmin, jmax  
          DO i = imin, imax  
 #endif  
             shsq(i,j,Nr) = p0  
          END DO  
       END DO  
   
 #ifndef FRUGAL_KPP  
 c     set array edges to zero  
       DO k = 1, Nr  
          DO j = jmin, jmax  
             DO i = 1-OLx, imin-1  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
             DO i = imax+1, sNx+OLx  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
          END DO  
          DO i = 1-OLx, sNx+OLx  
             DO j = 1-OLy, jmin-1  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
             DO j = jmax+1, sNy+OLy  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
          END DO  
       END DO  
 #endif  
   
583  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
584  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
585  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
586    
587  #ifdef FRUGAL_KPP        DO j = jbot, jtop
588        DO j = 1, sNy           DO i = ibot, itop
          DO i = 1, sNx  
 #else  
       DO j = 1-OLy, sNy+OLy  
          DO i = 1-OLx, sNx+OLx  
 #endif  
589              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
590              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
591           END DO           END DO
592        END DO        END DO
593        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)
594        CALL KPPMIX (        CALL KPPMIX (
595       I       lri, work1, shsq, dVsq, ustar       I       mytime, mythid
596         I     , work1, shsq, dVsq, ustar
597       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
598       I     , ikey       I     , ikey
599       O     , vddiff       O     , vddiff
600       U     , ghat       U     , ghat
601       O     , hbl       O     , hbl )
      &        )  
602    
603        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
604    
       IF (KPPmixingMaps) THEN  
 #ifdef FRUGAL_KPP  
        CALL PRINT_MAPRS(  
      I      hbl, 'hbl', PRINT_MAP_XY,  
      I      1, sNx, 1, sNy, 1, 1, 1, 1,  
      I      1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #else  
        CALL PRINT_MAPRS(  
      I      hbl, 'hbl', PRINT_MAP_XY,  
      I      1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I      1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #endif  
        ENDIF  
   
605  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
606  CADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey  cph( storing not necessary
607    cphCADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey
608    cph)
609  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
610    
611  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
612  c     zero out land values,  c     zero out land values and transfer to global variables
 c     make sure coefficients are within reasonable bounds,  
 c     and transfer to global variables  
613  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
614    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
615        DO j = jmin, jmax        DO j = jmin, jmax
616           DO i = imin, imax         DO i = imin, imax
617  #endif          DO k = 1, Nr
618              DO k = 1, Nr           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)
619  c KPPviscAz           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)
620                 tempVar = min( maxKPPviscAz(k), vddiff(i,j,k-1,1) )           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)
621                 tempVar = max( minKPPviscAz, tempVar )           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)
622                 KPPviscAz(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)          END DO
623  c KPPdiffKzS          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)
624                 tempVar = min( maxKPPdiffKzS, vddiff(i,j,k-1,2) )         END DO
                tempVar = max( minKPPdiffKzS, tempVar )  
                KPPdiffKzS(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)  
 c KPPdiffKzT  
                tempVar = min( maxKPPdiffKzT, vddiff(i,j,k-1,3) )  
                tempVar = max( minKPPdiffKzT, tempVar )  
                KPPdiffKzT(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)  
 c KPPghat  
                tempVar = min( maxKPPghat, ghat(i,j,k) )  
                tempVar = max( minKPPghat, tempVar )  
                KPPghat(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)  
             END DO  
 c KPPhbl: set to -zgrid(1) over land  
             KPPhbl(i,j,bi,bj) = (hbl(i,j) + zgrid(1))  
      &           * pMask(i,j,1,bi,bj) -  
      &           zgrid(1)  
          END DO  
625        END DO        END DO
626  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
627        _EXCH_XYZ_R8(KPPviscAz  , myThid )        _EXCH_XYZ_R8(KPPviscAz  , myThid )
# Line 829  c KPPhbl: set to -zgrid(1) over land Line 633  c KPPhbl: set to -zgrid(1) over land
633    
634  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
635  c     horizontal smoothing of vertical viscosity  c     horizontal smoothing of vertical viscosity
 c     as coded requires FRUGAL_KPP and OLx=4, OLy=4  
 c     alternatively could recode with OLx=5, OLy=5  
   
636        DO k = 1, Nr        DO k = 1, Nr
637           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
638       I        k, bi, bj,       I        k, bi, bj,
639       I        KPPviscAz(1-OLx,1-OLy,k,bi,bj),       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj) )
      O        KPPviscAz(1-OLx,1-OLy,k,bi,bj) )  
640        END DO        END DO
641          _EXCH_XYZ_R8(KPPviscAz  , myThid )
642  #endif /* KPP_SMOOTH_VISC */  #endif /* KPP_SMOOTH_VISC */
643    
644  #ifdef KPP_SMOOTH_DIFF  #ifdef KPP_SMOOTH_DIFF
645  c     horizontal smoothing of vertical diffusivity  c     horizontal smoothing of vertical diffusivity
 c     as coded requires FRUGAL_KPP and OLx=4, OLy=4  
 c     alternatively could recode with OLx=5, OLy=5  
   
646        DO k = 1, Nr        DO k = 1, Nr
647           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
648       I        k, bi, bj,       I        k, bi, bj,
649       I        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )
650       O        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )           CALL SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RL (  
651       I        k, bi, bj,       I        k, bi, bj,
652       I        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) )
      O        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) )  
653        END DO        END DO
654          _EXCH_XYZ_R8(KPPdiffKzS , myThid )
655          _EXCH_XYZ_R8(KPPdiffKzT , myThid )
656  #endif /* KPP_SMOOTH_DIFF */  #endif /* KPP_SMOOTH_DIFF */
657    
   
658  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
659  C     the bottom of the mixing layer.  C     the bottom of the mixing layer.
660        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 866  C     the bottom of the mixing layer. Line 663  C     the bottom of the mixing layer.
663           ENDDO           ENDDO
664        ENDDO        ENDDO
665        CALL SWFRAC(        CALL SWFRAC(
666       I     (sNx+2*OLx)*(sNy+2*OLy), m1, worka,       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
667       O     workb )       I     mytime, mythid,
668         U     worka )
669        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
670           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
671              KPPfrac(i,j,bi,bj) = workb(i,j)              KPPfrac(i,j,bi,bj) = worka(i,j)
672           ENDDO           ENDDO
673        ENDDO        ENDDO
674    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22