/[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.1 by adcroft, Wed Jun 21 19:45:47 2000 UTC revision 1.11 by heimbach, Sat Jul 13 03:12:30 2002 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
5    
# Line 83  c     KPPfrac     - Fraction of short-wa Line 84  c     KPPfrac     - Fraction of short-wa
84    
85  c--   KPP_CALC computes vertical viscosity and diffusivity for region  c--   KPP_CALC computes vertical viscosity and diffusivity for region
86  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires
87  c     values of uVel, vVel, fu, fv in the region (-2:sNx+4,-2:sNy+4).  c     values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the
88    c     region (-2:sNx+4,-2:sNy+4).
89  c     Hence overlap region needs to be set OLx=4, OLy=4.  c     Hence overlap region needs to be set OLx=4, OLy=4.
90  c     When option FRUGAL_KPP is used, computation in overlap regions  c     When option FRUGAL_KPP is used, computation in overlap regions
91  c     is replaced with exchange calls hence reducing overlap requirements  c     is replaced with exchange calls hence reducing overlap requirements
# Line 118  c     myTime - Current time in simulatio Line 120  c     myTime - Current time in simulatio
120    
121  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
122    
123    c Local constants
124    c     minusone, p0, p5, p25, p125, p0625
125    c     imin, imax, jmin, jmax  - array computation indices
126    
127          _RL        minusone
128          parameter( minusone=-1.0)
129          _KPP_RL    p0    , p5    , p25     , p125      , p0625
130          parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
131          integer    imin      , imax        , jmin      , jmax
132    #ifdef FRUGAL_KPP
133          parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )
134    #else
135          parameter( imin=-2   , imax=sNx+3  , jmin=-2   , jmax=sNy+3   )
136    #endif
137    
138  c Local arrays and variables  c Local arrays and variables
139  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
140  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)
# Line 142  c     zRef   (nx,ny)       - Reference d Line 159  c     zRef   (nx,ny)       - Reference d
159  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
160  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
161    
162        _RS worka (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
163        _RS workb (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        integer work1 ( ibot:itop    , jbot:jtop                    )
164  #ifdef FRUGAL_KPP        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )
165        integer work1(sNx,sNy)        _KPP_RL work3 ( ibot:itop    , jbot:jtop                    )
166        _RS work2    (sNx,sNy)        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )
167        _RS ustar    (sNx,sNy)        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )
168        _RS bo       (sNx,sNy)        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )
169        _RS bosol    (sNx,sNy)        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )
170        _RS shsq     (sNx,sNy,Nr)        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )
171        _RS dVsq     (sNx,sNy,Nr)        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )
172        _RS dbloc    (sNx,sNy,Nr)        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )
173        _RS Ritop    (sNx,sNy,Nr)        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )
174        _RS vddiff   (sNx,sNy,0:Nrp1,mdiff)        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )
175        _RS ghat     (sNx,sNy,Nr)        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )
       _RS hbl      (sNx,sNy)  
176  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
177        _RS z0       (sNx,sNy)        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )
178        _RS zRef     (sNx,sNy)        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )
179        _RS uRef     (sNx,sNy)        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )
180        _RS vRef     (sNx,sNy)        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )
181  #endif /* KPP_ESTIMATE_UREF */  #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)  
 #ifdef KPP_ESTIMATE_UREF  
       _RS z0       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS zRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS uRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS vRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 #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        p0    , p5    , p25     , p125      , p0625  
       parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )  
182                
183        _RS     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    
       IF (use_KPPmixing) THEN  
   
 CADJ STORE fu     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE fv     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
   
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?
191        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.
192       1     myTime .EQ. startTime ) THEN       1     myTime .EQ. startTime ) THEN
# 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+1, bi, bj,
266         U        dbloc (ibot,jbot,k) )
267             CALL KPP_SMOOTH_HORIZ (
268       I        k, bi, bj,       I        k, bi, bj,
269       I        dbloc (1-OLx,1-OLy,k)  ,       U        Ritop (ibot,jbot,k)  )
270       O        dbloc (1-OLx,1-OLy,k)   )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
271       I        k, bi, bj,       I        k, bi, bj,
272       I        Ritop (1-OLx,1-OLy,k)  ,       U        vddiff(ibot,jbot,k,1) )
273       O        Ritop (1-OLx,1-OLy,k)   )           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,1),       U        vddiff(ibot,jbot,k,2) )
      O        vddiff(1-OLx,1-OLy,k,1) )  
          CALL SMOOTH_HORIZ_RS (  
      I        k, bi, bj,  
      I        vddiff(1-OLx,1-OLy,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 340  c     so that the subroutine "bldepth" w Line 300  c     so that the subroutine "bldepth" w
300           END DO           END DO
301        END DO        END DO
302    
303    cph(
304    cph  this avoids a single or double recomp./call of statekpp
305    CADJ store work2              = comlev1_kpp, key = ikey
306    #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE
307    CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikey
308    CADJ store vddiff             = comlev1_kpp, key = ikey
309    #endif
310    cph)
311    
312  c------------------------------------------------------------------------  c------------------------------------------------------------------------
313  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
314  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
315  c     taux / rho = fu * delZ(1)                                   (N/m^2)  c     taux / rho = SurfaceTendencyU * drF(1)                     (N/m^2)
316  c     tauy / rho = fv * delZ(1)                                   (N/m^2)  c     tauy / rho = SurfaceTendencyV * drF(1)                     (N/m^2)
317  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
318  c     bo    = - g * (alpha*Qnet + beta*EmPmR) * delZ(1) / rho   (m^2/s^3)  c     bo    = - g * ( alpha*SurfaceTendencyT +
319  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c                     beta *SurfaceTendencyS ) * drF(1) / rho  (m^2/s^3)
320    c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
321  c------------------------------------------------------------------------  c------------------------------------------------------------------------
322    
323  #ifdef FRUGAL_KPP  c initialize arrays to zero
324        DO j = 1, sNy        DO j = jbot, jtop
325           jp1 = j + 1           DO i = ibot, itop
          DO i = 1, sNx  
 #else  
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
 #endif  
             ip1 = i+1  
             ustar(i,j) = SQRT( SQRT(  
      &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) *  
      &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) +  
      &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) *  
      &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) ))  
             bo(I,J) = - gravity *  
      &           ( vddiff(I,J,1,1) * Qnet(i,j,bi,bj) +  
      &             vddiff(I,J,1,2) * EmPmR(i,j,bi,bj)  
      &           ) *  
      &           delZ(1) / work2(I,J)  
             bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *  
      &           delZ(1) / work2(I,J)  
          END DO  
       END DO  
   
 #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  
326              ustar(i,j) = p0              ustar(i,j) = p0
327              bo   (I,J) = p0              bo   (I,J) = p0
328              bosol(I,J) = p0              bosol(I,J) = p0
329           END DO           END DO
330        END DO        END DO
331  #endif  
332          DO j = jmin, jmax
333           jp1 = j + 1
334           DO i = imin, imax
335            ip1 = i+1
336            work3(i,j) =
337         &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *
338         &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +
339         &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *
340         &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))
341           END DO
342          END DO
343    cph(
344    CADJ store work3 = comlev1_kpp, key = ikey
345    cph)
346          DO j = jmin, jmax
347           jp1 = j + 1
348           DO i = imin, imax
349            ip1 = i+1
350            if ( work3(i,j) .lt. (phepsi*phepsi) ) then
351               ustar(i,j) = SQRT( phepsi * p5 * drF(1) )
352            else
353               tempVar2 =  SQRT( work3(i,j) ) * p5 * drF(1)
354               ustar(i,j) = SQRT( tempVar2 )
355            endif
356            bo(I,J) = - gravity *
357         &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +
358         &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)
359         &       ) *
360         &       drF(1) / work2(I,J)
361            bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *
362         &       recip_Cp*recip_rhoNil*recip_dRf(1) *
363         &       drF(1) / work2(I,J)
364           END DO
365          END DO
366    
367    cph(
368    CADJ store ustar = comlev1_kpp, key = ikey
369    cph)
370    
371  c------------------------------------------------------------------------  c------------------------------------------------------------------------
372  c     velocity shear  c     velocity shear
# Line 412  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k)) Line 377  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))
377  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
378  c------------------------------------------------------------------------  c------------------------------------------------------------------------
379    
380    c initialize arrays to zero
381          DO k = 1, Nr
382             DO j = jbot, jtop
383                DO i = ibot, itop
384                   shsq(i,j,k) = p0
385                   dVsq(i,j,k) = p0
386                END DO
387             END DO
388          END DO
389    
390  c     dVsq computation  c     dVsq computation
391    
392  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
# Line 422  c     thickness in the model.  First det Line 397  c     thickness in the model.  First det
397  c     Second zRef = espilon * hMix.  Third determine roughness length  c     Second zRef = espilon * hMix.  Third determine roughness length
398  c     scale z0.  Third estimate reference velocity.  c     scale z0.  Third estimate reference velocity.
399    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          jp1 = j + 1  
          DO i = 1, sNx  
 #else  
400        DO j = jmin, jmax        DO j = jmin, jmax
401           jp1 = j + 1           jp1 = j + 1
402           DO i = imin, imax           DO i = imin, imax
 #endif /* FRUGAL_KPP */  
403              ip1 = i + 1              ip1 = i + 1
404    
405  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 459  c     Linearly interpolate to find hMix. Line 428  c     Linearly interpolate to find hMix.
428              ENDIF              ENDIF
429    
430  c     Compute roughness length scale z0 subject to 0 < z0  c     Compute roughness length scale z0 subject to 0 < z0
431                 tempVar = SQRT ( p5 * (                 tempVar1 = p5 * (
432       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
433       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
434       &              (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 467  c     Compute roughness length scale z0 Line 436  c     Compute roughness length scale z0
436       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *
437       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
438       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
439       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) ) )       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
440                   if ( tempVar1 .lt. (epsln*epsln) ) then
441                      tempVar2 = epsln
442                   else
443                      tempVar2 = SQRT ( tempVar1 )
444                   endif
445                 z0(i,j) = rF(2) *                 z0(i,j) = rF(2) *
446       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /
447       &                     ( rF(3) - rF(2) ) -       &                     ( rF(3) - rF(2) ) -
448       &                     tempVar * vonK /       &                     tempVar2 * vonK /
449       &                     MAX ( ustar(i,j), phepsi ) )       &                     MAX ( ustar(i,j), phepsi ) )
450                 z0(i,j) = MAX ( z0(i,j), phepsi )                 z0(i,j) = MAX ( z0(i,j), phepsi )
451    
# Line 485  c     Estimate reference velocity uRef a Line 459  c     Estimate reference velocity uRef a
459                 vRef(i,j) = p5 *                 vRef(i,j) = p5 *
460       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
461                 IF ( zRef(i,j) .LT. drF(1) ) THEN                 IF ( zRef(i,j) .LT. drF(1) ) THEN
462                    ustarX = ( fu(i,j,bi,bj) + fu(ip1,j,bi,bj) ) * p5                    ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +
463                    ustarY = ( fv(i,j,bi,bj) + fu(i,jp1,bi,bj) ) * p5       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5
464                    tempVar = MAX ( phepsi,                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +
465       &                 SQRT ( ustarX * ustarX + ustarY * ustarY ) )       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5
466                    tempVar = ustar(i,j) *                    tempVar1 = ustarX * ustarX + ustarY * ustarY
467                      if ( tempVar1 .lt. (epsln*epsln) ) then
468                         tempVar2 = epsln
469                      else
470                         tempVar2 = SQRT ( tempVar1 )
471                      endif
472                      tempVar2 = ustar(i,j) *
473       &                 ( LOG ( zRef(i,j) / rF(2) ) +       &                 ( LOG ( zRef(i,j) / rF(2) ) +
474       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /
475       &                 vonK / tempVar       &                 vonK / tempVar2
476                    uRef(i,j) = uRef(i,j) + ustarX * tempVar                    uRef(i,j) = uRef(i,j) + ustarX * tempVar2
477                    vRef(i,j) = vRef(i,j) + ustarY * tempVar                    vRef(i,j) = vRef(i,j) + ustarY * tempVar2
478                 ENDIF                 ENDIF
479    
480           END DO           END DO
481        END DO        END DO
482    
       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  
   
483        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  
484           DO j = jmin, jmax           DO j = jmin, jmax
485              jm1 = j - 1              jm1 = j - 1
486              jp1 = j + 1              jp1 = j + 1
487              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
488                 im1 = i - 1                 im1 = i - 1
489                 ip1 = i + 1                 ip1 = i + 1
490                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 587  c     Estimate reference velocity uRef a Line 522  c     Estimate reference velocity uRef a
522  #else /* KPP_ESTIMATE_UREF */  #else /* KPP_ESTIMATE_UREF */
523    
524        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  
525           DO j = jmin, jmax           DO j = jmin, jmax
526              jm1 = j - 1              jm1 = j - 1
527              jp1 = j + 1              jp1 = j + 1
528              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
529                 im1 = i - 1                 im1 = i - 1
530                 ip1 = i + 1                 ip1 = i + 1
531                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 637  c     Estimate reference velocity uRef a Line 565  c     Estimate reference velocity uRef a
565  c     shsq computation  c     shsq computation
566        DO k = 1, Nrm1        DO k = 1, Nrm1
567           kp1 = k + 1           kp1 = k + 1
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
568           DO j = jmin, jmax           DO j = jmin, jmax
569              jm1 = j - 1              jm1 = j - 1
570              jp1 = j + 1              jp1 = j + 1
571              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
572                 im1 = i - 1                 im1 = i - 1
573                 ip1 = i + 1                 ip1 = i + 1
574                 shsq(i,j,k) = p5 * (                 shsq(i,j,k) = p5 * (
# Line 682  c     shsq computation Line 603  c     shsq computation
603           END DO           END DO
604        END DO        END DO
605    
606  c     shsq @ Nr computation  cph(
607  #ifdef FRUGAL_KPP  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE
608        DO j = 1, sNy  CADJ store dvsq, shsq = comlev1_kpp, key = ikey
          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  
609  #endif  #endif
610    cph)
611    
612  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
613  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
614  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
615    
616  #ifdef FRUGAL_KPP        DO j = jbot, jtop
617        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  
618              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
619              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
620           END DO           END DO
621        END DO        END DO
622        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)
623        CALL KPPMIX (        CALL KPPMIX (
624       I       lri, work1, shsq, dVsq, ustar       I       mytime, mythid
625         I     , work1, shsq, dVsq, ustar
626       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
627       I     , ikey       I     , ikey
628       O     , vddiff       O     , vddiff
629       U     , ghat       U     , ghat
630       O     , hbl       O     , hbl )
      &        )  
631    
632        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
633    
       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  
   
 CADJ STORE vddiff, ghat  = uvtape, key = ikey  
   
634  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
635  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  
636  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
637    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
638        DO j = jmin, jmax        DO j = jmin, jmax
639           DO i = imin, imax         DO i = imin, imax
640  #endif          DO k = 1, Nr
641              DO k = 1, Nr           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)
642  c KPPviscAz           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)
643                 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)
644                 tempVar = max( minKPPviscAz, tempVar )           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)
645                 KPPviscAz(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)          END DO
646  c KPPdiffKzS          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)
647                 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  
648        END DO        END DO
649  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
650        _EXCH_XYZ_R8(KPPviscAz  , myThid )        _EXCH_XYZ_R8(KPPviscAz  , myThid )
# Line 810  c KPPhbl: set to -zgrid(1) over land Line 656  c KPPhbl: set to -zgrid(1) over land
656    
657  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
658  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  
   
659        DO k = 1, Nr        DO k = 1, Nr
660           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
661       I        k, bi, bj,       I        k, bi, bj,
662       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) )  
663        END DO        END DO
664          _EXCH_XYZ_R8(KPPviscAz  , myThid )
665  #endif /* KPP_SMOOTH_VISC */  #endif /* KPP_SMOOTH_VISC */
666    
667  #ifdef KPP_SMOOTH_DIFF  #ifdef KPP_SMOOTH_DIFF
668  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  
   
669        DO k = 1, Nr        DO k = 1, Nr
670           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
671       I        k, bi, bj,       I        k, bi, bj,
672       I        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )
673       O        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )           CALL SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RL (  
674       I        k, bi, bj,       I        k, bi, bj,
675       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) )  
676        END DO        END DO
677          _EXCH_XYZ_R8(KPPdiffKzS , myThid )
678          _EXCH_XYZ_R8(KPPdiffKzT , myThid )
679  #endif /* KPP_SMOOTH_DIFF */  #endif /* KPP_SMOOTH_DIFF */
680    
681    cph(
682    cph  crucial: this avoids full recomp./call of kppmix
683    CADJ store KPPhbl = comlev1_kpp, key = ikey
684    cph)
685    
686  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
687  C     the bottom of the mixing layer.  C     the bottom of the mixing layer.
# Line 847  C     the bottom of the mixing layer. Line 691  C     the bottom of the mixing layer.
691           ENDDO           ENDDO
692        ENDDO        ENDDO
693        CALL SWFRAC(        CALL SWFRAC(
694       I     (sNx+2*OLx)*(sNy+2*OLy), -1., worka,       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
695       O     workb )       I     mytime, mythid,
696         U     worka )
697        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
698           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
699              KPPfrac(i,j,bi,bj) = workb(i,j)              KPPfrac(i,j,bi,bj) = worka(i,j)
700           ENDDO           ENDDO
701        ENDDO        ENDDO
           
 CADJ STORE KPPhbl    (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPviscAz (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzT(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzS(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
702    
703        ENDIF        ENDIF
       ENDIF  
704    
705  #endif ALLOW_KPP  #endif /* ALLOW_KPP */
706    
707          RETURN
708          END
709    
710          subroutine KPP_CALC_DUMMY(
711         I     bi, bj, myTime, myThid )
712    C     /==========================================================\
713    C     | SUBROUTINE KPP_CALC_DUMMY                                |
714    C     | o Compute all KPP fields defined in KPP.h                |
715    C     | o Dummy routine for TAMC
716    C     |==========================================================|
717    C     | This subroutine serves as an interface between MITGCMUV  |
718    C     | code and NCOM 1-D routines in kpp_routines.F             |
719    C     \==========================================================/
720          IMPLICIT NONE
721    
722    #include "SIZE.h"
723    #include "EEPARAMS.h"
724    #include "PARAMS.h"
725    #include "KPP.h"
726    #include "KPP_PARAMS.h"
727    #include "GRID.h"
728    
729    c Routine arguments
730    c     bi, bj - array indices on which to apply calculations
731    c     myTime - Current time in simulation
732    
733          INTEGER bi, bj
734          INTEGER myThid
735          _RL     myTime
736    
737    #ifdef ALLOW_KPP
738    
739    c Local constants
740          integer i, j, k
741    
742          DO j=1-OLy,sNy+OLy
743             DO i=1-OLx,sNx+OLx
744                KPPhbl (i,j,bi,bj) = 1.0
745                KPPfrac(i,j,bi,bj) = 0.0
746                DO k = 1,Nr
747                   KPPghat   (i,j,k,bi,bj) = 0.0
748                   KPPviscAz (i,j,k,bi,bj) = viscAz
749                   KPPdiffKzT(i,j,k,bi,bj) = diffKzT
750                   KPPdiffKzS(i,j,k,bi,bj) = diffKzS
751                ENDDO
752             ENDDO
753          ENDDO
754          
755    #endif
756        RETURN        RETURN
757        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22