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