pkg/kpp/kpp_routines.F

2fa42a6013 Alis*0001 #include "KPP_OPTIONS.h"
853c5e0e2c Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
1f89baba18 Patr*0005 #ifdef ALLOW_SALT_PLUME
                0006 #include "SALT_PLUME_OPTIONS.h"
                0007 #endif
853c5e0e2c Jean*0008 #if (defined ALLOW_AUTODIFF_TAMC) && (defined KPP_AUTODIFF_EXCESSIVE_STORE)
                0009 # define KPP_AUTODIFF_MORE_STORE
                0010 #endif
2fa42a6013 Alis*0011
                0012 C-- File kpp_routines.F: subroutines needed to implement
                0013 C--                      KPP vertical mixing scheme
                0014 C--  Contents
                0015 C--  o KPPMIX      - Main driver and interface routine.
                0016 C--  o BLDEPTH     - Determine oceanic planetary boundary layer depth.
                0017 C--  o WSCALE      - Compute turbulent velocity scales.
                0018 C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
                0019 C--  o Z121        - Apply 121 vertical smoothing.
fe66051ebd Dimi*0020 C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
2fa42a6013 Alis*0021 C--  o BLMIX       - Boundary layer mixing coefficients.
                0022 C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
                0023 C--  o STATEKPP    - Compute buoyancy-related input arrays.
e750a5e49e Mart*0024 C--  o KPP_DOUBLEDIFF - Compute double diffusive contribution to diffusivities
2fa42a6013 Alis*0025
956c0a5824 Patr*0026 c***********************************************************************
2fa42a6013 Alis*0027
                0028       SUBROUTINE KPPMIX (
edb6656069 Mart*0029      I       kmtj, shsq, dvsq, ustar, msk,
                0030      I       bo, bosol,
63ceaaa79c Dimi*0031 #ifdef ALLOW_SALT_PLUME
edb6656069 Mart*0032      I       boplume, SPDepth,
1f89baba18 Patr*0033 #ifdef SALT_PLUME_SPLIT_BASIN
edb6656069 Mart*0034      I       lon, lat,
1f89baba18 Patr*0035 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0036 #endif /* ALLOW_SALT_PLUME */
edb6656069 Mart*0037      I       dbloc, Ritop, coriol,
                0038      I       diffusKzS, diffusKzT,
                0039      I       ikey,
                0040      O       diffus,
                0041      U       ghat,
                0042      O       hbl,
                0043      I       bi, bj, myTime, myIter, myThid )
2fa42a6013 Alis*0044
956c0a5824 Patr*0045 c-----------------------------------------------------------------------
2fa42a6013 Alis*0046 c
                0047 c     Main driver subroutine for kpp vertical mixing scheme and
                0048 c     interface to greater ocean model
                0049 c
                0050 c     written  by: bill large,    june  6, 1994
                0051 c     modified by: jan morzel,    june 30, 1994
                0052 c                  bill large,  august 11, 1994
                0053 c                  bill large, january 25, 1995 : "dVsq" and 1d code
                0054 c                  detlef stammer,  august 1997 : for use with MIT GCM Classic
                0055 c                  d. menemenlis,     june 1998 : for use with MIT GCM UV
                0056 c
                0057 c-----------------------------------------------------------------------
                0058
                0059       IMPLICIT NONE
                0060
d2175d6909 Jean*0061 #include "SIZE.h"
                0062 #include "EEPARAMS.h"
                0063 #include "PARAMS.h"
                0064 #include "KPP_PARAMS.h"
853c5e0e2c Jean*0065 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*0066 # include "tamc.h"
                0067 #endif
d2175d6909 Jean*0068
2fa42a6013 Alis*0069 c input
15fd724cec Dimi*0070 c   bi, bj :: Array indices on which to apply calculations
eeda3a3aa1 Jean*0071 c   myTime :: Current time in simulation
                0072 c   myIter :: Current iteration number in simulation
                0073 c   myThid :: My Thread Id. number
edb6656069 Mart*0074 C     ikey :: tape key TAF-AD simulations (depends on tiles)
059d9fc14f Dimi*0075 c     kmtj   (imt)     - number of vertical layers on this row
25c3477f99 Mart*0076 c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
059d9fc14f Dimi*0077 c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
                0078 c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
                0079 c     ustar  (imt)     - surface friction velocity                        (m/s)
                0080 c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
                0081 c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
1f89baba18 Patr*0082 c     boplume(imt,Nrp1)- haline buoyancy forcing                      (m^2/s^3)
059d9fc14f Dimi*0083 c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
                0084 c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
                0085 c                          stored in ghat to save space
                0086 c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
                0087 c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
                0088 c     coriol (imt)     - Coriolis parameter                               (1/s)
                0089 c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
                0090 c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
2fa42a6013 Alis*0091 c     note: there is a conversion from 2-D to 1-D for input output variables,
                0092 c           e.g., hbl(sNx,sNy) -> hbl(imt),
                0093 c           where hbl(i,j) -> hbl((j-1)*sNx+i)
15fd724cec Dimi*0094       INTEGER bi, bj
eeda3a3aa1 Jean*0095       _RL     myTime
2b4c90c108 Mart*0096       INTEGER myIter
                0097       INTEGER myThid
                0098       INTEGER kmtj (imt   )
fe66051ebd Dimi*0099       _RL shsq     (imt,Nr)
                0100       _RL dvsq     (imt,Nr)
                0101       _RL ustar    (imt   )
                0102       _RL bo       (imt   )
                0103       _RL bosol    (imt   )
63ceaaa79c Dimi*0104 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0105       _RL boplume  (imt,Nrp1)
63ceaaa79c Dimi*0106       _RL SPDepth  (imt   )
1f89baba18 Patr*0107 #ifdef SALT_PLUME_SPLIT_BASIN
                0108       _RL lon  (imt   )
                0109       _RL lat  (imt   )
                0110 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0111 #endif /* ALLOW_SALT_PLUME */
fe66051ebd Dimi*0112       _RL dbloc    (imt,Nr)
                0113       _RL Ritop    (imt,Nr)
                0114       _RL coriol   (imt   )
25c3477f99 Mart*0115       _RS msk      (imt   )
fe66051ebd Dimi*0116       _RL diffusKzS(imt,Nr)
                0117       _RL diffusKzT(imt,Nr)
2fa42a6013 Alis*0118
edb6656069 Mart*0119       INTEGER ikey
2fa42a6013 Alis*0120
                0121 c output
                0122 c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)
                0123 c     diffus (imt,2)  - vertical scalar diffusivity                     (m^2/s)
                0124 c     diffus (imt,3)  - vertical temperature diffusivity                (m^2/s)
                0125 c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
                0126 c     hbl    (imt)    - mixing layer depth                                  (m)
                0127
fe66051ebd Dimi*0128       _RL diffus(imt,0:Nrp1,mdiff)
                0129       _RL ghat  (imt,Nr)
                0130       _RL hbl   (imt)
2fa42a6013 Alis*0131
                0132 #ifdef ALLOW_KPP
                0133
                0134 c local
                0135 c     kbl    (imt         ) - index of first grid level below hbl
                0136 c     bfsfc  (imt         ) - surface buoyancy forcing                (m^2/s^3)
                0137 c     casea  (imt         ) - 1 in case A; 0 in case B
                0138 c     stable (imt         ) - 1 in stable forcing; 0 if unstable
                0139 c     dkm1   (imt,   mdiff) - boundary layer diffusivity at kbl-1 level
                0140 c     blmc   (imt,Nr,mdiff) - boundary layer mixing coefficients
                0141 c     sigma  (imt         ) - normalized depth (d / hbl)
                0142 c     Rib    (imt,Nr      ) - bulk Richardson number
                0143
2b4c90c108 Mart*0144       INTEGER kbl(imt         )
fe66051ebd Dimi*0145       _RL bfsfc  (imt         )
                0146       _RL casea  (imt         )
                0147       _RL stable (imt         )
                0148       _RL dkm1   (imt,   mdiff)
                0149       _RL blmc   (imt,Nr,mdiff)
                0150       _RL sigma  (imt         )
                0151       _RL Rib    (imt,Nr      )
2fa42a6013 Alis*0152
2b4c90c108 Mart*0153       INTEGER i, k, md
2fa42a6013 Alis*0154
                0155 c-----------------------------------------------------------------------
                0156 c compute interior mixing coefficients everywhere, due to constant
                0157 c internal wave activity, static instability, and local shear
                0158 c instability.
                0159 c (ghat is temporary storage for horizontally smoothed dbloc)
                0160 c-----------------------------------------------------------------------
                0161
6d45c3b90d Patr*0162 cph(
                0163 cph these storings avoid recomp. of Ri_iwmix
853c5e0e2c Jean*0164 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0165 CADJ STORE ghat  = comlev1_kpp, key=ikey, kind=isbyte
                0166 CADJ STORE dbloc = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*0167 #endif
6d45c3b90d Patr*0168 cph)
2b4c90c108 Mart*0169       CALL Ri_iwmix (
edb6656069 Mart*0170      I       kmtj, shsq, dbloc, ghat,
                0171      I       diffusKzS, diffusKzT,
                0172      I       ikey,
                0173      O       diffus, myThid )
2fa42a6013 Alis*0174
6d45c3b90d Patr*0175 cph(
                0176 cph these storings avoid recomp. of Ri_iwmix
                0177 cph DESPITE TAFs 'not necessary' warning!
853c5e0e2c Jean*0178 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0179 CADJ STORE dbloc  = comlev1_kpp, key=ikey, kind=isbyte
                0180 CADJ STORE shsq   = comlev1_kpp, key=ikey, kind=isbyte
                0181 CADJ STORE ghat   = comlev1_kpp, key=ikey, kind=isbyte
                0182 CADJ STORE diffus = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*0183 #endif
6d45c3b90d Patr*0184 cph)
                0185
2fa42a6013 Alis*0186 c-----------------------------------------------------------------------
                0187 c set seafloor values to zero and fill extra "Nrp1" coefficients
                0188 c for blmix
                0189 c-----------------------------------------------------------------------
                0190
2b4c90c108 Mart*0191       DO md = 1, mdiff
                0192        DO k=1,Nrp1
                0193         DO i = 1,imt
                0194          IF (k.GE.kmtj(i)) diffus(i,k,md) = 0.0
                0195         ENDDO
                0196        ENDDO
                0197       ENDDO
2fa42a6013 Alis*0198
                0199 c-----------------------------------------------------------------------
                0200 c compute boundary layer mixing coefficients:
                0201 c
                0202 c diagnose the new boundary layer depth
                0203 c-----------------------------------------------------------------------
                0204
2b4c90c108 Mart*0205       CALL bldepth (
edb6656069 Mart*0206      I       kmtj,
                0207      I       dvsq, dbloc, Ritop, ustar, bo, bosol,
63ceaaa79c Dimi*0208 #ifdef ALLOW_SALT_PLUME
edb6656069 Mart*0209      I       boplume, SPDepth,
1f89baba18 Patr*0210 #ifdef SALT_PLUME_SPLIT_BASIN
edb6656069 Mart*0211      I       lon, lat,
1f89baba18 Patr*0212 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0213 #endif /* ALLOW_SALT_PLUME */
edb6656069 Mart*0214      I       coriol,
                0215      I       ikey,
                0216      O       hbl, bfsfc, stable, casea, kbl, Rib, sigma,
                0217      I       bi, bj, myTime, myIter, myThid )
2fa42a6013 Alis*0218
853c5e0e2c Jean*0219 #ifdef ALLOW_AUTODIFF_TAMC
a9d2e4c565 Jean*0220 CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
edb6656069 Mart*0221 CADJ &     key=ikey, kind=isbyte
853c5e0e2c Jean*0222 #endif
2fa42a6013 Alis*0223
                0224 c-----------------------------------------------------------------------
                0225 c compute boundary layer diffusivities
                0226 c-----------------------------------------------------------------------
                0227
2b4c90c108 Mart*0228       CALL blmix (
edb6656069 Mart*0229      I       ustar, bfsfc, hbl, stable, casea, diffus, kbl,
                0230      O       dkm1, blmc, ghat, sigma, ikey,
                0231      I       myThid )
6d45c3b90d Patr*0232 cph(
853c5e0e2c Jean*0233 #ifdef ALLOW_AUTODIFF_TAMC
a9d2e4c565 Jean*0234 CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
edb6656069 Mart*0235 CADJ &     key=ikey, kind=isbyte
a9d2e4c565 Jean*0236 CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
edb6656069 Mart*0237 CADJ &     key=ikey, kind=isbyte
853c5e0e2c Jean*0238 #endif
6d45c3b90d Patr*0239 cph)
2fa42a6013 Alis*0240
                0241 c-----------------------------------------------------------------------
                0242 c enhance diffusivity at interface kbl - 1
                0243 c-----------------------------------------------------------------------
                0244
2b4c90c108 Mart*0245       CALL enhance (
edb6656069 Mart*0246      I       dkm1, hbl, kbl, diffus, casea,
                0247      U       ghat,
                0248      O       blmc,
                0249      I       myThid )
2fa42a6013 Alis*0250
6d45c3b90d Patr*0251 cph(
                0252 cph avoids recomp. of enhance
853c5e0e2c Jean*0253 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0254 CADJ STORE blmc = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*0255 #endif
6d45c3b90d Patr*0256 cph)
                0257
2fa42a6013 Alis*0258 c-----------------------------------------------------------------------
                0259 c combine interior and boundary layer coefficients and nonlocal term
6060ec2938 Dimi*0260 c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
04522666de Ed H*0261 c (< 1 level), diffusivity blmc can become negative.  The max-s below
6060ec2938 Dimi*0262 c are a hack until this problem is properly diagnosed and fixed.
2fa42a6013 Alis*0263 c-----------------------------------------------------------------------
2b4c90c108 Mart*0264       DO k = 1, Nr
                0265          DO i = 1, imt
                0266             IF (k .LT. kbl(i)) THEN
25c3477f99 Mart*0267 #ifdef ALLOW_SHELFICE
                0268 C     when there is shelfice on top (msk(i)=0), reset the boundary layer
                0269 C     mixing coefficients blmc to pure Ri-number based mixing
2b4c90c108 Mart*0270                blmc(i,k,1) = MAX ( blmc(i,k,1)*msk(i),
25c3477f99 Mart*0271      &              diffus(i,k,1) )
2b4c90c108 Mart*0272                blmc(i,k,2) = MAX ( blmc(i,k,2)*msk(i),
25c3477f99 Mart*0273      &              diffus(i,k,2) )
2b4c90c108 Mart*0274                blmc(i,k,3) = MAX ( blmc(i,k,3)*msk(i),
25c3477f99 Mart*0275      &              diffus(i,k,3) )
                0276 #endif /* not ALLOW_SHELFICE */
2b4c90c108 Mart*0277                diffus(i,k,1) = MAX ( blmc(i,k,1), viscArNr(1) )
                0278                diffus(i,k,2) = MAX ( blmc(i,k,2), diffusKzS(i,Nr) )
                0279                diffus(i,k,3) = MAX ( blmc(i,k,3), diffusKzT(i,Nr) )
                0280             ELSE
b5ecce171d Davi*0281                ghat(i,k) = 0. _d 0
2b4c90c108 Mart*0282             ENDIF
                0283          ENDDO
                0284       ENDDO
2fa42a6013 Alis*0285
                0286 #endif /* ALLOW_KPP */
                0287
2b4c90c108 Mart*0288       RETURN
                0289       END
2fa42a6013 Alis*0290
                0291 c*************************************************************************
                0292
                0293       subroutine bldepth (
edb6656069 Mart*0294      I       kmtj,
                0295      I       dvsq, dbloc, Ritop, ustar, bo, bosol,
63ceaaa79c Dimi*0296 #ifdef ALLOW_SALT_PLUME
edb6656069 Mart*0297      I       boplume, SPDepth,
1f89baba18 Patr*0298 #ifdef SALT_PLUME_SPLIT_BASIN
edb6656069 Mart*0299      I       lon, lat,
1f89baba18 Patr*0300 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0301 #endif /* ALLOW_SALT_PLUME */
edb6656069 Mart*0302      I       coriol,
                0303      I       ikey,
                0304      O       hbl, bfsfc, stable, casea, kbl, Rib, sigma,
                0305      I       bi, bj, myTime, myIter, myThid )
2fa42a6013 Alis*0306
                0307 c     the oceanic planetary boundary layer depth, hbl, is determined as
                0308 c     the shallowest depth where the bulk Richardson number is
                0309 c     equal to the critical value, Ricr.
                0310 c
                0311 c     bulk Richardson numbers are evaluated by computing velocity and
                0312 c     buoyancy differences between values at zgrid(kl) < 0 and surface
                0313 c     reference values.
                0314 c     in this configuration, the reference values are equal to the
                0315 c     values in the surface layer.
                0316 c     when using a very fine vertical grid, these values should be
                0317 c     computed as the vertical average of velocity and buoyancy from
                0318 c     the surface down to epsilon*zgrid(kl).
                0319 c
                0320 c     when the bulk Richardson number at k exceeds Ricr, hbl is
                0321 c     linearly interpolated between grid levels zgrid(k) and zgrid(k-1).
                0322 c
                0323 c     The water column and the surface forcing are diagnosed for
                0324 c     stable/ustable forcing conditions, and where hbl is relative
                0325 c     to grid points (caseA), so that conditional branches can be
                0326 c     avoided in later subroutines.
a9d2e4c565 Jean*0327 c
2fa42a6013 Alis*0328       IMPLICIT NONE
                0329
                0330 #include "SIZE.h"
30c6f5b1cd An T*0331 #include "EEPARAMS.h"
                0332 #include "PARAMS.h"
2fa42a6013 Alis*0333 #include "KPP_PARAMS.h"
853c5e0e2c Jean*0334 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*0335 # include "tamc.h"
                0336 #endif
2fa42a6013 Alis*0337
                0338 c input
                0339 c------
15fd724cec Dimi*0340 c   bi, bj :: Array indices on which to apply calculations
eeda3a3aa1 Jean*0341 c   myTime :: Current time in simulation
                0342 c   myIter :: Current iteration number in simulation
                0343 c   myThid :: My Thread Id. number
2fa42a6013 Alis*0344 c kmtj      : number of vertical layers
                0345 c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
                0346 c dbloc     : local delta buoyancy across interfaces  (m/s^2)
                0347 c Ritop     : numerator of bulk Richardson Number
                0348 c             =(z-zref)*dbsfc, where dbsfc=delta
                0349 c             buoyancy with respect to surface      ((m/s)^2)
                0350 c ustar     : surface friction velocity                 (m/s)
                0351 c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
                0352 c bosol     : radiative buoyancy forcing            (m^2/s^3)
63ceaaa79c Dimi*0353 c boplume   : haline buoyancy forcing               (m^2/s^3)
2fa42a6013 Alis*0354 c coriol    : Coriolis parameter                        (1/s)
15fd724cec Dimi*0355       INTEGER bi, bj
eeda3a3aa1 Jean*0356       _RL     myTime
2b4c90c108 Mart*0357       INTEGER myIter
                0358       INTEGER myThid
                0359       INTEGER kmtj(imt)
fe66051ebd Dimi*0360       _RL dvsq    (imt,Nr)
                0361       _RL dbloc   (imt,Nr)
                0362       _RL Ritop   (imt,Nr)
                0363       _RL ustar   (imt)
                0364       _RL bo      (imt)
                0365       _RL bosol   (imt)
                0366       _RL coriol  (imt)
edb6656069 Mart*0367       INTEGER ikey
63ceaaa79c Dimi*0368 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0369       _RL boplume (imt,Nrp1)
63ceaaa79c Dimi*0370       _RL SPDepth (imt)
1f89baba18 Patr*0371 #ifdef SALT_PLUME_SPLIT_BASIN
                0372       _RL lon (imt)
                0373       _RL lat (imt)
                0374 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0375 #endif /* ALLOW_SALT_PLUME */
2fa42a6013 Alis*0376
                0377 c  output
                0378 c--------
                0379 c hbl       : boundary layer depth                        (m)
                0380 c bfsfc     : Bo+radiation absorbed to d=hbf*hbl    (m^2/s^3)
                0381 c stable    : =1 in stable forcing; =0 unstable
                0382 c casea     : =1 in case A, =0 in case B
                0383 c kbl       : -1 of first grid level below hbl
                0384 c Rib       : Bulk Richardson number
                0385 c sigma     : normalized depth (d/hbl)
fe66051ebd Dimi*0386       _RL hbl    (imt)
                0387       _RL bfsfc  (imt)
                0388       _RL stable (imt)
                0389       _RL casea  (imt)
2b4c90c108 Mart*0390       INTEGER kbl(imt)
fe66051ebd Dimi*0391       _RL Rib    (imt,Nr)
                0392       _RL sigma  (imt)
2fa42a6013 Alis*0393
                0394 #ifdef ALLOW_KPP
                0395
                0396 c  local
                0397 c-------
                0398 c wm, ws    : turbulent velocity scales         (m/s)
fe66051ebd Dimi*0399       _RL wm(imt), ws(imt)
                0400       _RL worka(imt)
                0401       _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
2b4c90c108 Mart*0402       INTEGER i, kl
2fa42a6013 Alis*0403
fe66051ebd Dimi*0404       _RL         p5    , eins
2b4c90c108 Mart*0405       PARAMETER ( p5=0.5, eins=1.0 )
956c0a5824 Patr*0406       _RL         minusone
2b4c90c108 Mart*0407       PARAMETER ( minusone=-1.0 )
29e16c9d38 Jean*0408 #ifdef SALT_PLUME_VOLUME
2b4c90c108 Mart*0409       INTEGER km, km1
29e16c9d38 Jean*0410       _RL temp
                0411 #endif
88b144ae49 Jean*0412 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0413       INTEGER kkey
88b144ae49 Jean*0414 #endif
2fa42a6013 Alis*0415
30c6f5b1cd An T*0416 #ifdef ALLOW_DIAGNOSTICS
                0417 c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
                0418       _RL KPPBFSFC(imt,Nr)
                0419 #endif /* ALLOW_DIAGNOSTICS */
                0420
2fa42a6013 Alis*0421 c find bulk Richardson number at every grid level until > Ricr
                0422 c
                0423 c note: the reference depth is -epsilon/2.*zgrid(k), but the reference
                0424 c       u,v,t,s values are simply the surface layer values,
                0425 c       and not the averaged values from 0 to 2*ref.depth,
                0426 c       which is necessary for very fine grids(top layer < 2m thickness)
                0427 c note: max values when Ricr never satisfied are
                0428 c       kbl(i)=kmtj(i) and hbl(i)=-zgrid(kmtj(i))
                0429
                0430 c     initialize hbl and kbl to bottomed out values
                0431
2b4c90c108 Mart*0432       DO i = 1, imt
b5ecce171d Davi*0433          Rib(i,1) = 0. _d 0
80b2748a09 Patr*0434          kbl(i) = kmtj(i)
9be8f21400 Mart*0435          IF (kmtj(i).LT.1) kbl(i) = 1
80b2748a09 Patr*0436          kl     = kbl(i)
                0437          hbl(i) = -zgrid(kl)
2b4c90c108 Mart*0438       ENDDO
2fa42a6013 Alis*0439
30c6f5b1cd An T*0440 #ifdef ALLOW_DIAGNOSTICS
2b4c90c108 Mart*0441       DO kl = 1, Nr
                0442          DO i = 1, imt
b5ecce171d Davi*0443             KPPBFSFC(i,kl) = 0. _d 0
2b4c90c108 Mart*0444          ENDDO
                0445       ENDDO
30c6f5b1cd An T*0446 #endif /* ALLOW_DIAGNOSTICS */
                0447
2b4c90c108 Mart*0448       DO kl = 2, Nr
2fa42a6013 Alis*0449
a9d2e4c565 Jean*0450 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0451          kkey = (ikey-1)*Nr + kl
3067745c9b Patr*0452 #endif
                0453
2fa42a6013 Alis*0454 c     compute bfsfc = sw fraction at hbf * zgrid
                0455
2b4c90c108 Mart*0456          DO i = 1, imt
956c0a5824 Patr*0457             worka(i) = zgrid(kl)
2b4c90c108 Mart*0458          ENDDO
853c5e0e2c Jean*0459 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0460 CADJ store worka = comlev1_kpp_k, key = kkey, kind=isbyte
853c5e0e2c Jean*0461 #endif
2b4c90c108 Mart*0462          CALL SWFRAC(
a9d2e4c565 Jean*0463      I       imt, hbf,
eeda3a3aa1 Jean*0464      U       worka,
                0465      I       myTime, myIter, myThid )
853c5e0e2c Jean*0466 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0467 CADJ store worka = comlev1_kpp_k, key = kkey, kind=isbyte
853c5e0e2c Jean*0468 #endif
3067745c9b Patr*0469
2b4c90c108 Mart*0470          DO i = 1, imt
2fa42a6013 Alis*0471
                0472 c     use caseA as temporary array
                0473
                0474             casea(i) = -zgrid(kl)
                0475
                0476 c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
                0477
956c0a5824 Patr*0478             bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
2fa42a6013 Alis*0479
2b4c90c108 Mart*0480          ENDDO
63ceaaa79c Dimi*0481 #ifdef ALLOW_SALT_PLUME
                0482 c     compute bfsfc = plume fraction at hbf * zgrid
30c6f5b1cd An T*0483          IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0484 #ifndef SALT_PLUME_VOLUME
2b4c90c108 Mart*0485            DO i = 1, imt
30c6f5b1cd An T*0486               worka(i) = zgrid(kl)
2b4c90c108 Mart*0487            ENDDO
1f89baba18 Patr*0488 Ccatn: in original way: accumulate all fractions of boplume above zgrid(kl)
2b4c90c108 Mart*0489            CALL SALT_PLUME_FRAC(
30c6f5b1cd An T*0490      I         imt, hbf,SPDepth,
1f89baba18 Patr*0491 #ifdef SALT_PLUME_SPLIT_BASIN
                0492      I         lon,lat,
                0493 #endif /* SALT_PLUME_SPLIT_BASIN */
30c6f5b1cd An T*0494      U         worka,
                0495      I         myTime, myIter, myThid)
2b4c90c108 Mart*0496            DO i = 1, imt
1f89baba18 Patr*0497               bfsfc(i) = bfsfc(i) + boplume(i,1)*(worka(i))
2b4c90c108 Mart*0498 C            km=MAX(1,kbl(i)-1)
1f89baba18 Patr*0499 C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
                0500 C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
2b4c90c108 Mart*0501            ENDDO
1f89baba18 Patr*0502 #else /* def SALT_PLUME_VOLUME */
                0503 catn: in vol way: need to integrate down to hbl, so first locate
                0504 c     k level associated with this hbl, then sum up all SPforc[T,S]
                0505            DO i = 1, imt
2b4c90c108 Mart*0506             km =MAX(1,kbl(i)-1)
                0507             km1=MAX(1,kbl(i))
29e16c9d38 Jean*0508             temp = (boplume(i,km)+boplume(i,km1))*p5
1f89baba18 Patr*0509             bfsfc(i) = bfsfc(i) + temp
                0510            ENDDO
                0511 #endif /* ndef SALT_PLUME_VOLUME */
30c6f5b1cd An T*0512          ENDIF
                0513 #endif /* ALLOW_SALT_PLUME */
                0514
                0515 #ifdef ALLOW_DIAGNOSTICS
2b4c90c108 Mart*0516          DO i = 1, imt
30c6f5b1cd An T*0517             KPPBFSFC(i,kl) = bfsfc(i)
2b4c90c108 Mart*0518          ENDDO
30c6f5b1cd An T*0519 #endif /* ALLOW_DIAGNOSTICS */
63ceaaa79c Dimi*0520
2b4c90c108 Mart*0521          DO i = 1, imt
63ceaaa79c Dimi*0522             stable(i) = p5 + sign(p5,bfsfc(i))
                0523             sigma(i) = stable(i) + (1. - stable(i)) * epsilon
2b4c90c108 Mart*0524          ENDDO
2fa42a6013 Alis*0525
                0526 c-----------------------------------------------------------------------
                0527 c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
                0528 c-----------------------------------------------------------------------
a9d2e4c565 Jean*0529
2b4c90c108 Mart*0530          CALL wscale (
2fa42a6013 Alis*0531      I        sigma, casea, ustar, bfsfc,
25c3477f99 Mart*0532      O        wm, ws, myThid )
853c5e0e2c Jean*0533 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0534 CADJ store ws = comlev1_kpp_k, key = kkey, kind=isbyte
853c5e0e2c Jean*0535 #endif
2fa42a6013 Alis*0536
2b4c90c108 Mart*0537          DO i = 1, imt
2fa42a6013 Alis*0538
                0539 c-----------------------------------------------------------------------
                0540 c     compute the turbulent shear contribution to Rib
                0541 c-----------------------------------------------------------------------
                0542
                0543             bvsq = p5 *
                0544      1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
                0545      2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
2b4c90c108 Mart*0546 CMLC     Van Roekel et al 2018 suggest this, but our solution is
                0547 CMLC     probably OK, too:
                0548 CML            bvsq = MAX
                0549 CML     1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  )),
                0550 CML     2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
2fa42a6013 Alis*0551
2b4c90c108 Mart*0552             IF (bvsq .EQ. 0. _d 0) THEN
b5ecce171d Davi*0553               vtsq = 0. _d 0
2b4c90c108 Mart*0554             ELSE
                0555               vtsq = -zgrid(kl) * ws(i) * SQRT(ABS(bvsq)) * Vtc
                0556             ENDIF
2fa42a6013 Alis*0557
                0558 c     compute bulk Richardson number at new level
                0559 c     note: Ritop needs to be zero on land and ocean bottom
                0560 c     points so that the following if statement gets triggered
                0561 c     correctly; otherwise, hbl might get set to (big) negative
                0562 c     values, that might exceed the limit for the "exp" function
                0563 c     in "SWFRAC"
                0564
                0565 c
                0566 c     rg: assignment to double precision variable to avoid overflow
                0567 c     ph: test for zero nominator
                0568 c
956c0a5824 Patr*0569
a9d2e4c565 Jean*0570             tempVar1  = dvsq(i,kl) + vtsq
2b4c90c108 Mart*0571 #ifdef KPP_SMOOTH_REGULARISATION
                0572             tempVar2 = tempVar1 + phepsi
                0573 #else
                0574             tempVar2 = MAX(tempVar1, phepsi)
                0575 #endif /* KPP_SMOOTH_REGULARISATION */
956c0a5824 Patr*0576             Rib(i,kl) = Ritop(i,kl) / tempVar2
                0577
2b4c90c108 Mart*0578          ENDDO
                0579       ENDDO
32e2940ee8 Patr*0580
30c6f5b1cd An T*0581 #ifdef ALLOW_DIAGNOSTICS
52b77e9710 Jean*0582       IF ( useDiagnostics ) THEN
15fd724cec Dimi*0583          CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,2,bi,bj,myThid)
2b4c90c108 Mart*0584          CALL DIAGNOSTICS_FILL(Rib     ,'KPPRi   ',0,Nr,2,bi,bj,myThid)
52b77e9710 Jean*0585       ENDIF
30c6f5b1cd An T*0586 #endif /* ALLOW_DIAGNOSTICS */
                0587
32e2940ee8 Patr*0588 cph(
04522666de Ed H*0589 cph  without this store, there is a recomputation error for
a9d2e4c565 Jean*0590 cph  rib in adbldepth (probably partial recomputation problem)
853c5e0e2c Jean*0591 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0592 CADJ store Rib = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0593 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
853c5e0e2c Jean*0594 #endif
32e2940ee8 Patr*0595 cph)
                0596
2b4c90c108 Mart*0597       DO kl = 2, Nr
                0598          DO i = 1, imt
                0599             IF (kbl(i).EQ.kmtj(i) .AND. Rib(i,kl).GT.Ricr) kbl(i) = kl
                0600          ENDDO
                0601       ENDDO
2fa42a6013 Alis*0602
853c5e0e2c Jean*0603 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0604 CADJ store kbl = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0605 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0606 #endif
2fa42a6013 Alis*0607
2b4c90c108 Mart*0608       DO i = 1, imt
2fa42a6013 Alis*0609          kl = kbl(i)
                0610 c     linearly interpolate to find hbl where Rib = Ricr
2b4c90c108 Mart*0611          IF (kl.GT.1 .AND. kl.LT.kmtj(i)) THEN
956c0a5824 Patr*0612             tempVar1 = (Rib(i,kl)-Rib(i,kl-1))
2fa42a6013 Alis*0613             hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
956c0a5824 Patr*0614      1           (Ricr - Rib(i,kl-1)) / tempVar1
2b4c90c108 Mart*0615 C     this is the MOM5 formulation, (nearly) identical results
                0616 CML            hbl(i) = ((Rib(i,kl-1)-Ricr)*zgrid(kl) -
                0617 CML     1           (Rib(i,kl)-Ricr)*zgrid(kl-1))/tempVar1
                0618          ENDIF
                0619       ENDDO
2fa42a6013 Alis*0620
853c5e0e2c Jean*0621 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0622 CADJ store hbl = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0623 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0624 #endif
2fa42a6013 Alis*0625
                0626 c-----------------------------------------------------------------------
                0627 c     find stability and buoyancy forcing for boundary layer
                0628 c-----------------------------------------------------------------------
                0629
2b4c90c108 Mart*0630       DO i = 1, imt
956c0a5824 Patr*0631          worka(i) = hbl(i)
2b4c90c108 Mart*0632       ENDDO
853c5e0e2c Jean*0633 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0634 CADJ store worka = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0635 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0636 #endif
2b4c90c108 Mart*0637       CALL SWFRAC(
a9d2e4c565 Jean*0638      I       imt, minusone,
eeda3a3aa1 Jean*0639      U       worka,
                0640      I       myTime, myIter, myThid )
853c5e0e2c Jean*0641 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0642 CADJ store worka = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0643 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0644 #endif
2fa42a6013 Alis*0645
2b4c90c108 Mart*0646       DO i = 1, imt
956c0a5824 Patr*0647          bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
2b4c90c108 Mart*0648       ENDDO
63ceaaa79c Dimi*0649
                0650 #ifdef ALLOW_SALT_PLUME
7448700841 Mart*0651 # ifdef ALLOW_AUTODIFF_TAMC
                0652 C     suppress a TAF warning
                0653 CADJ incomplete worka
                0654 # endif
30c6f5b1cd An T*0655       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0656 #ifndef SALT_PLUME_VOLUME
2b4c90c108 Mart*0657         DO i = 1, imt
30c6f5b1cd An T*0658            worka(i) = hbl(i)
2b4c90c108 Mart*0659         ENDDO
                0660         CALL SALT_PLUME_FRAC(
30c6f5b1cd An T*0661      I         imt,minusone,SPDepth,
1f89baba18 Patr*0662 #ifdef SALT_PLUME_SPLIT_BASIN
                0663      I         lon,lat,
                0664 #endif /* SALT_PLUME_SPLIT_BASIN */
30c6f5b1cd An T*0665      U         worka,
                0666      I         myTime, myIter, myThid )
2b4c90c108 Mart*0667         DO i = 1, imt
1f89baba18 Patr*0668            bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
2b4c90c108 Mart*0669 C            km=MAX(1,kbl(i)-1)
1f89baba18 Patr*0670 C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
                0671 C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
2b4c90c108 Mart*0672         ENDDO
1f89baba18 Patr*0673 #else /* def SALT_PLUME_VOLUME */
                0674         DO i = 1, imt
2b4c90c108 Mart*0675             km =MAX(1,kbl(i)-1)
                0676             km1=MAX(1,kbl(i))
1f89baba18 Patr*0677             temp = (boplume(i,km)+boplume(i,km1))/2.0
                0678             bfsfc(i) = bfsfc(i) + temp
                0679         ENDDO
                0680 #endif /* ndef SALT_PLUME_VOLUME */
30c6f5b1cd An T*0681       ENDIF
63ceaaa79c Dimi*0682 #endif /* ALLOW_SALT_PLUME */
853c5e0e2c Jean*0683 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0684 CADJ store bfsfc = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0685 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0686 #endif
1d478690dc Patr*0687
956c0a5824 Patr*0688 c--   ensure bfsfc is never 0
2b4c90c108 Mart*0689       DO i = 1, imt
2fa42a6013 Alis*0690          stable(i) = p5 + sign( p5, bfsfc(i) )
2b4c90c108 Mart*0691          bfsfc(i) = sign(eins,bfsfc(i))*MAX(phepsi,ABS(bfsfc(i)))
                0692       ENDDO
2fa42a6013 Alis*0693
32e2940ee8 Patr*0694 cph(
                0695 cph  added stable to store list to avoid extensive recomp.
853c5e0e2c Jean*0696 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0697 CADJ store bfsfc, stable = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0698 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0699 #endif
32e2940ee8 Patr*0700 cph)
956c0a5824 Patr*0701
2fa42a6013 Alis*0702 c-----------------------------------------------------------------------
                0703 c check hbl limits for hekman or hmonob
                0704 c     ph: test for zero nominator
                0705 c-----------------------------------------------------------------------
                0706
35936beffd Davi*0707       IF ( LimitHblStable ) THEN
2b4c90c108 Mart*0708        DO i = 1, imt
                0709         IF (bfsfc(i) .GT. 0.0) THEN
                0710          hekman = cekman * ustar(i) / MAX(ABS(Coriol(i)),phepsi)
                0711          hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
                0712      &        / vonk / bfsfc(i)
                0713          hlimit = stable(i) * MIN(hekman,hmonob)
                0714      &        + (stable(i)-1.) * zgrid(Nr)
                0715          hbl(i) = MIN(hbl(i),hlimit)
                0716         ENDIF
                0717        ENDDO
35936beffd Davi*0718       ENDIF
                0719
853c5e0e2c Jean*0720 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0721 CADJ store hbl = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0722 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0723 #endif
956c0a5824 Patr*0724
2b4c90c108 Mart*0725       DO i = 1, imt
                0726          hbl(i) = MAX(hbl(i),minKPPhbl)
2fa42a6013 Alis*0727          kbl(i) = kmtj(i)
2b4c90c108 Mart*0728       ENDDO
2fa42a6013 Alis*0729
853c5e0e2c Jean*0730 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0731 CADJ store hbl = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0732 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0733 #endif
2fa42a6013 Alis*0734
                0735 c-----------------------------------------------------------------------
                0736 c      find new kbl
                0737 c-----------------------------------------------------------------------
                0738
2b4c90c108 Mart*0739       DO kl = 2, Nr
                0740          DO i = 1, imt
                0741             IF ( kbl(i).EQ.kmtj(i) .AND. (-zgrid(kl)).GT.hbl(i) ) THEN
2fa42a6013 Alis*0742                kbl(i) = kl
2b4c90c108 Mart*0743             ENDIF
                0744          ENDDO
                0745       ENDDO
2fa42a6013 Alis*0746
                0747 c-----------------------------------------------------------------------
                0748 c      find stability and buoyancy forcing for final hbl values
                0749 c-----------------------------------------------------------------------
                0750
2b4c90c108 Mart*0751       DO i = 1, imt
956c0a5824 Patr*0752          worka(i) = hbl(i)
2b4c90c108 Mart*0753       ENDDO
853c5e0e2c Jean*0754 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0755 CADJ store worka = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0756 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0757 #endif
2b4c90c108 Mart*0758       CALL SWFRAC(
a9d2e4c565 Jean*0759      I     imt, minusone,
eeda3a3aa1 Jean*0760      U     worka,
                0761      I     myTime, myIter, myThid )
853c5e0e2c Jean*0762 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0763 CADJ store worka = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0764 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0765 #endif
a9d2e4c565 Jean*0766
2b4c90c108 Mart*0767       DO i = 1, imt
956c0a5824 Patr*0768          bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
2b4c90c108 Mart*0769       ENDDO
63ceaaa79c Dimi*0770
                0771 #ifdef ALLOW_SALT_PLUME
30c6f5b1cd An T*0772       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0773 #ifndef SALT_PLUME_VOLUME
2b4c90c108 Mart*0774         DO i = 1, imt
30c6f5b1cd An T*0775            worka(i) = hbl(i)
2b4c90c108 Mart*0776         ENDDO
                0777         CALL SALT_PLUME_FRAC(
30c6f5b1cd An T*0778      I         imt,minusone,SPDepth,
1f89baba18 Patr*0779 #ifdef SALT_PLUME_SPLIT_BASIN
                0780      I         lon,lat,
                0781 #endif /* SALT_PLUME_SPLIT_BASIN */
30c6f5b1cd An T*0782      U         worka,
                0783      I         myTime, myIter, myThid )
2b4c90c108 Mart*0784         DO i = 1, imt
1f89baba18 Patr*0785            bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
2b4c90c108 Mart*0786 C            km=MAX(1,kbl(i)-1)
1f89baba18 Patr*0787 C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
                0788 C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
2b4c90c108 Mart*0789         ENDDO
1f89baba18 Patr*0790 #else /* def SALT_PLUME_VOLUME */
                0791         DO i = 1, imt
2b4c90c108 Mart*0792             km =MAX(1,kbl(i)-1)
                0793             km1=MAX(1,kbl(i)-0)
1f89baba18 Patr*0794             temp = (boplume(i,km)+boplume(i,km1))/2.0
                0795             bfsfc(i) = bfsfc(i) + temp
                0796         ENDDO
                0797 #endif /* ndef SALT_PLUME_VOLUME */
30c6f5b1cd An T*0798       ENDIF
63ceaaa79c Dimi*0799 #endif /* ALLOW_SALT_PLUME */
853c5e0e2c Jean*0800 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0801 CADJ store bfsfc = comlev1_kpp, key=ikey, kind=isbyte,
95c72ef3a1 Patr*0802 CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
853c5e0e2c Jean*0803 #endif
1d478690dc Patr*0804
956c0a5824 Patr*0805 c--   ensures bfsfc is never 0
2b4c90c108 Mart*0806       DO i = 1, imt
2fa42a6013 Alis*0807          stable(i) = p5 + sign( p5, bfsfc(i) )
2b4c90c108 Mart*0808          bfsfc(i) = sign(eins,bfsfc(i))*MAX(phepsi,ABS(bfsfc(i)))
                0809       ENDDO
2fa42a6013 Alis*0810
                0811 c-----------------------------------------------------------------------
                0812 c determine caseA and caseB
                0813 c-----------------------------------------------------------------------
                0814
2b4c90c108 Mart*0815       DO i = 1, imt
80b2748a09 Patr*0816          kl = kbl(i)
2fa42a6013 Alis*0817          casea(i) = p5 +
80b2748a09 Patr*0818      1        sign(p5, -zgrid(kl) - p5*hwide(kl) - hbl(i))
2b4c90c108 Mart*0819       ENDDO
2fa42a6013 Alis*0820
                0821 #endif /* ALLOW_KPP */
                0822
2b4c90c108 Mart*0823       RETURN
                0824       END
2fa42a6013 Alis*0825
                0826 c*************************************************************************
                0827
                0828       subroutine wscale (
                0829      I     sigma, hbl, ustar, bfsfc,
a9d2e4c565 Jean*0830      O     wm, ws,
25c3477f99 Mart*0831      I     myThid )
2fa42a6013 Alis*0832
                0833 c     compute turbulent velocity scales.
                0834 c     use a 2D-lookup table for wm and ws as functions of ustar and
                0835 c     zetahat (=vonk*sigma*hbl*bfsfc).
                0836 c
                0837 c     note: the lookup table is only used for unstable conditions
2b4c90c108 Mart*0838 c     (zehat.LE.0), in the stable domain wm (=ws) gets computed
2fa42a6013 Alis*0839 c     directly.
a9d2e4c565 Jean*0840 c
2fa42a6013 Alis*0841       IMPLICIT NONE
                0842
                0843 #include "SIZE.h"
                0844 #include "KPP_PARAMS.h"
                0845
                0846 c input
                0847 c------
                0848 c sigma   : normalized depth (d/hbl)
                0849 c hbl     : boundary layer depth (m)
                0850 c ustar   : surface friction velocity         (m/s)
                0851 c bfsfc   : total surface buoyancy flux       (m^2/s^3)
25c3477f99 Mart*0852 c myThid  : thread number for this instance of the routine
2b4c90c108 Mart*0853       INTEGER myThid
fe66051ebd Dimi*0854       _RL sigma(imt)
                0855       _RL hbl  (imt)
                0856       _RL ustar(imt)
                0857       _RL bfsfc(imt)
a9d2e4c565 Jean*0858
2fa42a6013 Alis*0859 c  output
                0860 c--------
                0861 c wm, ws  : turbulent velocity scales at sigma
fe66051ebd Dimi*0862       _RL wm(imt), ws(imt)
2fa42a6013 Alis*0863
                0864 #ifdef ALLOW_KPP
a9d2e4c565 Jean*0865
2fa42a6013 Alis*0866 c local
                0867 c------
                0868 c zehat   : = zeta *  ustar**3
fe66051ebd Dimi*0869       _RL zehat
2fa42a6013 Alis*0870
2b4c90c108 Mart*0871       INTEGER iz, izp1, ju, i, jup1
fe66051ebd Dimi*0872       _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
                0873       _RL wbm, was, wbs, u3, tempVar
2fa42a6013 Alis*0874
                0875 c-----------------------------------------------------------------------
                0876 c use lookup table for zehat < zmax only; otherwise use
                0877 c stable formulae
                0878 c-----------------------------------------------------------------------
                0879
2b4c90c108 Mart*0880       DO i = 1, imt
2fa42a6013 Alis*0881          zehat = vonk*sigma(i)*hbl(i)*bfsfc(i)
                0882
2b4c90c108 Mart*0883          IF (zehat .LE. zmax) THEN
2fa42a6013 Alis*0884
                0885             zdiff = zehat - zmin
36d05ba58c Mart*0886 C     For extremely negative buoyancy forcing bfsfc, zehat and hence
                0887 C     zdiff can become very negative (default value of zmin = 4.e-7) and
                0888 C     the extrapolation beyond the limit zmin of the lookup table can
                0889 C     give very bad values and may make the model crash. Here is a
                0890 C     simple fix (thanks to Dimitry Sidorenko) that effectively replaces
                0891 C     linear extrapolation with nearest neighbor extrapolation so that
                0892 C     only the lower limit values of the lookup tables wmt/wst are used.
                0893 C     Alternatively, one can get rid of the lookup table altogether
                0894 C     and compute the coefficients online (done in NEMO, for example).
2b4c90c108 Mart*0895 C           zdiff = MAX( 0. _d 0, zehat - zmin )
                0896             iz    = INT( zdiff / deltaz )
                0897             iz    = MIN( iz, nni )
                0898             iz    = MAX( iz, 0 )
2fa42a6013 Alis*0899             izp1  = iz + 1
956c0a5824 Patr*0900
2fa42a6013 Alis*0901             udiff = ustar(i) - umin
2b4c90c108 Mart*0902             ju    = INT( udiff / deltau )
                0903             ju    = MIN( ju, nnj )
                0904             ju    = MAX( ju, 0 )
2fa42a6013 Alis*0905             jup1  = ju + 1
956c0a5824 Patr*0906
2fa42a6013 Alis*0907             zfrac = zdiff / deltaz - float(iz)
                0908             ufrac = udiff / deltau - float(ju)
956c0a5824 Patr*0909
2fa42a6013 Alis*0910             fzfrac= 1. - zfrac
                0911             wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
                0912             wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )
                0913             wm(i) = (1.-ufrac) * wbm          + ufrac * wam
956c0a5824 Patr*0914
2fa42a6013 Alis*0915             was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)
                0916             wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )
                0917             ws(i) = (1.-ufrac) * wbs          + ufrac * was
956c0a5824 Patr*0918
2b4c90c108 Mart*0919          ELSE
956c0a5824 Patr*0920
                0921             u3 = ustar(i) * ustar(i) * ustar(i)
                0922             tempVar = u3 + conc1 * zehat
                0923             wm(i) = vonk * ustar(i) * u3 / tempVar
2fa42a6013 Alis*0924             ws(i) = wm(i)
956c0a5824 Patr*0925
2b4c90c108 Mart*0926          ENDIF
a9d2e4c565 Jean*0927
2b4c90c108 Mart*0928       ENDDO
2fa42a6013 Alis*0929
                0930 #endif /* ALLOW_KPP */
a9d2e4c565 Jean*0931
2b4c90c108 Mart*0932       RETURN
                0933       END
2fa42a6013 Alis*0934
                0935 c*************************************************************************
a9d2e4c565 Jean*0936
2fa42a6013 Alis*0937       subroutine Ri_iwmix (
eeda3a3aa1 Jean*0938      I       kmtj, shsq, dbloc, dblocSm,
                0939      I       diffusKzS, diffusKzT,
edb6656069 Mart*0940      I       ikey,
eeda3a3aa1 Jean*0941      O       diffus,
                0942      I       myThid )
2fa42a6013 Alis*0943
                0944 c     compute interior viscosity diffusivity coefficients due
                0945 c     to shear instability (dependent on a local Richardson number),
                0946 c     to background internal wave activity, and
                0947 c     to static instability (local Richardson number < 0).
                0948
                0949       IMPLICIT NONE
                0950
                0951 #include "SIZE.h"
                0952 #include "EEPARAMS.h"
                0953 #include "PARAMS.h"
                0954 #include "KPP_PARAMS.h"
95c72ef3a1 Patr*0955 #ifdef ALLOW_AUTODIFF
5127d1d91b Jean*0956 # include "AUTODIFF_PARAMS.h"
853c5e0e2c Jean*0957 #endif
                0958 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*0959 # include "tamc.h"
                0960 #endif
2fa42a6013 Alis*0961
                0962 c  input
                0963 c     kmtj   (imt)         number of vertical layers on this row
                0964 c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
                0965 c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
                0966 c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
059d9fc14f Dimi*0967 c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
                0968 c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
eeda3a3aa1 Jean*0969 c     myThid :: My Thread Id. number
2b4c90c108 Mart*0970       INTEGER kmtj (imt)
fe66051ebd Dimi*0971       _RL shsq     (imt,Nr)
                0972       _RL dbloc    (imt,Nr)
                0973       _RL dblocSm  (imt,Nr)
059d9fc14f Dimi*0974       _RL diffusKzS(imt,Nr)
                0975       _RL diffusKzT(imt,Nr)
edb6656069 Mart*0976       INTEGER ikey
2b4c90c108 Mart*0977       INTEGER myThid
a9d2e4c565 Jean*0978
2fa42a6013 Alis*0979 c output
                0980 c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
                0981 c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
                0982 c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
fe66051ebd Dimi*0983       _RL diffus(imt,0:Nrp1,3)
2fa42a6013 Alis*0984
                0985 #ifdef ALLOW_KPP
                0986
                0987 c local variables
                0988 c     Rig                   local Richardson number
                0989 c     fRi, fcon             function of Rig
fe66051ebd Dimi*0990       _RL Rig
                0991       _RL fRi, fcon
                0992       _RL ratio
2b4c90c108 Mart*0993       INTEGER i, ki, kp1
fe66051ebd Dimi*0994       _RL c1, c0
2fa42a6013 Alis*0995
8164aa1823 Patr*0996 #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
2b4c90c108 Mart*0997       INTEGER mr
8164aa1823 Patr*0998 CADJ INIT kpp_ri_tape_mr = common, 1
                0999 #endif
                1000
2fa42a6013 Alis*1001 c constants
b5ecce171d Davi*1002       c1 = 1. _d 0
                1003       c0 = 0. _d 0
2fa42a6013 Alis*1004
                1005 c-----------------------------------------------------------------------
                1006 c     compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)
                1007 c     use diffus(*,*,1) as temporary storage of Ri to be smoothed
                1008 c     use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared
                1009 c     set values at bottom and below to nearest value above bottom
853c5e0e2c Jean*1010 #ifdef ALLOW_AUTODIFF
8164aa1823 Patr*1011 C     break data flow dependence on diffus
                1012       diffus(1,1,1) = 0.0
2b4c90c108 Mart*1013       DO ki = 1, Nr
                1014          DO i = 1, imt
3067745c9b Patr*1015             diffus(i,ki,1) = 0.
                1016             diffus(i,ki,2) = 0.
                1017             diffus(i,ki,3) = 0.
2b4c90c108 Mart*1018          ENDDO
                1019       ENDDO
8164aa1823 Patr*1020 #endif
a9d2e4c565 Jean*1021
2b4c90c108 Mart*1022       DO ki = 1, Nr
                1023          DO i = 1, imt
                1024             IF     (kmtj(i) .LE. 1      ) THEN
2fa42a6013 Alis*1025                diffus(i,ki,1) = 0.
                1026                diffus(i,ki,2) = 0.
2b4c90c108 Mart*1027             ELSEIF (ki      .GE. kmtj(i)) THEN
2fa42a6013 Alis*1028                diffus(i,ki,1) = diffus(i,ki-1,1)
                1029                diffus(i,ki,2) = diffus(i,ki-1,2)
2b4c90c108 Mart*1030             ELSE
2fa42a6013 Alis*1031                diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))
2b4c90c108 Mart*1032 #ifdef KPP_SMOOTH_REGULARISATION
                1033      &            / ( Shsq(i,ki) + phepsi**2 )
                1034 #else
                1035      &            / MAX( Shsq(i,ki), phepsi )
                1036 #endif
2fa42a6013 Alis*1037                diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))
2b4c90c108 Mart*1038             ENDIF
                1039          ENDDO
                1040       ENDDO
853c5e0e2c Jean*1041 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1042 CADJ store diffus = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1043 #endif
2fa42a6013 Alis*1044
                1045 c-----------------------------------------------------------------------
                1046 c     vertically smooth Ri
8164aa1823 Patr*1047 #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
2b4c90c108 Mart*1048       DO mr = 1, num_v_smooth_Ri
1d478690dc Patr*1049
853c5e0e2c Jean*1050 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1051 CADJ store diffus(:,:,1) = kpp_ri_tape_mr, key=mr
                1052 CADJ &  , shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
                1053 CMLCADJ store diffus(:,:,2) = kpp_ri_tape_mr, key=mr
                1054 CMLCADJ &  , shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
853c5e0e2c Jean*1055 #endif
1d478690dc Patr*1056
2b4c90c108 Mart*1057          CALL z121 (
30ef4f2c65 Davi*1058      U     diffus(1,0,1),
25c3477f99 Mart*1059      I     myThid )
2b4c90c108 Mart*1060 C     it may also make sense to smooth buoyancy vertically
                1061 CML         CALL z121 (
                1062 CML     U     diffus(1,0,2),
                1063 CML     I     myThid )
                1064       ENDDO
8164aa1823 Patr*1065 #endif
2fa42a6013 Alis*1066
                1067 c-----------------------------------------------------------------------
                1068 c                           after smoothing loop
                1069
2b4c90c108 Mart*1070       DO ki = 1, Nr
                1071          DO i = 1, imt
a9d2e4c565 Jean*1072
2fa42a6013 Alis*1073 c  evaluate f of Brunt-Vaisala squared for convection, store in fcon
                1074
2b4c90c108 Mart*1075             Rig   = MAX ( diffus(i,ki,2) , BVSQcon )
                1076             ratio = MIN ( (BVSQcon - Rig) / BVSQcon, c1 )
2fa42a6013 Alis*1077             fcon  = c1 - ratio * ratio
                1078             fcon  = fcon * fcon * fcon
a9d2e4c565 Jean*1079
2fa42a6013 Alis*1080 c  evaluate f of smooth Ri for shear instability, store in fRi
                1081
2b4c90c108 Mart*1082             Rig  = MAX ( diffus(i,ki,1), c0 )
                1083             ratio = MIN ( Rig / Riinfty , c1 )
2fa42a6013 Alis*1084             fRi   = c1 - ratio * ratio
                1085             fRi   = fRi * fRi * fRi
2b4c90c108 Mart*1086 #ifdef KPP_SCALE_SHEARMIXING
                1087 C     reduce shear mixing when there is no shear (Polzin, 1996, JPO, 1409-1425)
                1088 C     importend from MOM5 code
                1089             fRi   = fRi * shsq(i,ki)*shsq(i,ki)
                1090      &           /(shsq(i,ki)*shsq(i,ki) + 1. _d -16)
                1091 #endif
2fa42a6013 Alis*1092 c ----------------------------------------------------------------------
                1093 c            evaluate diffusivities and viscosity
                1094 c    mixing due to internal waves, and shear and static instability
f7ab1e3823 Patr*1095
faff480f8c Davi*1096             kp1 = MIN(ki+1,Nr)
5127d1d91b Jean*1097 #ifdef EXCLUDE_KPP_SHEAR_MIX
a9d2e4c565 Jean*1098             diffus(i,ki,1) = viscArNr(1)
faff480f8c Davi*1099             diffus(i,ki,2) = diffusKzS(i,kp1)
                1100             diffus(i,ki,3) = diffusKzT(i,kp1)
5127d1d91b Jean*1101 #else /* EXCLUDE_KPP_SHEAR_MIX */
                1102 # ifdef ALLOW_AUTODIFF
2b4c90c108 Mart*1103             IF ( inAdMode .AND. .NOT. inAdExact ) THEN
5127d1d91b Jean*1104               diffus(i,ki,1) = viscArNr(1)
                1105               diffus(i,ki,2) = diffusKzS(i,kp1)
                1106               diffus(i,ki,3) = diffusKzT(i,kp1)
2b4c90c108 Mart*1107             ELSE
5127d1d91b Jean*1108 # else /* ALLOW_AUTODIFF */
2b4c90c108 Mart*1109             IF ( .TRUE. ) THEN
5127d1d91b Jean*1110 # endif /* ALLOW_AUTODIFF */
                1111               diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0
                1112               diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0
                1113               diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0
2b4c90c108 Mart*1114             ENDIF
5127d1d91b Jean*1115 #endif /* EXCLUDE_KPP_SHEAR_MIX */
2b4c90c108 Mart*1116          ENDDO
                1117       ENDDO
a9d2e4c565 Jean*1118
2fa42a6013 Alis*1119 c ------------------------------------------------------------------------
                1120 c         set surface values to 0.0
a9d2e4c565 Jean*1121
2b4c90c108 Mart*1122       DO i = 1, imt
2fa42a6013 Alis*1123          diffus(i,0,1) = c0
                1124          diffus(i,0,2) = c0
                1125          diffus(i,0,3) = c0
2b4c90c108 Mart*1126       ENDDO
2fa42a6013 Alis*1127
                1128 #endif /* ALLOW_KPP */
a9d2e4c565 Jean*1129
2b4c90c108 Mart*1130       RETURN
                1131       END
2fa42a6013 Alis*1132
                1133 c*************************************************************************
                1134
                1135       subroutine z121 (
25c3477f99 Mart*1136      U     v,
                1137      I     myThid )
2fa42a6013 Alis*1138
                1139 c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
                1140 c     top (0) value is used as a dummy
                1141 c     bottom (Nrp1) value is set to input value from above.
                1142
1d478690dc Patr*1143 c     Note that it is important to exclude from the smoothing any points
2fa42a6013 Alis*1144 c     that are outside the range of the K(Ri) scheme, ie.  >0.8, or <0.0.
1d478690dc Patr*1145 c     Otherwise, there is interference with other physics, especially
2fa42a6013 Alis*1146 c     penetrative convection.
                1147
                1148       IMPLICIT NONE
                1149 #include "SIZE.h"
                1150 #include "KPP_PARAMS.h"
                1151
a9d2e4c565 Jean*1152 c input/output
2fa42a6013 Alis*1153 c-------------
                1154 c v     : 2-D array to be smoothed in Nrp1 direction
25c3477f99 Mart*1155 c myThid: thread number for this instance of the routine
2b4c90c108 Mart*1156       INTEGER myThid
fe66051ebd Dimi*1157       _RL v(imt,0:Nrp1)
2fa42a6013 Alis*1158
                1159 #ifdef ALLOW_KPP
                1160
                1161 c local
fe66051ebd Dimi*1162       _RL zwork, zflag
                1163       _RL KRi_range(1:Nrp1)
2b4c90c108 Mart*1164       INTEGER i, k, km1, kp1
2fa42a6013 Alis*1165
fe66051ebd Dimi*1166       _RL         p0      , p25       , p5      , p2
2b4c90c108 Mart*1167       PARAMETER ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
1d478690dc Patr*1168
                1169       KRi_range(Nrp1) = p0
                1170
                1171 #ifdef ALLOW_AUTODIFF_TAMC
                1172 C--   dummy assignment to end declaration part for TAMC
                1173       i = 0
                1174 C--   HPF directive to help TAMC
                1175 CHPF$ INDEPENDENT
                1176 #endif /* ALLOW_AUTODIFF_TAMC */
956c0a5824 Patr*1177
2b4c90c108 Mart*1178       DO i = 1, imt
2fa42a6013 Alis*1179
956c0a5824 Patr*1180          k = 1
2fa42a6013 Alis*1181          v(i,Nrp1) = v(i,Nr)
                1182
2b4c90c108 Mart*1183          DO k = 1, Nr
2fa42a6013 Alis*1184             KRi_range(k) = p5 + SIGN(p5,v(i,k))
                1185             KRi_range(k) = KRi_range(k) *
                1186      &                     ( p5 + SIGN(p5,(Riinfty-v(i,k))) )
2b4c90c108 Mart*1187          ENDDO
2fa42a6013 Alis*1188
                1189          zwork  = KRi_range(1) * v(i,1)
                1190          v(i,1) = p2 * v(i,1) +
                1191      &            KRi_range(1) * KRi_range(2) * v(i,2)
                1192          zflag  = p2 + KRi_range(1) * KRi_range(2)
                1193          v(i,1) = v(i,1) / zflag
                1194
2b4c90c108 Mart*1195          DO k = 2, Nr
2fa42a6013 Alis*1196             km1 = k - 1
                1197             kp1 = k + 1
                1198             zflag = v(i,k)
                1199             v(i,k) = p2 * v(i,k) +
                1200      &           KRi_range(k) * KRi_range(kp1) * v(i,kp1) +
                1201      &           KRi_range(k) * zwork
                1202             zwork = KRi_range(k) * zflag
                1203             zflag = p2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1))
                1204             v(i,k) = v(i,k) / zflag
2b4c90c108 Mart*1205          ENDDO
2fa42a6013 Alis*1206
2b4c90c108 Mart*1207       ENDDO
2fa42a6013 Alis*1208
                1209 #endif /* ALLOW_KPP */
                1210
2b4c90c108 Mart*1211       RETURN
                1212       END
2fa42a6013 Alis*1213
                1214 c*************************************************************************
                1215
956c0a5824 Patr*1216       subroutine smooth_horiz (
2fa42a6013 Alis*1217      I     k, bi, bj,
25c3477f99 Mart*1218      U     fld,
                1219      I     myThid )
2fa42a6013 Alis*1220
956c0a5824 Patr*1221 c     Apply horizontal smoothing to global _RL 2-D array
2fa42a6013 Alis*1222
                1223       IMPLICIT NONE
                1224 #include "SIZE.h"
11b9c75340 Jean*1225 #include "GRID.h"
2fa42a6013 Alis*1226 #include "KPP_PARAMS.h"
                1227
                1228 c     input
956c0a5824 Patr*1229 c     bi, bj : array indices
                1230 c     k      : vertical index used for masking
25c3477f99 Mart*1231 c     myThid : thread number for this instance of the routine
                1232       INTEGER myThid
2b4c90c108 Mart*1233       INTEGER k, bi, bj
2fa42a6013 Alis*1234
956c0a5824 Patr*1235 c     input/output
                1236 c     fld    : 2-D array to be smoothed
                1237       _RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
2fa42a6013 Alis*1238
                1239 #ifdef ALLOW_KPP
                1240
                1241 c     local
2b4c90c108 Mart*1242       INTEGER i, j, im1, ip1, jm1, jp1
2fa42a6013 Alis*1243       _RL tempVar
956c0a5824 Patr*1244       _RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
2fa42a6013 Alis*1245
2b4c90c108 Mart*1246       INTEGER   iMin      , iMax          , jMin      , jMax
                1247       PARAMETER(iMin=2-OLx, iMax=sNx+OLx-1, jMin=2-OLy, jMax=sNy+OLy-1)
2fa42a6013 Alis*1248
956c0a5824 Patr*1249       _RL        p0    , p5    , p25     , p125      , p0625
2b4c90c108 Mart*1250       PARAMETER( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
2fa42a6013 Alis*1251
2b4c90c108 Mart*1252       DO j = jMin, jMax
2fa42a6013 Alis*1253          jm1 = j-1
                1254          jp1 = j+1
2b4c90c108 Mart*1255          DO i = iMin, iMax
2fa42a6013 Alis*1256             im1 = i-1
                1257             ip1 = i+1
                1258             tempVar =
11b9c75340 Jean*1259      &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
                1260      &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
                1261      &                     maskC(ip1,j  ,k,bi,bj)   +
                1262      &                     maskC(i  ,jm1,k,bi,bj)   +
                1263      &                     maskC(i  ,jp1,k,bi,bj) ) +
                1264      &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
                1265      &                     maskC(im1,jp1,k,bi,bj)   +
                1266      &                     maskC(ip1,jm1,k,bi,bj)   +
                1267      &                     maskC(ip1,jp1,k,bi,bj) )
2fa42a6013 Alis*1268             IF ( tempVar .GE. p25 ) THEN
                1269                fld_tmp(i,j) = (
11b9c75340 Jean*1270      &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
                1271      &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
                1272      &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
                1273      &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
                1274      &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
                1275      &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
                1276      &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
                1277      &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
                1278      &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
2fa42a6013 Alis*1279      &              / tempVar
                1280             ELSE
956c0a5824 Patr*1281                fld_tmp(i,j) = fld(i,j)
2fa42a6013 Alis*1282             ENDIF
                1283          ENDDO
                1284       ENDDO
                1285
                1286 c     transfer smoothed field to output array
2b4c90c108 Mart*1287       DO j = jMin, jMax
                1288          DO i = iMin, iMax
956c0a5824 Patr*1289             fld(i,j) = fld_tmp(i,j)
2fa42a6013 Alis*1290          ENDDO
                1291       ENDDO
                1292
                1293 #endif /* ALLOW_KPP */
                1294
2b4c90c108 Mart*1295       RETURN
                1296       END
2fa42a6013 Alis*1297
                1298 c*************************************************************************
                1299
                1300       subroutine blmix (
edb6656069 Mart*1301      I       ustar, bfsfc, hbl, stable, casea, diffus, kbl,
                1302      O       dkm1, blmc, ghat, sigma, ikey,
                1303      I       myThid )
2fa42a6013 Alis*1304
                1305 c     mixing coefficients within boundary layer depend on surface
                1306 c     forcing and the magnitude and gradient of interior mixing below
                1307 c     the boundary layer ("matching").
                1308 c
                1309 c     caution: if mixing bottoms out at hbl = -zgrid(Nr) then
                1310 c     fictitious layer at Nrp1 is needed with small but finite width
                1311 c     hwide(Nrp1) (eg. epsln = 1.e-20).
                1312 c
                1313       IMPLICIT NONE
                1314
                1315 #include "SIZE.h"
                1316 #include "KPP_PARAMS.h"
853c5e0e2c Jean*1317 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*1318 # include "tamc.h"
                1319 #endif
2fa42a6013 Alis*1320
                1321 c input
                1322 c     ustar (imt)                 surface friction velocity             (m/s)
                1323 c     bfsfc (imt)                 surface buoyancy forcing          (m^2/s^3)
                1324 c     hbl   (imt)                 boundary layer depth                    (m)
                1325 c     stable(imt)                 = 1 in stable forcing
                1326 c     casea (imt)                 = 1 in case A
                1327 c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
fe66051ebd Dimi*1328 c     kbl   (imt)                 -1 of first grid level below hbl
25c3477f99 Mart*1329 c     myThid               thread number for this instance of the routine
2b4c90c108 Mart*1330       INTEGER myThid
fe66051ebd Dimi*1331       _RL ustar (imt)
                1332       _RL bfsfc (imt)
                1333       _RL hbl   (imt)
                1334       _RL stable(imt)
                1335       _RL casea (imt)
                1336       _RL diffus(imt,0:Nrp1,mdiff)
2b4c90c108 Mart*1337       INTEGER kbl(imt)
2fa42a6013 Alis*1338
                1339 c output
                1340 c     dkm1 (imt,mdiff)            boundary layer difs at kbl-1 level
                1341 c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
                1342 c     ghat (imt,Nr)               nonlocal scalar transport
                1343 c     sigma(imt)                  normalized depth (d / hbl)
fe66051ebd Dimi*1344       _RL dkm1 (imt,mdiff)
                1345       _RL blmc (imt,Nr,mdiff)
                1346       _RL ghat (imt,Nr)
                1347       _RL sigma(imt)
edb6656069 Mart*1348       INTEGER ikey
2fa42a6013 Alis*1349
                1350 #ifdef ALLOW_KPP
                1351
                1352 c  local
1d478690dc Patr*1353 c     gat1*(imt)                 shape function at sigma = 1
                1354 c     dat1*(imt)                 derivative of shape function at sigma = 1
2fa42a6013 Alis*1355 c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
a9d2e4c565 Jean*1356       _RL gat1m(imt), gat1s(imt), gat1t(imt)
fe66051ebd Dimi*1357       _RL dat1m(imt), dat1s(imt), dat1t(imt)
                1358       _RL ws(imt), wm(imt)
2b4c90c108 Mart*1359       INTEGER i, kn, ki, kl
                1360 #ifndef KPP_DO_NOT_MATCH_DIFFUSIVITIES
                1361 # ifndef KPP_DO_NOT_MATCH_DERIVATIVES
                1362       _RL R, dvdzup, dvdzdn
                1363 # endif
                1364       _RL delhat
                1365 #endif
                1366       _RL viscp, difsp, diftp, visch, difsh, difth
                1367       _RL f1, sig, a1, a2, a3
fe66051ebd Dimi*1368       _RL Gm, Gs, Gt
                1369       _RL tempVar
2fa42a6013 Alis*1370
fe66051ebd Dimi*1371       _RL    p0    , eins
2b4c90c108 Mart*1372       PARAMETER (p0=0.0, eins=1.0)
88b144ae49 Jean*1373 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1374       INTEGER kkey
88b144ae49 Jean*1375 #endif
2fa42a6013 Alis*1376
                1377 c-----------------------------------------------------------------------
                1378 c compute velocity scales at hbl
                1379 c-----------------------------------------------------------------------
                1380
2b4c90c108 Mart*1381       DO i = 1, imt
2fa42a6013 Alis*1382          sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
2b4c90c108 Mart*1383       ENDDO
2fa42a6013 Alis*1384
853c5e0e2c Jean*1385 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1386 CADJ STORE sigma = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1387 #endif
2b4c90c108 Mart*1388       CALL wscale (
2fa42a6013 Alis*1389      I        sigma, hbl, ustar, bfsfc,
25c3477f99 Mart*1390      O        wm, ws, myThid )
853c5e0e2c Jean*1391 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1392 CADJ STORE wm = comlev1_kpp, key=ikey, kind=isbyte
                1393 CADJ STORE ws = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1394 #endif
2fa42a6013 Alis*1395
2b4c90c108 Mart*1396       DO i = 1, imt
                1397          wm(i) = sign(eins,wm(i))*MAX(phepsi,ABS(wm(i)))
                1398          ws(i) = sign(eins,ws(i))*MAX(phepsi,ABS(ws(i)))
                1399       ENDDO
853c5e0e2c Jean*1400 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1401 CADJ STORE wm = comlev1_kpp, key=ikey, kind=isbyte
                1402 CADJ STORE ws = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1403 #endif
956c0a5824 Patr*1404
2b4c90c108 Mart*1405       DO i = 1, imt
2fa42a6013 Alis*1406
2b4c90c108 Mart*1407          kn = INT(caseA(i)+phepsi) *(kbl(i) -1) +
                1408      $        (1 - INT(caseA(i)+phepsi)) * kbl(i)
2fa42a6013 Alis*1409
                1410 c-----------------------------------------------------------------------
                1411 c find the interior viscosities and derivatives at hbl(i)
                1412 c-----------------------------------------------------------------------
                1413
2b4c90c108 Mart*1414 C     initialise diffusivities (*h) and derivatives (*p)
                1415 #ifdef KPP_DO_NOT_MATCH_DIFFUSIVITIES
                1416          visch  = 0.
                1417          difsh  = 0.
                1418          difth  = 0.
                1419          viscp  = 0.
                1420          difsp  = 0.
                1421          diftp  = 0.
                1422 #else /* DO_MATCH_DIFFUSIVITIES */
                1423 # ifdef KPP_DO_NOT_MATCH_DERIVATIVES
                1424          viscp  = 0.
                1425          difsp  = 0.
                1426          diftp  = 0.
                1427          delhat = 0.
                1428 # else /*  DO_MATCH_DERIVATIVES */
2fa42a6013 Alis*1429          delhat = 0.5*hwide(kn) - zgrid(kn) - hbl(i)
                1430          R      = 1.0 - delhat / hwide(kn)
                1431          dvdzup = (diffus(i,kn-1,1) - diffus(i,kn  ,1)) / hwide(kn)
                1432          dvdzdn = (diffus(i,kn  ,1) - diffus(i,kn+1,1)) / hwide(kn+1)
2b4c90c108 Mart*1433 CML   This is the same as:
                1434 CML      viscp  = (1.-R) * MAX(dvdzup,0.) + R  * MAX(dvdzdn,0.)
                1435 CML   why do we need that here?
                1436          viscp  = 0.5 * ( (1.-R) * (dvdzup + ABS(dvdzup)) +
                1437      1                        R  * (dvdzdn + ABS(dvdzdn))  )
2fa42a6013 Alis*1438
                1439          dvdzup = (diffus(i,kn-1,2) - diffus(i,kn  ,2)) / hwide(kn)
                1440          dvdzdn = (diffus(i,kn  ,2) - diffus(i,kn+1,2)) / hwide(kn+1)
2b4c90c108 Mart*1441          difsp  = 0.5 * ( (1.-R) * (dvdzup + ABS(dvdzup)) +
                1442      1                        R  * (dvdzdn + ABS(dvdzdn))  )
2fa42a6013 Alis*1443
                1444          dvdzup = (diffus(i,kn-1,3) - diffus(i,kn  ,3)) / hwide(kn)
                1445          dvdzdn = (diffus(i,kn  ,3) - diffus(i,kn+1,3)) / hwide(kn+1)
2b4c90c108 Mart*1446          diftp  = 0.5 * ( (1.-R) * (dvdzup + ABS(dvdzup)) +
                1447      1                        R  * (dvdzdn + ABS(dvdzdn))  )
                1448 # endif /* KPP_DO_NOT_MATCH_DERIVATIVES */
2fa42a6013 Alis*1449          visch  = diffus(i,kn,1) + viscp * delhat
                1450          difsh  = diffus(i,kn,2) + difsp * delhat
                1451          difth  = diffus(i,kn,3) + diftp * delhat
2b4c90c108 Mart*1452 #endif /* KPP_DO_NOT_MATCH_DIFFUSIVITIES */
2fa42a6013 Alis*1453
a9d2e4c565 Jean*1454          f1 = stable(i) * conc1 * bfsfc(i) /
2b4c90c108 Mart*1455 #ifdef KPP_SMOOTH_REGULARISATION
                1456      &        (ustar(i)**4 + phepsi)
                1457 #else
                1458      &        MAX(ustar(i)**4,phepsi)
                1459 #endif
1d478690dc Patr*1460          gat1m(i) = visch / hbl(i) / wm(i)
                1461          dat1m(i) = -viscp / wm(i) + f1 * visch
a9d2e4c565 Jean*1462
1d478690dc Patr*1463          gat1s(i) = difsh  / hbl(i) / ws(i)
a9d2e4c565 Jean*1464          dat1s(i) = -difsp / ws(i) + f1 * difsh
                1465
1d478690dc Patr*1466          gat1t(i) = difth /  hbl(i) / ws(i)
a9d2e4c565 Jean*1467          dat1t(i) = -diftp / ws(i) + f1 * difth
                1468
2b4c90c108 Mart*1469       ENDDO
853c5e0e2c Jean*1470 #ifdef KPP_AUTODIFF_MORE_STORE
edb6656069 Mart*1471 CADJ STORE gat1m = comlev1_kpp, key=ikey, kind=isbyte
                1472 CADJ STORE gat1s = comlev1_kpp, key=ikey, kind=isbyte
                1473 CADJ STORE gat1t = comlev1_kpp, key=ikey, kind=isbyte
                1474 CADJ STORE dat1m = comlev1_kpp, key=ikey, kind=isbyte
                1475 CADJ STORE dat1s = comlev1_kpp, key=ikey, kind=isbyte
                1476 CADJ STORE dat1t = comlev1_kpp, key=ikey, kind=isbyte
3067745c9b Patr*1477 #endif
2b4c90c108 Mart*1478       DO i = 1, imt
                1479          dat1m(i) = MIN(dat1m(i),p0)
                1480          dat1s(i) = MIN(dat1s(i),p0)
                1481          dat1t(i) = MIN(dat1t(i),p0)
                1482       ENDDO
853c5e0e2c Jean*1483 #ifdef KPP_AUTODIFF_MORE_STORE
edb6656069 Mart*1484 CADJ STORE dat1m = comlev1_kpp, key=ikey, kind=isbyte
                1485 CADJ STORE dat1s = comlev1_kpp, key=ikey, kind=isbyte
                1486 CADJ STORE dat1t = comlev1_kpp, key=ikey, kind=isbyte
3067745c9b Patr*1487 #endif
2fa42a6013 Alis*1488
2b4c90c108 Mart*1489       DO ki = 1, Nr
2fa42a6013 Alis*1490
a9d2e4c565 Jean*1491 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1492          kkey = (ikey-1)*Nr + ki
3067745c9b Patr*1493 #endif
                1494
2fa42a6013 Alis*1495 c-----------------------------------------------------------------------
                1496 c     compute turbulent velocity scales on the interfaces
                1497 c-----------------------------------------------------------------------
                1498
2b4c90c108 Mart*1499          DO i = 1, imt
2fa42a6013 Alis*1500             sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
2b4c90c108 Mart*1501             sigma(i) = stable(i)*sig + (1.-stable(i))*MIN(sig,epsilon)
                1502          ENDDO
853c5e0e2c Jean*1503 #ifdef KPP_AUTODIFF_MORE_STORE
edb6656069 Mart*1504 CADJ STORE wm = comlev1_kpp_k, key = kkey
                1505 CADJ STORE ws = comlev1_kpp_k, key = kkey
3067745c9b Patr*1506 #endif
853c5e0e2c Jean*1507 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1508 CADJ STORE sigma = comlev1_kpp_k, key = kkey
853c5e0e2c Jean*1509 #endif
2b4c90c108 Mart*1510          CALL wscale (
2fa42a6013 Alis*1511      I        sigma, hbl, ustar, bfsfc,
25c3477f99 Mart*1512      O        wm, ws, myThid )
853c5e0e2c Jean*1513 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1514 CADJ STORE wm = comlev1_kpp_k, key = kkey
                1515 CADJ STORE ws = comlev1_kpp_k, key = kkey
853c5e0e2c Jean*1516 #endif
2fa42a6013 Alis*1517
                1518 c-----------------------------------------------------------------------
                1519 c     compute the dimensionless shape functions at the interfaces
                1520 c-----------------------------------------------------------------------
                1521
2b4c90c108 Mart*1522          DO i = 1, imt
2fa42a6013 Alis*1523             sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
                1524             a1 = sig - 2.
                1525             a2 = 3. - 2. * sig
                1526             a3 = sig - 1.
                1527
1d478690dc Patr*1528             Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
                1529             Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
                1530             Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
2fa42a6013 Alis*1531
                1532 c-----------------------------------------------------------------------
                1533 c     compute boundary layer diffusivities at the interfaces
                1534 c-----------------------------------------------------------------------
                1535
1d478690dc Patr*1536             blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
2fa42a6013 Alis*1537             blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
                1538             blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
                1539
                1540 c-----------------------------------------------------------------------
                1541 c     nonlocal transport term = ghat * <ws>o
                1542 c-----------------------------------------------------------------------
956c0a5824 Patr*1543
                1544             tempVar = ws(i) * hbl(i)
2b4c90c108 Mart*1545 #ifdef KPP_SMOOTH_REGULARISATION
                1546             ghat(i,ki) = (1.-stable(i)) * cg / (phepsi+tempVar)
                1547 #else
                1548             ghat(i,ki) = (1.-stable(i)) * cg / MAX(phepsi,tempVar)
                1549 #endif
                1550          ENDDO
                1551       ENDDO
2fa42a6013 Alis*1552
                1553 c-----------------------------------------------------------------------
                1554 c find diffusivities at kbl-1 grid level
                1555 c-----------------------------------------------------------------------
                1556
2b4c90c108 Mart*1557       DO i = 1, imt
80b2748a09 Patr*1558          kl = kbl(i)
                1559          sig      = -zgrid(kl-1) / hbl(i)
a9d2e4c565 Jean*1560          sigma(i) = stable(i) * sig
2b4c90c108 Mart*1561      &            + (1. - stable(i)) * MIN(sig,epsilon)
                1562       ENDDO
2fa42a6013 Alis*1563
853c5e0e2c Jean*1564 #ifdef KPP_AUTODIFF_MORE_STORE
edb6656069 Mart*1565 CADJ STORE wm = comlev1_kpp, key=ikey, kind=isbyte
                1566 CADJ STORE ws = comlev1_kpp, key=ikey, kind=isbyte
3067745c9b Patr*1567 #endif
853c5e0e2c Jean*1568 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1569 CADJ STORE sigma = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1570 #endif
2b4c90c108 Mart*1571       CALL wscale (
2fa42a6013 Alis*1572      I        sigma, hbl, ustar, bfsfc,
25c3477f99 Mart*1573      O        wm, ws, myThid )
853c5e0e2c Jean*1574 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*1575 CADJ STORE wm = comlev1_kpp, key=ikey, kind=isbyte
                1576 CADJ STORE ws = comlev1_kpp, key=ikey, kind=isbyte
853c5e0e2c Jean*1577 #endif
2fa42a6013 Alis*1578
2b4c90c108 Mart*1579       DO i = 1, imt
80b2748a09 Patr*1580          kl = kbl(i)
                1581          sig = -zgrid(kl-1) / hbl(i)
2fa42a6013 Alis*1582          a1 = sig - 2.
                1583          a2 = 3. - 2. * sig
                1584          a3 = sig - 1.
1d478690dc Patr*1585          Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
                1586          Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
                1587          Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
                1588          dkm1(i,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
2fa42a6013 Alis*1589          dkm1(i,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
                1590          dkm1(i,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
2b4c90c108 Mart*1591       ENDDO
2fa42a6013 Alis*1592
                1593 #endif /* ALLOW_KPP */
                1594
2b4c90c108 Mart*1595       RETURN
                1596       END
2fa42a6013 Alis*1597
                1598 c*************************************************************************
                1599
a9d2e4c565 Jean*1600       subroutine enhance (
edb6656069 Mart*1601      I       dkm1, hbl, kbl, diffus, casea,
                1602      U       ghat,
                1603      O       blmc,
                1604      &       myThid )
2fa42a6013 Alis*1605
                1606 c enhance the diffusivity at the kbl-.5 interface
                1607
                1608       IMPLICIT NONE
                1609
                1610 #include "SIZE.h"
                1611 #include "KPP_PARAMS.h"
                1612
                1613 c input
                1614 c     dkm1(imt,mdiff)          bl diffusivity at kbl-1 grid level
                1615 c     hbl(imt)                  boundary layer depth                 (m)
                1616 c     kbl(imt)                  grid above hbl
                1617 c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
                1618 c     casea(imt)                = 1 in caseA, = 0 in case B
25c3477f99 Mart*1619 c     myThid                    thread number for this instance of the routine
2b4c90c108 Mart*1620       INTEGER   myThid
fe66051ebd Dimi*1621       _RL dkm1  (imt,mdiff)
                1622       _RL hbl   (imt)
2b4c90c108 Mart*1623       INTEGER kbl   (imt)
fe66051ebd Dimi*1624       _RL diffus(imt,0:Nrp1,mdiff)
                1625       _RL casea (imt)
2fa42a6013 Alis*1626
                1627 c input/output
                1628 c     nonlocal transport, modified ghat at kbl(i)-1 interface    (s/m**2)
fe66051ebd Dimi*1629       _RL ghat (imt,Nr)
2fa42a6013 Alis*1630
                1631 c output
                1632 c     enhanced bound. layer mixing coeff.
fe66051ebd Dimi*1633       _RL blmc  (imt,Nr,mdiff)
2fa42a6013 Alis*1634
                1635 #ifdef ALLOW_KPP
                1636
                1637 c local
                1638 c     fraction hbl lies beteen zgrid neighbors
fe66051ebd Dimi*1639       _RL delta
2b4c90c108 Mart*1640       INTEGER ki, i, md
fe66051ebd Dimi*1641       _RL dkmp5, dstar
2fa42a6013 Alis*1642
2b4c90c108 Mart*1643       DO i = 1, imt
2fa42a6013 Alis*1644          ki = kbl(i)-1
2b4c90c108 Mart*1645          IF ((ki .ge. 1) .AND. (ki .LT. Nr)) THEN
2fa42a6013 Alis*1646             delta = (hbl(i) + zgrid(ki)) / (zgrid(ki) - zgrid(ki+1))
2b4c90c108 Mart*1647             DO md = 1, mdiff
2fa42a6013 Alis*1648                dkmp5         =      casea(i)  * diffus(i,ki,md) +
                1649      1                         (1.- casea(i)) * blmc  (i,ki,md)
2b4c90c108 Mart*1650 C     I think that this is meant here, but I cannot be sure because there
                1651 C     there is no reference for this. In MOM6/CVMix, the square is outside
                1652 C     of the parentheses
                1653 CML               dstar         = (1.- delta**2) * dkm1(i,md)
2fa42a6013 Alis*1654                dstar         = (1.- delta)**2 * dkm1(i,md)
                1655      &                       + delta**2 * dkmp5
                1656                blmc(i,ki,md) = (1.- delta)*diffus(i,ki,md)
                1657      &                       + delta*dstar
2b4c90c108 Mart*1658             ENDDO
2fa42a6013 Alis*1659             ghat(i,ki) = (1.- casea(i)) * ghat(i,ki)
2b4c90c108 Mart*1660          ENDIF
                1661       ENDDO
2fa42a6013 Alis*1662
                1663 #endif /* ALLOW_KPP */
                1664
2b4c90c108 Mart*1665       RETURN
                1666       END
2fa42a6013 Alis*1667
                1668 c*************************************************************************
                1669
                1670       SUBROUTINE STATEKPP (
25c3477f99 Mart*1671      O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
edb6656069 Mart*1672      I     ikey, bi, bj, myThid )
2fa42a6013 Alis*1673 c
                1674 c-----------------------------------------------------------------------
                1675 c     "statekpp" computes all necessary input arrays
                1676 c     for the kpp mixing scheme
                1677 c
                1678 c     input:
                1679 c      bi, bj = array indices on which to apply calculations
                1680 c
                1681 c     output:
                1682 c      rho1   = potential density of surface layer                     (kg/m^3)
                1683 c      dbloc  = local buoyancy gradient at Nr interfaces
                1684 c               g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ]          (m/s^2)
                1685 c      dbsfc  = buoyancy difference with respect to the surface
                1686 c               g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ]         (m/s^2)
                1687 c      ttalpha= thermal expansion coefficient without 1/rho factor
                1688 c               d(rho) / d(potential temperature)                    (kg/m^3/C)
                1689 c      ssbeta = salt expansion coefficient without 1/rho factor
ba0b047096 Mart*1690 c               d(rho) / d(salinity)                          ((kg/m^3)/(g/kg))
2fa42a6013 Alis*1691 c
                1692 c     see also subroutines find_rho.F find_alpha.F find_beta.F
                1693 c
                1694 c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
                1695 c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV
                1696 c
7e819019d5 Dimi*1697
2fa42a6013 Alis*1698 c-----------------------------------------------------------------------
                1699
                1700       IMPLICIT NONE
                1701
                1702 #include "SIZE.h"
                1703 #include "EEPARAMS.h"
                1704 #include "PARAMS.h"
                1705 #include "KPP_PARAMS.h"
42c525bfb4 Alis*1706 #include "DYNVARS.h"
7e819019d5 Dimi*1707 #include "GRID.h"
853c5e0e2c Jean*1708 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*1709 # include "tamc.h"
                1710 #endif
2fa42a6013 Alis*1711
                1712 c-------------- Routine arguments -----------------------------------------
                1713       INTEGER bi, bj, myThid
fe66051ebd Dimi*1714       _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
                1715       _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
                1716       _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
                1717       _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
                1718       _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
2fa42a6013 Alis*1719
                1720 #ifdef ALLOW_KPP
                1721
                1722 c--------------------------------------------------------------------------
                1723 c
                1724 c     local arrays:
                1725 c
                1726 c     rhok         - density of t(k  ) & s(k  ) at depth k
                1727 c     rhokm1       - density of t(k-1) & s(k-1) at depth k
                1728 c     rho1k        - density of t(1  ) & s(1  ) at depth k
7e819019d5 Dimi*1729 c     work1,2,3    - work arrays for holding horizontal slabs
2fa42a6013 Alis*1730
                1731       _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1732       _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1733       _RL RHO1K (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1734       _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1735       _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1736       _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7e819019d5 Dimi*1737
2b4c90c108 Mart*1738       INTEGER i, j, k
edb6656069 Mart*1739       INTEGER ikey, kkey
2fa42a6013 Alis*1740
                1741 c calculate density, alpha, beta in surface layer, and set dbsfc to zero
                1742
edb6656069 Mart*1743       kkey = (ikey-1)*Nr + 1
3067745c9b Patr*1744
2b4c90c108 Mart*1745       k = 1
853c5e0e2c Jean*1746 #ifdef KPP_AUTODIFF_MORE_STORE
2b4c90c108 Mart*1747 CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1748 CADJ &     key=kkey, kind=isbyte
2b4c90c108 Mart*1749 CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1750 CADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1751 #endif /* KPP_AUTODIFF_MORE_STORE */
94c8eb5701 Jean*1752       CALL FIND_RHO_2D(
                1753      I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
2b4c90c108 Mart*1754      I     theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
2fa42a6013 Alis*1755      O     WORK1,
2b4c90c108 Mart*1756      I     k, bi, bj, myThid )
853c5e0e2c Jean*1757 #ifdef KPP_AUTODIFF_MORE_STORE
2b4c90c108 Mart*1758 CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1759 CADJ &     key=kkey, kind=isbyte
2b4c90c108 Mart*1760 CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1761 CADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1762 #endif /* KPP_AUTODIFF_MORE_STORE */
2fa42a6013 Alis*1763
2b4c90c108 Mart*1764       CALL FIND_ALPHA(
4fa51dc730 Dimi*1765      I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
25c3477f99 Mart*1766      O     WORK2, myThid )
2fa42a6013 Alis*1767
2b4c90c108 Mart*1768       CALL FIND_BETA(
4fa51dc730 Dimi*1769      I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
25c3477f99 Mart*1770      O     WORK3, myThid )
2fa42a6013 Alis*1771
2b4c90c108 Mart*1772       DO j = 1-OLy, sNy+OLy
                1773        DO i = 1-OLx, sNx+OLx
                1774         RHO1(i,j)      = WORK1(i,j) + rhoConst
                1775         TTALPHA(i,j,1) = WORK2(i,j)
                1776         SSBETA(i,j,1)  = WORK3(i,j)
                1777         DBSFC(i,j,1)   = 0.
                1778        ENDDO
                1779       ENDDO
2fa42a6013 Alis*1780
                1781 c calculate alpha, beta, and gradients in interior layers
                1782
1d478690dc Patr*1783 CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
2b4c90c108 Mart*1784       DO k = 2, Nr
2fa42a6013 Alis*1785
edb6656069 Mart*1786          kkey = (ikey-1)*Nr + k
3067745c9b Patr*1787
853c5e0e2c Jean*1788 #ifdef KPP_AUTODIFF_MORE_STORE
95c72ef3a1 Patr*1789 CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1790 CADJ &     key=kkey, kind=isbyte
95c72ef3a1 Patr*1791 CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1792 CADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1793 #endif /* KPP_AUTODIFF_MORE_STORE */
94c8eb5701 Jean*1794          CALL FIND_RHO_2D(
                1795      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
                1796      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
2fa42a6013 Alis*1797      O        RHOK,
94c8eb5701 Jean*1798      I        k, bi, bj, myThid )
2fa42a6013 Alis*1799
853c5e0e2c Jean*1800 #ifdef KPP_AUTODIFF_MORE_STORE
95c72ef3a1 Patr*1801 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1802 CADJ &     key=kkey, kind=isbyte
95c72ef3a1 Patr*1803 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1804 CADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1805 #endif /* KPP_AUTODIFF_MORE_STORE */
94c8eb5701 Jean*1806          CALL FIND_RHO_2D(
                1807      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
                1808      I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
2fa42a6013 Alis*1809      O        RHOKM1,
94c8eb5701 Jean*1810      I        k-1, bi, bj, myThid )
2fa42a6013 Alis*1811
853c5e0e2c Jean*1812 #ifdef KPP_AUTODIFF_MORE_STORE
95c72ef3a1 Patr*1813 CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1814 CADJ &     key=kkey, kind=isbyte
95c72ef3a1 Patr*1815 CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1816 CADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1817 #endif /* KPP_AUTODIFF_MORE_STORE */
94c8eb5701 Jean*1818          CALL FIND_RHO_2D(
                1819      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
                1820      I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
2fa42a6013 Alis*1821      O        RHO1K,
94c8eb5701 Jean*1822      I        1, bi, bj, myThid )
2fa42a6013 Alis*1823
853c5e0e2c Jean*1824 #ifdef KPP_AUTODIFF_MORE_STORE
edb6656069 Mart*1825 CADJ STORE rhok  (:,:) = comlev1_kpp_k, key=kkey, kind=isbyte
                1826 CADJ STORE rhokm1(:,:) = comlev1_kpp_k, key=kkey, kind=isbyte
                1827 CADJ STORE rho1k (:,:) = comlev1_kpp_k, key=kkey, kind=isbyte
853c5e0e2c Jean*1828 #endif /* KPP_AUTODIFF_MORE_STORE */
3067745c9b Patr*1829
2b4c90c108 Mart*1830          CALL FIND_ALPHA(
                1831      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k, k,
25c3477f99 Mart*1832      O        WORK1, myThid )
2fa42a6013 Alis*1833
2b4c90c108 Mart*1834          CALL FIND_BETA(
                1835      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k, k,
25c3477f99 Mart*1836      O        WORK2, myThid )
2fa42a6013 Alis*1837
2b4c90c108 Mart*1838          DO j = 1-OLy, sNy+OLy
                1839             DO i = 1-OLx, sNx+OLx
                1840                TTALPHA(i,j,k) = WORK1 (i,j)
                1841                SSBETA(i,j,k)  = WORK2 (i,j)
                1842                DBLOC(i,j,k-1) = gravity * (RHOK(i,j) - RHOKM1(i,j)) /
                1843      &                                    (RHOK(i,j) + rhoConst)
                1844                DBSFC(i,j,k)   = gravity * (RHOK(i,j) - RHO1K (i,j)) /
                1845      &                                    (RHOK(i,j) + rhoConst)
                1846             ENDDO
                1847          ENDDO
                1848
                1849       ENDDO
                1850
                1851 c     compute arrays for k = Nrp1
                1852       DO j = 1-OLy, sNy+OLy
                1853        DO i = 1-OLx, sNx+OLx
                1854         TTALPHA(i,j,Nrp1) = TTALPHA(i,j,Nr)
                1855         SSBETA(i,j,Nrp1)  = SSBETA(i,j,Nr)
                1856         DBLOC(i,j,Nr)     = 0.
                1857        ENDDO
                1858       ENDDO
2fa42a6013 Alis*1859
7e819019d5 Dimi*1860 #ifdef ALLOW_DIAGNOSTICS
                1861       IF ( useDiagnostics ) THEN
2b4c90c108 Mart*1862        CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,2,bi,bj,myThid)
                1863        CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,2,bi,bj,myThid)
7e819019d5 Dimi*1864       ENDIF
                1865 #endif /* ALLOW_DIAGNOSTICS */
                1866
2fa42a6013 Alis*1867 #endif /* ALLOW_KPP */
                1868
                1869       RETURN
                1870       END
e750a5e49e Mart*1871
                1872 c*************************************************************************
                1873
                1874       SUBROUTINE KPP_DOUBLEDIFF (
                1875      I     TTALPHA, SSBETA,
a9d2e4c565 Jean*1876      U     kappaRT,
e750a5e49e Mart*1877      U     kappaRS,
edb6656069 Mart*1878      I     ikey, iMin, iMax, jMin, jMax, bi, bj, myThid )
e750a5e49e Mart*1879 c
                1880 c-----------------------------------------------------------------------
a9d2e4c565 Jean*1881 c     "KPP_DOUBLEDIFF" adds the double diffusive contributions
                1882 C     as Rrho-dependent parameterizations to kappaRT and kappaRS
e750a5e49e Mart*1883 c
                1884 c     input:
                1885 c     bi, bj  = array indices on which to apply calculations
2b4c90c108 Mart*1886 c     iMin, iMax, jMin, jMax = array boundaries
edb6656069 Mart*1887 c     ikey = key for TAMC/TAF automatic differentiation
e750a5e49e Mart*1888 c     myThid  = thread id
                1889 c
                1890 c      ttalpha= thermal expansion coefficient without 1/rho factor
                1891 c               d(rho) / d(potential temperature)                    (kg/m^3/C)
                1892 c      ssbeta = salt expansion coefficient without 1/rho factor
ba0b047096 Mart*1893 c               d(rho) / d(salinity)                          ((kg/m^3)/(g/kg))
e750a5e49e Mart*1894 c     output: updated
                1895 c     kappaRT/S :: background diffusivities for temperature and salinity
                1896 c
a9d2e4c565 Jean*1897 c     written  by: martin losch,   sept. 15, 2009
e750a5e49e Mart*1898 c
                1899
                1900 c-----------------------------------------------------------------------
                1901
                1902       IMPLICIT NONE
                1903
                1904 #include "SIZE.h"
                1905 #include "EEPARAMS.h"
                1906 #include "PARAMS.h"
                1907 #include "KPP_PARAMS.h"
                1908 #include "DYNVARS.h"
                1909 #include "GRID.h"
853c5e0e2c Jean*1910 #ifdef ALLOW_AUTODIFF_TAMC
e750a5e49e Mart*1911 # include "tamc.h"
                1912 #endif
                1913
                1914 c-------------- Routine arguments -----------------------------------------
edb6656069 Mart*1915       INTEGER ikey, iMin, iMax, jMin, jMax, bi, bj, myThid
e750a5e49e Mart*1916
                1917       _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
                1918       _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
52b77e9710 Jean*1919       _RL KappaRT( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
                1920       _RL KappaRS( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
e750a5e49e Mart*1921
                1922 #ifdef ALLOW_KPP
                1923
                1924 C--------------------------------------------------------------------------
                1925 C
                1926 C     local variables
2b4c90c108 Mart*1927 C     i,j,k :: loop indices
edb6656069 Mart*1928 C     kkey  :: key for TAMC/TAF automatic differentiation
e750a5e49e Mart*1929 C
2b4c90c108 Mart*1930       INTEGER i, j, k
edb6656069 Mart*1931       INTEGER kkey
e750a5e49e Mart*1932 C     alphaDT   :: d\rho/d\theta * d\theta
                1933 C     betaDS    :: d\rho/dsalt * dsalt
                1934 C     Rrho      :: "density ratio" R_{\rho} = \alpha dT/dz / \beta dS/dz
                1935 C     nuddt/s   :: double diffusive diffusivities
                1936 C     numol     :: molecular diffusivity
                1937 C     rFac      :: abbreviation for 1/(R_{\rho0}-1)
                1938
                1939       _RL alphaDT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
                1940       _RL betaDS  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
                1941       _RL nuddt   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
                1942       _RL nudds   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
                1943       _RL Rrho
                1944       _RL numol, rFac, nutmp
                1945       INTEGER Km1
                1946
                1947 C     set some constants here
                1948       numol = 1.5 _d -06
                1949       rFac  = 1. _d 0 / (Rrho0 - 1. _d 0 )
                1950 C
edb6656069 Mart*1951       kkey = (ikey-1)*Nr + 1
e750a5e49e Mart*1952
853c5e0e2c Jean*1953 CML#ifdef KPP_AUTODIFF_MORE_STORE
e750a5e49e Mart*1954 CMLCADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1955 CMLCADJ &     key=kkey, kind=isbyte
e750a5e49e Mart*1956 CMLCADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
edb6656069 Mart*1957 CMLCADJ &     key=kkey, kind=isbyte
853c5e0e2c Jean*1958 CML#endif /* KPP_AUTODIFF_MORE_STORE */
e750a5e49e Mart*1959
2b4c90c108 Mart*1960       DO k = 1, Nr
                1961        Km1 = MAX(k-1,1)
                1962        DO j = 1-OLy, sNy+OLy
                1963         DO i = 1-OLx, sNx+OLx
                1964          alphaDT(i,j) = ( theta(i,j,km1,bi,bj)-theta(i,j,k,bi,bj) )
                1965      &        * 0.5 _d 0 * ABS( TTALPHA(i,j,km1) + TTALPHA(i,j,k) )
                1966          betaDS(i,j)  = ( salt(i,j,km1,bi,bj)-salt(i,j,k,bi,bj) )
                1967      &        * 0.5 _d 0 * ( SSBETA(i,j,km1) + SSBETA(i,j,k) )
                1968          nuddt(i,j) = 0. _d 0
                1969          nudds(i,j) = 0. _d 0
e750a5e49e Mart*1970         ENDDO
                1971        ENDDO
2b4c90c108 Mart*1972        IF ( k .GT. 1 ) THEN
                1973         DO j = jMin, jMax
                1974          DO i = iMin, iMax
e750a5e49e Mart*1975           Rrho  = 0. _d 0
                1976 C     Now we have many different cases
a9d2e4c565 Jean*1977 C     a. alphaDT > 0 and betaDS > 0 => salt fingering
e750a5e49e Mart*1978 C        (salinity destabilizes)
2b4c90c108 Mart*1979           IF (      alphaDT(i,j) .GT. betaDS(i,j)
                1980      &         .AND. betaDS(i,j) .GT. 0. _d 0 ) THEN
                1981            Rrho = MIN( alphaDT(i,j)/betaDS(i,j), Rrho0 )
e750a5e49e Mart*1982 C     Large et al. 1994, eq. 31a
2b4c90c108 Mart*1983 C          nudds(i,j) = dsfmax * ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )**3
e750a5e49e Mart*1984            nutmp      =          ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )
2b4c90c108 Mart*1985            nudds(i,j) = dsfmax * nutmp * nutmp * nutmp
e750a5e49e Mart*1986 C     Large et al. 1994, eq. 31c
2b4c90c108 Mart*1987            nuddt(i,j) = 0.7 _d 0 * nudds(i,j)
                1988           ELSEIF (   alphaDT(i,j) .LT. 0. _d 0
                1989      &          .AND. betaDS(i,j) .LT. 0. _d 0
                1990      &          .AND.alphaDT(i,j) .GT. betaDS(i,j) ) THEN
e750a5e49e Mart*1991 C     b. alphaDT < 0 and betaDS < 0 => semi-convection, diffusive convection
                1992 C        (temperature destabilizes)
                1993 C     for Rrho >= 1 the water column is statically unstable and we never
                1994 C     reach this point
2b4c90c108 Mart*1995            Rrho = alphaDT(i,j)/betaDS(i,j)
e750a5e49e Mart*1996 C     Large et al. 1994, eq. 32
2b4c90c108 Mart*1997            nuddt(i,j) = numol * 0.909 _d 0
a9d2e4c565 Jean*1998      &          * exp ( 4.6 _d 0 * exp (
e750a5e49e Mart*1999      &          - 5.4 _d 0 * ( 1. _d 0/Rrho - 1. _d 0 ) ) )
                2000 CMLC     or
                2001 CMLC     Large et al. 1994, eq. 33
2b4c90c108 Mart*2002 CML         nuddt(i,j) = numol * 8.7 _d 0 * Rrho**1.1
e750a5e49e Mart*2003 C     Large et al. 1994, eqs. 34
2b4c90c108 Mart*2004            nudds(i,j) = nuddt(i,j) * MAX( 0.15 _d 0 * Rrho,
e750a5e49e Mart*2005      &          1.85 _d 0 * Rrho - 0.85 _d 0 )
                2006           ELSE
a9d2e4c565 Jean*2007 C     Do nothing, because in this case the water colume is unstable
                2008 C     => double diffusive processes are negligible and mixing due
                2009 C     to shear instability will dominate
e750a5e49e Mart*2010           ENDIF
                2011          ENDDO
                2012         ENDDO
2b4c90c108 Mart*2013 C     ENDIF ( k .GT. 1 )
e750a5e49e Mart*2014        ENDIF
                2015 C
2b4c90c108 Mart*2016        DO j = 1-OLy, sNy+OLy
                2017         DO i = 1-OLx, sNx+OLx
                2018          kappaRT(i,j,k) = kappaRT(i,j,k) + nuddt(i,j)
                2019          kappaRS(i,j,k) = kappaRS(i,j,k) + nudds(i,j)
e750a5e49e Mart*2020         ENDDO
                2021        ENDDO
                2022 #ifdef ALLOW_DIAGNOSTICS
                2023        IF ( useDiagnostics ) THEN
                2024         CALL DIAGNOSTICS_FILL(nuddt,'KPPnuddt',k,1,2,bi,bj,myThid)
                2025         CALL DIAGNOSTICS_FILL(nudds,'KPPnudds',k,1,2,bi,bj,myThid)
                2026        ENDIF
                2027 #endif /* ALLOW_DIAGNOSTICS */
2b4c90c108 Mart*2028 C     end of k-loop
e750a5e49e Mart*2029       ENDDO
                2030 #endif /* ALLOW_KPP */
                2031
                2032       RETURN
                2033       END