/[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.15 by heimbach, Fri Mar 7 04:19:01 2003 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 99  c     to OLx=1, OLy=1. Line 100  c     to OLx=1, OLy=1.
100  #include "FFIELDS.h"  #include "FFIELDS.h"
101  #include "GRID.h"  #include "GRID.h"
102    
103    #ifdef ALLOW_SEAICE
104    #include "SEAICE_EXTERNAL.h"
105    #endif ALLOW_SEAICE
106    
107  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
108  #include "tamc.h"  #include "tamc.h"
109  #include "tamc_keys.h"  #include "tamc_keys.h"
       INTEGER    isbyte  
       PARAMETER( isbyte = 4 )  
110  #else /* ALLOW_AUTODIFF_TAMC */  #else /* ALLOW_AUTODIFF_TAMC */
111        integer ikey        integer ikey
112  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 121  c     myTime - Current time in simulatio Line 124  c     myTime - Current time in simulatio
124    
125  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
126    
127    c Local constants
128    c     minusone, p0, p5, p25, p125, p0625
129    c     imin, imax, jmin, jmax  - array computation indices
130    
131          _RL        minusone
132          parameter( minusone=-1.0)
133          _KPP_RL    p0    , p5    , p25     , p125      , p0625
134          parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
135          integer    imin      , imax        , jmin      , jmax
136    #ifdef FRUGAL_KPP
137          parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )
138    #else
139          parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
140    #endif
141    
142  c Local arrays and variables  c Local arrays and variables
143  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
144  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 163  c     zRef   (nx,ny)       - Reference d
163  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
164  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
165    
166        _RS worka (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
167        _RS workb (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        integer work1 ( ibot:itop    , jbot:jtop                    )
168  #ifdef FRUGAL_KPP        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )
169        integer work1(sNx,sNy)        _KPP_RL work3 ( ibot:itop    , jbot:jtop                    )
170        _RS work2    (sNx,sNy)        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )
171        _RS ustar    (sNx,sNy)        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )
172        _RS bo       (sNx,sNy)        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )
173        _RS bosol    (sNx,sNy)        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )
174        _RS shsq     (sNx,sNy,Nr)        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )
175        _RS dVsq     (sNx,sNy,Nr)        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )
176        _RS dbloc    (sNx,sNy,Nr)        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )
177        _RS Ritop    (sNx,sNy,Nr)        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )
178        _RS vddiff   (sNx,sNy,0:Nrp1,mdiff)        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )
179        _RS ghat     (sNx,sNy,Nr)        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )
       _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)  
180  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
181        _RS z0       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )
182        _RS zRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )
183        _RS uRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )
184        _RS vRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )
185  #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 )  
186                
187        _RL     tempVar        _KPP_RL tempvar2
188        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, im1, ip1, jm1, jp1
189    
190  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
191        _RS     dBdz1, dBdz2, ustarX, ustarY        _KPP_RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
192  #endif  #endif
193    
194  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 233  c---------------------------------------
233        CALL STATEKPP(        CALL STATEKPP(
234       I       bi, bj, myThid       I       bi, bj, myThid
235       O     , work2, dbloc, Ritop       O     , work2, dbloc, Ritop
236  #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  
237       &     )       &     )
238        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
239    
 #ifdef KPP_SMOOTH_DBLOC  
 c     horizontally smooth dbloc with a 121 filter  
 c     (stored in ghat to save space)  
   
       DO k = 1, Nr  
          CALL SMOOTH_HORIZ_RS (  
      I        k, bi, bj,  
      I        dbloc(1-OLx,1-OLy,k),  
      O        ghat (1-OLx,1-OLy,k) )  
       ENDDO  
   
 #else /* KPP_SMOOTH_DBLOC */  
   
240        DO k = 1, Nr        DO k = 1, Nr
241  #ifdef FRUGAL_KPP           DO j = jbot, jtop
242           DO j = 1, sNy              DO i = ibot, itop
             DO i = 1, sNx  
 #else  
          DO j = 1-OLy, sNy+OLy  
             DO i = imin, imax  
 #endif  
243                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
244              ENDDO              ENDDO
245           ENDDO           ENDDO
246        ENDDO        ENDDO
247    
248    #ifdef KPP_SMOOTH_DBLOC
249    c     horizontally smooth dbloc with a 121 filter
250    c     smooth dbloc stored in ghat to save space
251    c     dbloc(k) is buoyancy gradientnote between k and k+1
252    c     levels therefore k+1 mask must be used
253    
254          DO k = 1, Nr-1
255             CALL KPP_SMOOTH_HORIZ (
256         I        k+1, bi, bj,
257         U        ghat (ibot,jbot,k) )
258          ENDDO
259    
260  #endif /* KPP_SMOOTH_DBLOC */  #endif /* KPP_SMOOTH_DBLOC */
261    
262  #ifdef KPP_SMOOTH_DENS  #ifdef KPP_SMOOTH_DENS
263  c     horizontally smooth density related quantities with 121 filters  c     horizontally smooth density related quantities with 121 filters
264        CALL SMOOTH_HORIZ_RS (        CALL KPP_SMOOTH_HORIZ (
265       I     k, bi, bj,       I     1, bi, bj,
266       I     work2,       U     work2 )
      O     work2 )  
267        DO k = 1, Nr        DO k = 1, Nr
268           CALL SMOOTH_HORIZ_RS (           CALL KPP_SMOOTH_HORIZ (
269       I        k, bi, bj,       I        k+1, bi, bj,
270       I        dbloc (1-OLx,1-OLy,k)  ,       U        dbloc (ibot,jbot,k) )
271       O        dbloc (1-OLx,1-OLy,k)   )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
272       I        k, bi, bj,       I        k, bi, bj,
273       I        Ritop (1-OLx,1-OLy,k)  ,       U        Ritop (ibot,jbot,k)  )
274       O        Ritop (1-OLx,1-OLy,k)   )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
275       I        k, bi, bj,       I        k, bi, bj,
276       I        vddiff(1-OLx,1-OLy,k,1),       U        vddiff(ibot,jbot,k,1) )
277       O        vddiff(1-OLx,1-OLy,k,1) )           CALL KPP_SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RS (  
278       I        k, bi, bj,       I        k, bi, bj,
279       I        vddiff(1-OLx,1-OLy,k,2),       U        vddiff(ibot,jbot,k,2) )
      O        vddiff(1-OLx,1-OLy,k,2) )  
280        ENDDO        ENDDO
281  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
282    
283        DO k = 1, Nr        DO k = 1, Nr
284  #ifdef FRUGAL_KPP           DO j = jbot, jtop
285           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  
286    
287  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
288  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 304  c     so that the subroutine "bldepth" w
304           END DO           END DO
305        END DO        END DO
306    
307    cph(
308    cph  this avoids a single or double recomp./call of statekpp
309    CADJ store work2              = comlev1_kpp, key = ikey
310    #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE
311    CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikey
312    CADJ store vddiff             = comlev1_kpp, key = ikey
313    #endif
314    cph)
315    
316  c------------------------------------------------------------------------  c------------------------------------------------------------------------
317  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
318  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
319  c     taux / rho = SurfaceTendencyU * delZ(1)                     (N/m^2)  c     taux / rho = SurfaceTendencyU * drF(1)                     (N/m^2)
320  c     tauy / rho = SurfaceTendencyV * delZ(1)                     (N/m^2)  c     tauy / rho = SurfaceTendencyV * drF(1)                     (N/m^2)
321  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
322  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*SurfaceTendencyT +
323  c                     beta *SurfaceTendencyS ) * delZ(1) / rho  (m^2/s^3)  c                     beta *SurfaceTendencyS ) * drF(1) / rho  (m^2/s^3)
324  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
325  c------------------------------------------------------------------------  c------------------------------------------------------------------------
326    
327  #ifdef FRUGAL_KPP  c initialize arrays to zero
328        DO j = 1, sNy        DO j = jbot, jtop
329           jp1 = j + 1           DO i = ibot, itop
330           DO i = 1, sNx              ustar(i,j) = p0
331  #else              bo   (I,J) = p0
332                bosol(I,J) = p0
333             END DO
334          END DO
335    
336        DO j = jmin, jmax        DO j = jmin, jmax
337         jp1 = j + 1         jp1 = j + 1
338         DO i = imin, imax         DO i = imin, imax
 #endif  
339          ip1 = i+1          ip1 = i+1
340          tempVar =          work3(i,j) =
341       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *
342       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +
343       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *
344       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))
345          if ( tempVar .lt. (epsln*epsln) ) then         END DO
346             ustar(i,j) = SQRT( epsln * p5 * delZ(1) )        END DO
347    cph(
348    CADJ store work3 = comlev1_kpp, key = ikey
349    cph)
350          DO j = jmin, jmax
351           jp1 = j + 1
352           DO i = imin, imax
353            ip1 = i+1
354    
355            if ( work3(i,j) .lt. (phepsi*phepsi) ) then
356               ustar(i,j) = SQRT( phepsi * p5 * drF(1) )
357          else          else
358             ustar(i,j) = SQRT( SQRT( tempVar ) * p5 * delZ(1) )             tempVar2 =  SQRT( work3(i,j) ) * p5 * drF(1)
359               ustar(i,j) = SQRT( tempVar2 )
360          endif          endif
361          bo(I,J) = - gravity *  
362            if ( .NOT. useSEAICE )
363         &       bo(I,J) = - gravity *
364       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +
365       &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)       &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)
366       &       ) *       &       ) *
367       &       delZ(1) / work2(I,J)       &       drF(1) / work2(I,J)
368          bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *  
369       &       delZ(1) / work2(I,J)  #ifdef ALLOW_SEAICE
370            if ( useSEAICE )
371         &       bo(I,J) = - gravity *
372         &       ( vddiff(I,J,1,1) * (SurfaceTendencyT(i,j,bi,bj)+
373         &                            SurfaceTendencyTice(i,j,bi,bj)) +
374         &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)
375         &       ) *
376         &       drF(1) / work2(I,J)
377    #endif ALLOW_SEAICE
378    
379            bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *
380         &       recip_Cp*recip_rhoConst*recip_dRf(1) *
381         &       drF(1) / work2(I,J)
382    
383         END DO         END DO
384        END DO        END DO
385    
386  #ifndef FRUGAL_KPP  cph(
387  c     set array edges to zero  CADJ store ustar = comlev1_kpp, key = ikey
388        DO j = jmin, jmax  cph)
          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  
389    
390  c------------------------------------------------------------------------  c------------------------------------------------------------------------
391  c     velocity shear  c     velocity shear
# Line 418  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k)) Line 396  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))
396  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
397  c------------------------------------------------------------------------  c------------------------------------------------------------------------
398    
399    c initialize arrays to zero
400          DO k = 1, Nr
401             DO j = jbot, jtop
402                DO i = ibot, itop
403                   shsq(i,j,k) = p0
404                   dVsq(i,j,k) = p0
405                END DO
406             END DO
407          END DO
408    
409  c     dVsq computation  c     dVsq computation
410    
411  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
# Line 428  c     thickness in the model.  First det Line 416  c     thickness in the model.  First det
416  c     Second zRef = espilon * hMix.  Third determine roughness length  c     Second zRef = espilon * hMix.  Third determine roughness length
417  c     scale z0.  Third estimate reference velocity.  c     scale z0.  Third estimate reference velocity.
418    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          jp1 = j + 1  
          DO i = 1, sNx  
 #else  
419        DO j = jmin, jmax        DO j = jmin, jmax
420           jp1 = j + 1           jp1 = j + 1
421           DO i = imin, imax           DO i = imin, imax
 #endif /* FRUGAL_KPP */  
422              ip1 = i + 1              ip1 = i + 1
423    
424  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 447  c     Linearly interpolate to find hMix.
447              ENDIF              ENDIF
448    
449  c     Compute roughness length scale z0 subject to 0 < z0  c     Compute roughness length scale z0 subject to 0 < z0
450                 tempVar = p5 * (                 tempVar1 = p5 * (
451       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
452       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
453       &              (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 456  c     Compute roughness length scale z0
456       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
457       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
458       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
459                 if ( tempVar .lt. (epsln*epsln) ) then                 if ( tempVar1 .lt. (epsln*epsln) ) then
460                    tempVar = epsln                    tempVar2 = epsln
461                 else                 else
462                    tempVar = SQRT ( tempVar )                    tempVar2 = SQRT ( tempVar1 )
463                 endif                 endif
464                 z0(i,j) = rF(2) *                 z0(i,j) = rF(2) *
465       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /
466       &                     ( rF(3) - rF(2) ) -       &                     ( rF(3) - rF(2) ) -
467       &                     tempVar * vonK /       &                     tempVar2 * vonK /
468       &                     MAX ( ustar(i,j), phepsi ) )       &                     MAX ( ustar(i,j), phepsi ) )
469                 z0(i,j) = MAX ( z0(i,j), phepsi )                 z0(i,j) = MAX ( z0(i,j), phepsi )
470    
# Line 500  c     Estimate reference velocity uRef a Line 482  c     Estimate reference velocity uRef a
482       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5
483                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +
484       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5
485                    tempVar = ustarX * ustarX + ustarY * ustarY                    tempVar1 = ustarX * ustarX + ustarY * ustarY
486                    if ( tempVar .lt. (epsln*epsln) ) then                    if ( tempVar1 .lt. (epsln*epsln) ) then
487                       tempVar = epsln                       tempVar2 = epsln
488                    else                    else
489                       tempVar = SQRT ( tempVar )                       tempVar2 = SQRT ( tempVar1 )
490                    endif                    endif
491                    tempVar = ustar(i,j) *                    tempVar2 = ustar(i,j) *
492       &                 ( LOG ( zRef(i,j) / rF(2) ) +       &                 ( LOG ( zRef(i,j) / rF(2) ) +
493       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /
494       &                 vonK / tempVar       &                 vonK / tempVar2
495                    uRef(i,j) = uRef(i,j) + ustarX * tempVar                    uRef(i,j) = uRef(i,j) + ustarX * tempVar2
496                    vRef(i,j) = vRef(i,j) + ustarY * tempVar                    vRef(i,j) = vRef(i,j) + ustarY * tempVar2
497                 ENDIF                 ENDIF
498    
499           END DO           END DO
500        END DO        END DO
501    
       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  
   
502        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  
503           DO j = jmin, jmax           DO j = jmin, jmax
504              jm1 = j - 1              jm1 = j - 1
505              jp1 = j + 1              jp1 = j + 1
506              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
507                 im1 = i - 1                 im1 = i - 1
508                 ip1 = i + 1                 ip1 = i + 1
509                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 604  c     Estimate reference velocity uRef a Line 541  c     Estimate reference velocity uRef a
541  #else /* KPP_ESTIMATE_UREF */  #else /* KPP_ESTIMATE_UREF */
542    
543        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  
544           DO j = jmin, jmax           DO j = jmin, jmax
545              jm1 = j - 1              jm1 = j - 1
546              jp1 = j + 1              jp1 = j + 1
547              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
548                 im1 = i - 1                 im1 = i - 1
549                 ip1 = i + 1                 ip1 = i + 1
550                 dVsq(i,j,k) = p5 * (                 dVsq(i,j,k) = p5 * (
# Line 654  c     Estimate reference velocity uRef a Line 584  c     Estimate reference velocity uRef a
584  c     shsq computation  c     shsq computation
585        DO k = 1, Nrm1        DO k = 1, Nrm1
586           kp1 = k + 1           kp1 = k + 1
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
587           DO j = jmin, jmax           DO j = jmin, jmax
588              jm1 = j - 1              jm1 = j - 1
589              jp1 = j + 1              jp1 = j + 1
590              DO i = imin, imax              DO i = imin, imax
 #endif /* FRUGAL_KPP */  
591                 im1 = i - 1                 im1 = i - 1
592                 ip1 = i + 1                 ip1 = i + 1
593                 shsq(i,j,k) = p5 * (                 shsq(i,j,k) = p5 * (
# Line 699  c     shsq computation Line 622  c     shsq computation
622           END DO           END DO
623        END DO        END DO
624    
625  c     shsq @ Nr computation  cph(
626  #ifdef FRUGAL_KPP  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE
627        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  
628  #endif  #endif
629    cph)
630    
631  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
632  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
633  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
634    
635  #ifdef FRUGAL_KPP        DO j = jbot, jtop
636        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  
637              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
638              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
639           END DO           END DO
640        END DO        END DO
641        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)
642        CALL KPPMIX (        CALL KPPMIX (
643       I       lri, work1, shsq, dVsq, ustar       I       mytime, mythid
644         I     , work1, shsq, dVsq, ustar
645       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
646       I     , ikey       I     , ikey
647       O     , vddiff       O     , vddiff
648       U     , ghat       U     , ghat
649       O     , hbl       O     , hbl )
      &        )  
650    
651        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
652    
       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  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
653  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
654  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  
655  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
656    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
657        DO j = jmin, jmax        DO j = jmin, jmax
658           DO i = imin, imax         DO i = imin, imax
659  #endif          DO k = 1, Nr
660              DO k = 1, Nr           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)
661  c KPPviscAz           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)
662                 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)
663                 tempVar = max( minKPPviscAz, tempVar )           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)
664                 KPPviscAz(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)          END DO
665  c KPPdiffKzS          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)
666                 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  
667        END DO        END DO
668  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
669        _EXCH_XYZ_R8(KPPviscAz  , myThid )        _EXCH_XYZ_R8(KPPviscAz  , myThid )
# Line 829  c KPPhbl: set to -zgrid(1) over land Line 675  c KPPhbl: set to -zgrid(1) over land
675    
676  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
677  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  
   
678        DO k = 1, Nr        DO k = 1, Nr
679           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
680       I        k, bi, bj,       I        k, bi, bj,
681       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) )  
682        END DO        END DO
683          _EXCH_XYZ_R8(KPPviscAz  , myThid )
684  #endif /* KPP_SMOOTH_VISC */  #endif /* KPP_SMOOTH_VISC */
685    
686  #ifdef KPP_SMOOTH_DIFF  #ifdef KPP_SMOOTH_DIFF
687  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  
   
688        DO k = 1, Nr        DO k = 1, Nr
689           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
690       I        k, bi, bj,       I        k, bi, bj,
691       I        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )
692       O        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )           CALL SMOOTH_HORIZ (
          CALL SMOOTH_HORIZ_RL (  
693       I        k, bi, bj,       I        k, bi, bj,
694       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) )  
695        END DO        END DO
696          _EXCH_XYZ_R8(KPPdiffKzS , myThid )
697          _EXCH_XYZ_R8(KPPdiffKzT , myThid )
698  #endif /* KPP_SMOOTH_DIFF */  #endif /* KPP_SMOOTH_DIFF */
699    
700    cph(
701    cph  crucial: this avoids full recomp./call of kppmix
702    CADJ store KPPhbl = comlev1_kpp, key = ikey
703    cph)
704    
705  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
706  C     the bottom of the mixing layer.  C     the bottom of the mixing layer.
# Line 866  C     the bottom of the mixing layer. Line 710  C     the bottom of the mixing layer.
710           ENDDO           ENDDO
711        ENDDO        ENDDO
712        CALL SWFRAC(        CALL SWFRAC(
713       I     (sNx+2*OLx)*(sNy+2*OLy), m1, worka,       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
714       O     workb )       I     mytime, mythid,
715         U     worka )
716        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
717           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
718              KPPfrac(i,j,bi,bj) = workb(i,j)              KPPfrac(i,j,bi,bj) = worka(i,j)
719           ENDDO           ENDDO
720        ENDDO        ENDDO
721    
722        ENDIF        ENDIF
723    
724  #endif ALLOW_KPP  #endif /* ALLOW_KPP */
725    
726          RETURN
727          END
728    
729          subroutine KPP_CALC_DUMMY(
730         I     bi, bj, myTime, myThid )
731    C     /==========================================================\
732    C     | SUBROUTINE KPP_CALC_DUMMY                                |
733    C     | o Compute all KPP fields defined in KPP.h                |
734    C     | o Dummy routine for TAMC
735    C     |==========================================================|
736    C     | This subroutine serves as an interface between MITGCMUV  |
737    C     | code and NCOM 1-D routines in kpp_routines.F             |
738    C     \==========================================================/
739          IMPLICIT NONE
740    
741    #include "SIZE.h"
742    #include "EEPARAMS.h"
743    #include "PARAMS.h"
744    #include "KPP.h"
745    #include "KPP_PARAMS.h"
746    #include "GRID.h"
747    
748    c Routine arguments
749    c     bi, bj - array indices on which to apply calculations
750    c     myTime - Current time in simulation
751    
752          INTEGER bi, bj
753          INTEGER myThid
754          _RL     myTime
755    
756    #ifdef ALLOW_KPP
757    
758    c Local constants
759          integer i, j, k
760    
761          DO j=1-OLy,sNy+OLy
762             DO i=1-OLx,sNx+OLx
763                KPPhbl (i,j,bi,bj) = 1.0
764                KPPfrac(i,j,bi,bj) = 0.0
765                DO k = 1,Nr
766                   KPPghat   (i,j,k,bi,bj) = 0.0
767                   KPPviscAz (i,j,k,bi,bj) = viscAz
768                   KPPdiffKzT(i,j,k,bi,bj) = diffKzT
769                   KPPdiffKzS(i,j,k,bi,bj) = diffKzS
770                ENDDO
771             ENDDO
772          ENDDO
773          
774    #endif
775        RETURN        RETURN
776        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22