/[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.11 by heimbach, Sat Jul 13 03:12:30 2002 UTC revision 1.41 by mlosch, Thu May 3 14:51:05 2007 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
5    
6    CBOP
7    C !ROUTINE: KPP_CALC
8    
9    C !INTERFACE: ==========================================================
10        subroutine KPP_CALC(        subroutine KPP_CALC(
11       I     bi, bj, myTime, myThid )       I     bi, bj, myTime, myThid )
12    
13    C !DESCRIPTION: \bv
14  C     /==========================================================\  C     /==========================================================\
15  C     | SUBROUTINE KPP_CALC                                      |  C     | SUBROUTINE KPP_CALC                                      |
16  C     | o Compute all KPP fields defined in KPP.h                |  C     | o Compute all KPP fields defined in KPP.h                |
# Line 84  c     KPPfrac     - Fraction of short-wa Line 90  c     KPPfrac     - Fraction of short-wa
90    
91  c--   KPP_CALC computes vertical viscosity and diffusivity for region  c--   KPP_CALC computes vertical viscosity and diffusivity for region
92  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
93  c     values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the  c     values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
94  c     region (-2:sNx+4,-2:sNy+4).  c     region (-2:sNx+4,-2:sNy+4).
95  c     Hence overlap region needs to be set OLx=4, OLy=4.  c     Hence overlap region needs to be set OLx=4, OLy=4.
96  c     When option FRUGAL_KPP is used, computation in overlap regions  c \ev
 c     is replaced with exchange calls hence reducing overlap requirements  
 c     to OLx=1, OLy=1.  
97    
98    C !USES: ===============================================================
99  #include "SIZE.h"  #include "SIZE.h"
100  #include "EEPARAMS.h"  #include "EEPARAMS.h"
101  #include "PARAMS.h"  #include "PARAMS.h"
# Line 99  c     to OLx=1, OLy=1. Line 104  c     to OLx=1, OLy=1.
104  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
105  #include "FFIELDS.h"  #include "FFIELDS.h"
106  #include "GRID.h"  #include "GRID.h"
107    #include "GAD.h"
108    #ifdef ALLOW_SHELFICE
109    # include "SHELFICE.h"
110    #endif /* ALLOW_SHELFICE */
111  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
112  #include "tamc.h"  #include "tamc.h"
113  #include "tamc_keys.h"  #include "tamc_keys.h"
114  #else /* ALLOW_AUTODIFF_TAMC */  #else /* ALLOW_AUTODIFF_TAMC */
115        integer ikey        integer ikppkey
116  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
117    
118        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
119        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
120    
121    C !INPUT PARAMETERS: ===================================================
122  c Routine arguments  c Routine arguments
123  c     bi, bj - array indices on which to apply calculations  c     bi, bj - array indices on which to apply calculations
124  c     myTime - Current time in simulation  c     myTime - Current time in simulation
# Line 120  c     myTime - Current time in simulatio Line 129  c     myTime - Current time in simulatio
129    
130  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
131    
132    C !LOCAL VARIABLES: ====================================================
133  c Local constants  c Local constants
134  c     minusone, p0, p5, p25, p125, p0625  c     minusone, p0, p5, p25, p125, p0625
135  c     imin, imax, jmin, jmax  - array computation indices  c     imin, imax, jmin, jmax  - array computation indices
136    
137        _RL        minusone        _RL        minusone
138        parameter( minusone=-1.0)        parameter( minusone=-1.0)
139        _KPP_RL    p0    , p5    , p25     , p125      , p0625        _RL        p0    , p5    , p25     , p125      , p0625
140        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
141        integer    imin      , imax        , jmin      , jmax        integer   imin      ,imax          ,jmin      ,jmax
142  #ifdef FRUGAL_KPP        parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
       parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )  
 #else  
       parameter( imin=-2   , imax=sNx+3  , jmin=-2   , jmax=sNy+3   )  
 #endif  
143    
144  c Local arrays and variables  c Local arrays and variables
145  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
# Line 149  c                            for ri_iwmi Line 155  c                            for ri_iwmi
155  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number
156  c                            at grid levels for bldepth  c                            at grid levels for bldepth
157  c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)  c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)
158  c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for temperature  (m^2/s)  c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for salt&tracers (m^2/s)
159  c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for salt&tracers (m^2/s)  c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for temperature  (m^2/s)
160  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)
161  c     hbl    (nx,ny)       - mixing layer depth                           (m)  c     hbl    (nx,ny)       - mixing layer depth                           (m)
162  c     kmtj   (nx,ny)       - maximum number of wet levels in each column  c     kmtj   (nx,ny)       - maximum number of wet levels in each column
# Line 159  c     zRef   (nx,ny)       - Reference d Line 165  c     zRef   (nx,ny)       - Reference d
165  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
166  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
167    
168        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )        integer work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy            )
169        integer work1 ( ibot:itop    , jbot:jtop                    )        _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
170        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )        _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
171        _KPP_RL work3 ( ibot:itop    , jbot:jtop                    )        _RL work3 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
172        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )        _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
173        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )        _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
174        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )        _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
175        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
176        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
177        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )        _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
178        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )        _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
179        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )        _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
180        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )        _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
181        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )        _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
182    cph(
183          _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
184          _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
185    cph)
186  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
187        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )        _RL z0    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
188        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )        _RL zRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
189        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )        _RL uRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
190        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )        _RL vRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
191  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
192                
193        _KPP_RL tempvar1, tempvar2        _RL tempvar2
194        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, km1, im1, ip1, jm1, jp1
195    
196  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
197        _KPP_RL dBdz1, dBdz2, ustarX, ustarY        _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
198  #endif  #endif
199    
200    #ifdef ALLOW_AUTODIFF_TAMC
201              act1 = bi - myBxLo(myThid)
202              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
203              act2 = bj - myByLo(myThid)
204              max2 = myByHi(myThid) - myByLo(myThid) + 1
205              act3 = myThid - 1
206              max3 = nTx*nTy
207              act4 = ikey_dynamics - 1
208              ikppkey = (act1 + 1) + act2*max1
209         &                      + act3*max1*max2
210         &                      + act4*max1*max2*max3
211    #endif /* ALLOW_AUTODIFF_TAMC */
212    CEOP
213    
214  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?
215        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
216       1     myTime .EQ. startTime ) THEN       1     .OR. myTime .EQ. startTime ) THEN
217                    
218  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
219  c     prepare input arrays for subroutine "kppmix" to compute  c     prepare input arrays for subroutine "kppmix" to compute
# Line 227  c--------------------------------------- Line 251  c---------------------------------------
251    
252        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)
253        CALL STATEKPP(        CALL STATEKPP(
254       I       bi, bj, myThid       O     work2, dbloc, Ritop,
255       O     , work2, dbloc, Ritop       O     TTALPHA, SSBETA,
256       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)       I     ikppkey, bi, bj, myThid )
      &     )  
257        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
258    
259        DO k = 1, Nr        DO k = 1, Nr
260           DO j = jbot, jtop           DO j = 1-OLy, sNy+OLy
261              DO i = ibot, itop              DO i = 1-OLx, sNx+OLx
262                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
263              ENDDO              ENDDO
264           ENDDO           ENDDO
# Line 248  c     dbloc(k) is buoyancy gradientnote Line 271  c     dbloc(k) is buoyancy gradientnote
271  c     levels therefore k+1 mask must be used  c     levels therefore k+1 mask must be used
272    
273        DO k = 1, Nr-1        DO k = 1, Nr-1
274           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
275       I        k+1, bi, bj,       I        k+1, bi, bj,
276       U        ghat (ibot,jbot,k) )       U        ghat (1-OLx,1-OLy,k),
277         I        myThid )
278        ENDDO        ENDDO
279    
280  #endif /* KPP_SMOOTH_DBLOC */  #endif /* KPP_SMOOTH_DBLOC */
281    
282  #ifdef KPP_SMOOTH_DENS  #ifdef KPP_SMOOTH_DENS
283  c     horizontally smooth density related quantities with 121 filters  c     horizontally smooth density related quantities with 121 filters
284        CALL KPP_SMOOTH_HORIZ (        CALL SMOOTH_HORIZ (
285       I     1, bi, bj,       I     1, bi, bj,
286       U     work2 )       U     work2,
287         I     myThid )
288        DO k = 1, Nr        DO k = 1, Nr
289           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
290       I        k+1, bi, bj,       I        k+1, bi, bj,
291       U        dbloc (ibot,jbot,k) )       U        dbloc (1-OLx,1-OLy,k),
292           CALL KPP_SMOOTH_HORIZ (       I        myThid )
293             CALL SMOOTH_HORIZ (
294       I        k, bi, bj,       I        k, bi, bj,
295       U        Ritop (ibot,jbot,k)  )       U        Ritop (1-OLx,1-OLy,k),
296           CALL KPP_SMOOTH_HORIZ (       I        myThid )
297             CALL SMOOTH_HORIZ (
298       I        k, bi, bj,       I        k, bi, bj,
299       U        vddiff(ibot,jbot,k,1) )       U        TTALPHA(1-OLx,1-OLy,k),
300           CALL KPP_SMOOTH_HORIZ (       I        myThid )
301             CALL SMOOTH_HORIZ (
302       I        k, bi, bj,       I        k, bi, bj,
303       U        vddiff(ibot,jbot,k,2) )       U        SSBETA(1-OLx,1-OLy,k),
304         I        myThid )
305        ENDDO        ENDDO
306  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
307    
308        DO k = 1, Nr        DO k = 1, Nr
309           DO j = jbot, jtop           km1 = max(1,k-1)
310              DO i = ibot, itop           DO j = 1-OLy, sNy+OLy
311                DO i = 1-OLx, sNx+OLx
312    
313  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
314  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
315                 dbloc(i,j,k) = dbloc(i,j,k) * pMask(i,j,k,bi,bj)                 dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
316                 ghat(i,j,k)  = ghat(i,j,k)  * pMask(i,j,k,bi,bj)       &              * maskC(i,j,km1,bi,bj)
317                 Ritop(i,j,k) = Ritop(i,j,k) * pMask(i,j,k,bi,bj)                 ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
318         &              * maskC(i,j,km1,bi,bj)
319                   Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
320         &              * maskC(i,j,km1,bi,bj)
321                 if(k.eq.nzmax(i,j,bi,bj)) then                 if(k.eq.nzmax(i,j,bi,bj)) then
322                    dbloc(i,j,k) = p0                    dbloc(i,j,k) = p0
323                    ghat(i,j,k)  = p0                    ghat(i,j,k)  = p0
# Line 302  c     so that the subroutine "bldepth" w Line 335  c     so that the subroutine "bldepth" w
335    
336  cph(  cph(
337  cph  this avoids a single or double recomp./call of statekpp  cph  this avoids a single or double recomp./call of statekpp
338  CADJ store work2              = comlev1_kpp, key = ikey  CADJ store work2              = comlev1_kpp, key = ikppkey
339  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
340  CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikey  CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey
341  CADJ store vddiff             = comlev1_kpp, key = ikey  CADJ store vddiff             = comlev1_kpp, key = ikppkey
342    CADJ store TTALPHA, SSBETA    = comlev1_kpp, key = ikppkey
343  #endif  #endif
344  cph)  cph)
345    
346    CML#ifdef ALLOW_SHELFICE
347    CMLC     For the pbl parameterisation to work underneath the ice shelves
348    CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
349    CMLC     and indexing problems make this part of the code not work
350    CMLC     underneath the ice shelves and the following lines are only here
351    CMLC     to remind me that this still needs to be sorted out.
352    CML      shelfIceFac = 0. _d 0
353    CML      IF ( useShelfIce ) selfIceFac = 1. _d 0
354    CML      DO j = jmin, jmax
355    CML       DO i = imin, imax
356    CML        surfForcT = surfaceForcingT(i,j,bi,bj)
357    CML     &       + shelficeForcingT(i,j,bi,bj) * shelfIceFac
358    CML        surfForcS = surfaceForcingS(i,j,bi,bj)
359    CML     &       + shelficeForcingS(i,j,bi,bj) * shelfIceFac
360    CML       END DO
361    CML      END DO
362    CML#endif /* ALLOW_SHELFICE */
363    
364  c------------------------------------------------------------------------  c------------------------------------------------------------------------
365  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
366  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
367  c     taux / rho = SurfaceTendencyU * drF(1)                     (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
368  c     tauy / rho = SurfaceTendencyV * drF(1)                     (N/m^2)  c     tauy / rho = surfaceForcingV                               (N/m^2)
369  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
370  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*surfaceForcingT +
371  c                     beta *SurfaceTendencyS ) * drF(1) / rho  (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
372  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
373  c------------------------------------------------------------------------  c------------------------------------------------------------------------
   
 c initialize arrays to zero  
       DO j = jbot, jtop  
          DO i = ibot, itop  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
   
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
         work3(i,j) =  
      &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *  
      &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +  
      &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *  
      &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))  
        END DO  
       END DO  
 cph(  
 CADJ store work3 = comlev1_kpp, key = ikey  
 cph)  
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
         if ( work3(i,j) .lt. (phepsi*phepsi) ) then  
            ustar(i,j) = SQRT( phepsi * p5 * drF(1) )  
         else  
            tempVar2 =  SQRT( work3(i,j) ) * p5 * drF(1)  
            ustar(i,j) = SQRT( tempVar2 )  
         endif  
         bo(I,J) = - gravity *  
      &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +  
      &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)  
      &       ) *  
      &       drF(1) / work2(I,J)  
         bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *  
      &       recip_Cp*recip_rhoNil*recip_dRf(1) *  
      &       drF(1) / work2(I,J)  
        END DO  
       END DO  
   
 cph(  
 CADJ store ustar = comlev1_kpp, key = ikey  
 cph)  
   
 c------------------------------------------------------------------------  
374  c     velocity shear  c     velocity shear
375  c     --------------  c     --------------
376  c     Get velocity shear squared, averaged from "u,v-grid"  c     Get velocity shear squared, averaged from "u,v-grid"
377  c     onto "t-grid" (in (m/s)**2):  c     onto "t-grid" (in (m/s)**2):
378  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
379  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
380    c    
381    c     note: Vref can depend on the surface fluxes that is why we compute
382    c     dVsq in the subroutine that does the surface related stuff
383    c     (admittedly this is a bit messy)
384  c------------------------------------------------------------------------  c------------------------------------------------------------------------
385          
386          CALL KPP_FORCING_SURF(
387         I     work2, surfaceForcingU, surfaceForcingV,
388         I     surfaceForcingT, surfaceForcingS, surfaceForcingTice,
389         I     Qsw, ttalpha, ssbeta,  
390         O     ustar, bo, bosol, dVsq,
391         I     ikppkey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
392    
393    CMLcph(
394    CMLCADJ store ustar = comlev1_kpp, key = ikppkey
395    CMLcph)
396    
397  c initialize arrays to zero  c initialize arrays to zero
398        DO k = 1, Nr        DO k = 1, Nr
399           DO j = jbot, jtop         DO j = 1-OLy, sNy+OLy
400              DO i = ibot, itop          DO i = 1-OLx, sNx+OLx
401                 shsq(i,j,k) = p0           shsq(i,j,k) = p0
402                 dVsq(i,j,k) = p0          END DO
403              END DO         END DO
          END DO  
       END DO  
   
 c     dVsq computation  
   
 #ifdef KPP_ESTIMATE_UREF  
   
 c     Get rid of vertical resolution dependence of dVsq term by  
 c     estimating a surface velocity that is independent of first level  
 c     thickness in the model.  First determine mixed layer depth hMix.  
 c     Second zRef = espilon * hMix.  Third determine roughness length  
 c     scale z0.  Third estimate reference velocity.  
   
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
             ip1 = i + 1  
   
 c     Determine mixed layer depth hMix as the shallowest depth at which  
 c     dB/dz exceeds 5.2e-5 s^-2.  
             work1(i,j) = nzmax(i,j,bi,bj)  
             DO k = 1, Nr  
                IF ( k .LT. nzmax(i,j,bi,bj) .AND.  
      &              dbloc(i,j,k) / drC(k+1) .GT. dB_dz )  
      &              work1(i,j) = k  
             END DO  
   
 c     Linearly interpolate to find hMix.  
             k = work1(i,j)  
             IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN  
                zRef(i,j) = p0  
             ELSEIF ( k .EQ. 1) THEN  
                dBdz2 = dbloc(i,j,1) / drC(2)  
                zRef(i,j) = drF(1) * dB_dz / dBdz2  
             ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN  
                dBdz1 = dbloc(i,j,k-1) / drC(k  )  
                dBdz2 = dbloc(i,j,k  ) / drC(k+1)  
                zRef(i,j) = rF(k) + drF(k) * (dB_dz - dBdz1) /  
      &                     MAX ( phepsi, dBdz2 - dBdz1 )  
             ELSE  
                zRef(i,j) = rF(k+1)  
             ENDIF  
   
 c     Compute roughness length scale z0 subject to 0 < z0  
                tempVar1 = p5 * (  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) +  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )  
                if ( tempVar1 .lt. (epsln*epsln) ) then  
                   tempVar2 = epsln  
                else  
                   tempVar2 = SQRT ( tempVar1 )  
                endif  
                z0(i,j) = rF(2) *  
      &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /  
      &                     ( rF(3) - rF(2) ) -  
      &                     tempVar2 * vonK /  
      &                     MAX ( ustar(i,j), phepsi ) )  
                z0(i,j) = MAX ( z0(i,j), phepsi )  
   
 c     zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1)  
                zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) )  
                zRef(i,j) = MIN ( zRef(i,j), drF(1) )  
   
 c     Estimate reference velocity uRef and vRef.  
                uRef(i,j) = p5 *  
      &                     ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) )  
                vRef(i,j) = p5 *  
      &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )  
                IF ( zRef(i,j) .LT. drF(1) ) THEN  
                   ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +  
      &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5  
                   ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +  
      &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5  
                   tempVar1 = ustarX * ustarX + ustarY * ustarY  
                   if ( tempVar1 .lt. (epsln*epsln) ) then  
                      tempVar2 = epsln  
                   else  
                      tempVar2 = SQRT ( tempVar1 )  
                   endif  
                   tempVar2 = ustar(i,j) *  
      &                 ( LOG ( zRef(i,j) / rF(2) ) +  
      &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /  
      &                 vonK / tempVar2  
                   uRef(i,j) = uRef(i,j) + ustarX * tempVar2  
                   vRef(i,j) = vRef(i,j) + ustarY * tempVar2  
                ENDIF  
   
          END DO  
       END DO  
   
       DO k = 1, Nr  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
                im1 = i - 1  
                ip1 = i + 1  
                dVsq(i,j,k) = p5 * (  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
       END DO  
   
 #else /* KPP_ESTIMATE_UREF */  
   
       DO k = 1, Nr  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
                im1 = i - 1  
                ip1 = i + 1  
                dVsq(i,j,k) = p5 * (  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) *  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) +  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) *  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) *  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) +  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) *  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) *  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) +  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) *  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) +  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) *  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) +  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) *  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) +  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) *  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) +  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) *  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) +  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) *  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) *  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
404        END DO        END DO
405    
 #endif /* KPP_ESTIMATE_UREF */  
   
406  c     shsq computation  c     shsq computation
407        DO k = 1, Nrm1        DO k = 1, Nrm1
408           kp1 = k + 1         kp1 = k + 1
409           DO j = jmin, jmax         DO j = jmin, jmax
410              jm1 = j - 1          jm1 = j - 1
411              jp1 = j + 1          jp1 = j + 1
412              DO i = imin, imax          DO i = imin, imax
413                 im1 = i - 1           im1 = i - 1
414                 ip1 = i + 1           ip1 = i + 1
415                 shsq(i,j,k) = p5 * (           shsq(i,j,k) = p5 * (
416       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *       $        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
417       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +       $        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
418       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *       $        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
419       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +       $        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
420       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *       $        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
421       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +       $        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
422       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *       $        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
423       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )       $        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
424  #ifdef KPP_SMOOTH_SHSQ  #ifdef KPP_SMOOTH_SHSQ
425                 shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (           shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
426       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *       $        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *
427       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +       $        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +
428       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *       $        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
429       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +       $        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
430       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *       $        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *
431       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +       $        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +
432       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *       $        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
433       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +       $        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
434       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *       $        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *
435       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +       $        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +
436       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *       $        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
437       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +       $        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
438       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *       $        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *
439       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +       $        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +
440       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *       $        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
441       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )       $        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
442  #endif  #endif
443              END DO          END DO
444           END DO         END DO
445        END DO        END DO
446    
447  cph(  cph(
448  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
449  CADJ store dvsq, shsq = comlev1_kpp, key = ikey  CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
450  #endif  #endif
451  cph)  cph)
452    
# Line 613  c--------------------------------------- Line 454  c---------------------------------------
454  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
455  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
456    
457        DO j = jbot, jtop  c     precompute background vertical diffusivities, which are needed for
458           DO i = ibot, itop  c     matching diffusivities at bottom of KPP PBL
459          CALL CALC_3D_DIFFUSIVITY(
460         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
461         I        GAD_SALINITY, .FALSE., .FALSE.,
462         O        KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
463         I        myThid)
464          CALL CALC_3D_DIFFUSIVITY(
465         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
466         I        GAD_TEMPERATURE, .FALSE., .FALSE.,
467         O        KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
468         I        myThid)
469    
470          DO j = 1-OLy, sNy+OLy
471             DO i = 1-OLx, sNx+OLx
472              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
473              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
474           END DO           END DO
475        END DO        END DO
476        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)
477        CALL KPPMIX (        CALL KPPMIX (
478       I       mytime, mythid       I       work1, shsq, dVsq, ustar
479       I     , work1, shsq, dVsq, ustar       I     , maskC(1-Olx,1-Oly,1,bi,bj)
480       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
481       I     , ikey       I     , KPPdiffKzS(1-Olx,1-Oly,1,bi,bj)
482         I     , KPPdiffKzT(1-Olx,1-Oly,1,bi,bj)
483         I     , ikppkey
484       O     , vddiff       O     , vddiff
485       U     , ghat       U     , ghat
486       O     , hbl )       O     , hbl
487         I     , mytime, mythid )
488        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
489    
490  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 638  c--------------------------------------- Line 494  c---------------------------------------
494        DO j = jmin, jmax        DO j = jmin, jmax
495         DO i = imin, imax         DO i = imin, imax
496          DO k = 1, Nr          DO k = 1, Nr
497           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)           km1 = max(1,k-1)
498           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
499           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)       &        * maskC(i,j,km1,bi,bj)
500           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
501         &        * maskC(i,j,km1,bi,bj)
502             KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
503         &        * maskC(i,j,km1,bi,bj)
504             KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
505         &        * maskC(i,j,km1,bi,bj)
506          END DO          END DO
507          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)          k = 1
508    #ifdef ALLOW_SHELFICE
509            if ( useShelfIce ) k = kTopC(i,j,bi,bj)
510    #endif /* ALLOW_SHELFICE */
511            KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,k,bi,bj)
512    
513         END DO         END DO
514        END DO        END DO
 #ifdef FRUGAL_KPP  
       _EXCH_XYZ_R8(KPPviscAz  , myThid )  
       _EXCH_XYZ_R8(KPPdiffKzS , myThid )  
       _EXCH_XYZ_R8(KPPdiffKzT , myThid )  
       _EXCH_XYZ_R8(KPPghat    , myThid )  
       _EXCH_XY_R8 (KPPhbl     , myThid )  
 #endif  
515    
516  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
517  c     horizontal smoothing of vertical viscosity  c     horizontal smoothing of vertical viscosity
518        DO k = 1, Nr        DO k = 1, Nr
519           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
520       I        k, bi, bj,       I        k, bi, bj,
521       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj) )       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj),
522         I        myThid )
523        END DO        END DO
524        _EXCH_XYZ_R8(KPPviscAz  , myThid )        _EXCH_XYZ_R8(KPPviscAz  , myThid )
525  #endif /* KPP_SMOOTH_VISC */  #endif /* KPP_SMOOTH_VISC */
# Line 669  c     horizontal smoothing of vertical d Line 529  c     horizontal smoothing of vertical d
529        DO k = 1, Nr        DO k = 1, Nr
530           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
531       I        k, bi, bj,       I        k, bi, bj,
532       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
533         I        myThid )
534           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
535       I        k, bi, bj,       I        k, bi, bj,
536       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) )       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
537         I        myThid )
538        END DO        END DO
539        _EXCH_XYZ_R8(KPPdiffKzS , myThid )        _EXCH_XYZ_R8(KPPdiffKzS , myThid )
540        _EXCH_XYZ_R8(KPPdiffKzT , myThid )        _EXCH_XYZ_R8(KPPdiffKzT , myThid )
# Line 680  c     horizontal smoothing of vertical d Line 542  c     horizontal smoothing of vertical d
542    
543  cph(  cph(
544  cph  crucial: this avoids full recomp./call of kppmix  cph  crucial: this avoids full recomp./call of kppmix
545  CADJ store KPPhbl = comlev1_kpp, key = ikey  CADJ store KPPhbl = comlev1_kpp, key = ikppkey
546  cph)  cph)
547    
548  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
# Line 725  C     \================================= Line 587  C     \=================================
587  #include "KPP.h"  #include "KPP.h"
588  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
589  #include "GRID.h"  #include "GRID.h"
590    #include "GAD.h"
591    
592  c Routine arguments  c Routine arguments
593  c     bi, bj - array indices on which to apply calculations  c     bi, bj - array indices on which to apply calculations
# Line 745  c Local constants Line 608  c Local constants
608              KPPfrac(i,j,bi,bj) = 0.0              KPPfrac(i,j,bi,bj) = 0.0
609              DO k = 1,Nr              DO k = 1,Nr
610                 KPPghat   (i,j,k,bi,bj) = 0.0                 KPPghat   (i,j,k,bi,bj) = 0.0
611                 KPPviscAz (i,j,k,bi,bj) = viscAz                 KPPviscAz (i,j,k,bi,bj) = viscAr
                KPPdiffKzT(i,j,k,bi,bj) = diffKzT  
                KPPdiffKzS(i,j,k,bi,bj) = diffKzS  
612              ENDDO              ENDDO
613           ENDDO           ENDDO
614        ENDDO        ENDDO
615          
616          CALL CALC_3D_DIFFUSIVITY(
617         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
618         I     GAD_SALINITY, .FALSE., .FALSE.,
619         O     KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
620         I     myThid)
621          CALL CALC_3D_DIFFUSIVITY(
622         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
623         I     GAD_TEMPERATURE, .FALSE., .FALSE.,
624         O     KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
625         I     myThid)
626    
627  #endif  #endif
628        RETURN        RETURN
629        END        END
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22