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