C $Header: /u/gcmpack/MITgcm/pkg/kpp/kpp_routines.F,v 1.57 2016/11/30 00:14:14 jmc Exp $
C $Name: $
#include "KPP_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
#ifdef ALLOW_SALT_PLUME
#include "SALT_PLUME_OPTIONS.h"
#endif
#if (defined ALLOW_AUTODIFF_TAMC) (defined KPP_AUTODIFF_EXCESSIVE_STORE)
# define KPP_AUTODIFF_MORE_STORE
#endif
C-- File kpp_routines.F: subroutines needed to implement
C-- KPP vertical mixing scheme
C-- Contents
C-- o KPPMIX - Main driver and interface routine.
C-- o BLDEPTH - Determine oceanic planetary boundary layer depth.
C-- o WSCALE - Compute turbulent velocity scales.
C-- o RI_IWMIX - Compute interior viscosity diffusivity coefficients.
C-- o Z121 - Apply 121 vertical smoothing.
C-- o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
C-- o BLMIX - Boundary layer mixing coefficients.
C-- o ENHANCE - Enhance diffusivity at boundary layer interface.
C-- o STATEKPP - Compute buoyancy-related input arrays.
C-- o KPP_DOUBLEDIFF - Compute double diffusive contribution to diffusivities
c***********************************************************************
SUBROUTINE KPPMIX (
I kmtj, shsq, dvsq, ustar, msk
I , bo, bosol
#ifdef ALLOW_SALT_PLUME
I , boplume,SPDepth
#ifdef SALT_PLUME_SPLIT_BASIN
I , lon,lat
#endif /* SALT_PLUME_SPLIT_BASIN */
#endif /* ALLOW_SALT_PLUME */
I , dbloc, Ritop, coriol
I , diffusKzS, diffusKzT
I , ikppkey
O , diffus
U , ghat
O , hbl
I , bi, bj, myTime, myIter, myThid )
c-----------------------------------------------------------------------
c
c Main driver subroutine for kpp vertical mixing scheme and
c interface to greater ocean model
c
c written by: bill large, june 6, 1994
c modified by: jan morzel, june 30, 1994
c bill large, august 11, 1994
c bill large, january 25, 1995 : "dVsq" and 1d code
c detlef stammer, august 1997 : for use with MIT GCM Classic
c d. menemenlis, june 1998 : for use with MIT GCM UV
c
c-----------------------------------------------------------------------
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "KPP_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c input
c bi, bj :: Array indices on which to apply calculations
c myTime :: Current time in simulation
c myIter :: Current iteration number in simulation
c myThid :: My Thread Id. number
c kmtj (imt) - number of vertical layers on this row
c msk (imt) - surface mask (=1 if water, =0 otherwise)
c shsq (imt,Nr) - (local velocity shear)^2 ((m/s)^2)
c dvsq (imt,Nr) - (velocity shear re sfc)^2 ((m/s)^2)
c ustar (imt) - surface friction velocity (m/s)
c bo (imt) - surface turbulent buoy. forcing (m^2/s^3)
c bosol (imt) - radiative buoyancy forcing (m^2/s^3)
c boplume(imt,Nrp1)- haline buoyancy forcing (m^2/s^3)
c dbloc (imt,Nr) - local delta buoyancy across interfaces (m/s^2)
c dblocSm(imt,Nr) - horizontally smoothed dbloc (m/s^2)
c stored in ghat to save space
c Ritop (imt,Nr) - numerator of bulk Richardson Number
c (zref-z) * delta buoyancy w.r.t. surface ((m/s)^2)
c coriol (imt) - Coriolis parameter (1/s)
c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s)
c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s)
c note: there is a conversion from 2-D to 1-D for input output variables,
c e.g., hbl(sNx,sNy) -> hbl(imt),
c where hbl(i,j) -> hbl((j-1)*sNx+i)
INTEGER bi, bj
_RL myTime
integer myIter
integer myThid
integer kmtj (imt )
_RL shsq (imt,Nr)
_RL dvsq (imt,Nr)
_RL ustar (imt )
_RL bo (imt )
_RL bosol (imt )
#ifdef ALLOW_SALT_PLUME
_RL boplume (imt,Nrp1)
_RL SPDepth (imt )
#ifdef SALT_PLUME_SPLIT_BASIN
_RL lon (imt )
_RL lat (imt )
#endif /* SALT_PLUME_SPLIT_BASIN */
#endif /* ALLOW_SALT_PLUME */
_RL dbloc (imt,Nr)
_RL Ritop (imt,Nr)
_RL coriol (imt )
_RS msk (imt )
_RL diffusKzS(imt,Nr)
_RL diffusKzT(imt,Nr)
integer ikppkey
c output
c diffus (imt,1) - vertical viscosity coefficient (m^2/s)
c diffus (imt,2) - vertical scalar diffusivity (m^2/s)
c diffus (imt,3) - vertical temperature diffusivity (m^2/s)
c ghat (imt) - nonlocal transport coefficient (s/m^2)
c hbl (imt) - mixing layer depth (m)
_RL diffus(imt,0:Nrp1,mdiff)
_RL ghat (imt,Nr)
_RL hbl (imt)
#ifdef ALLOW_KPP
c local
c kbl (imt ) - index of first grid level below hbl
c bfsfc (imt ) - surface buoyancy forcing (m^2/s^3)
c casea (imt ) - 1 in case A; 0 in case B
c stable (imt ) - 1 in stable forcing; 0 if unstable
c dkm1 (imt, mdiff) - boundary layer diffusivity at kbl-1 level
c blmc (imt,Nr,mdiff) - boundary layer mixing coefficients
c sigma (imt ) - normalized depth (d / hbl)
c Rib (imt,Nr ) - bulk Richardson number
integer kbl(imt )
_RL bfsfc (imt )
_RL casea (imt )
_RL stable (imt )
_RL dkm1 (imt, mdiff)
_RL blmc (imt,Nr,mdiff)
_RL sigma (imt )
_RL Rib (imt,Nr )
integer i, k, md
c-----------------------------------------------------------------------
c compute interior mixing coefficients everywhere, due to constant
c internal wave activity, static instability, and local shear
c instability.
c (ghat is temporary storage for horizontally smoothed dbloc)
c-----------------------------------------------------------------------
cph(
cph these storings avoid recomp. of Ri_iwmix
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ghat = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
cph)
call RI_IWMIX (
I kmtj, shsq, dbloc, ghat
I , diffusKzS, diffusKzT
I , ikppkey
O , diffus, myThid )
cph(
cph these storings avoid recomp. of Ri_iwmix
cph DESPITE TAFs 'not necessary' warning!
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE shsq = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE ghat = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE diffus = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
cph)
c-----------------------------------------------------------------------
c set seafloor values to zero and fill extra "Nrp1" coefficients
c for blmix
c-----------------------------------------------------------------------
do md = 1, mdiff
do k=1,Nrp1
do i = 1,imt
if(k.ge.kmtj(i)) diffus(i,k,md) = 0.0
end
do
end
do
end
do
c-----------------------------------------------------------------------
c compute boundary layer mixing coefficients:
c
c diagnose the new boundary layer depth
c-----------------------------------------------------------------------
call BLDEPTH (
I kmtj
I , dvsq, dbloc, Ritop, ustar, bo, bosol
#ifdef ALLOW_SALT_PLUME
I , boplume,SPDepth
#ifdef SALT_PLUME_SPLIT_BASIN
I , lon,lat
#endif /* SALT_PLUME_SPLIT_BASIN */
#endif /* ALLOW_SALT_PLUME */
I , coriol
I , ikppkey
O , hbl, bfsfc, stable, casea, kbl, Rib, sigma
I , bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
CADJ & key=ikppkey, kind=isbyte
#endif
c-----------------------------------------------------------------------
c compute boundary layer diffusivities
c-----------------------------------------------------------------------
call BLMIX (
I ustar, bfsfc, hbl, stable, casea, diffus, kbl
O , dkm1, blmc, ghat, sigma, ikppkey
I , myThid )
cph(
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
CADJ & key=ikppkey, kind=isbyte
CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
CADJ & key=ikppkey, kind=isbyte
#endif
cph)
c-----------------------------------------------------------------------
c enhance diffusivity at interface kbl - 1
c-----------------------------------------------------------------------
call ENHANCE (
I dkm1, hbl, kbl, diffus, casea
U , ghat
O , blmc
I , myThid )
cph(
cph avoids recomp. of enhance
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE blmc = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
cph)
c-----------------------------------------------------------------------
c combine interior and boundary layer coefficients and nonlocal term
c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
c (< 1 level), diffusivity blmc can become negative. The max-s below
c are a hack until this problem is properly diagnosed and fixed.
c-----------------------------------------------------------------------
do k = 1, Nr
do i = 1, imt
if (k .lt. kbl(i)) then
#ifdef ALLOW_SHELFICE
C when there is shelfice on top (msk(i)=0), reset the boundary layer
C mixing coefficients blmc to pure Ri-number based mixing
blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
& diffus(i,k,1) )
blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
& diffus(i,k,2) )
blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
& diffus(i,k,3) )
#endif /* not ALLOW_SHELFICE */
diffus(i,k,1) = max ( blmc(i,k,1), viscArNr(1) )
diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
else
ghat(i,k) = 0. _d 0
endif
end
do
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine BLDEPTH (
I kmtj
I , dvsq, dbloc, Ritop, ustar, bo, bosol
#ifdef ALLOW_SALT_PLUME
I , boplume,SPDepth
#ifdef SALT_PLUME_SPLIT_BASIN
I , lon,lat
#endif /* SALT_PLUME_SPLIT_BASIN */
#endif /* ALLOW_SALT_PLUME */
I , coriol
I , ikppkey
O , hbl, bfsfc, stable, casea, kbl, Rib, sigma
I , bi, bj, myTime, myIter, myThid )
c the oceanic planetary boundary layer depth, hbl, is determined as
c the shallowest depth where the bulk Richardson number is
c equal to the critical value, Ricr.
c
c bulk Richardson numbers are evaluated by computing velocity and
c buoyancy differences between values at zgrid(kl) < 0 and surface
c reference values.
c in this configuration, the reference values are equal to the
c values in the surface layer.
c when using a very fine vertical grid, these values should be
c computed as the vertical average of velocity and buoyancy from
c the surface down to epsilon*zgrid(kl).
c
c when the bulk Richardson number at k exceeds Ricr, hbl is
c linearly interpolated between grid levels zgrid(k) and zgrid(k-1).
c
c The water column and the surface forcing are diagnosed for
c stable/ustable forcing conditions, and where hbl is relative
c to grid points (caseA), so that conditional branches can be
c avoided in later subroutines.
c
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "KPP_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c input
c------
c bi, bj :: Array indices on which to apply calculations
c myTime :: Current time in simulation
c myIter :: Current iteration number in simulation
c myThid :: My Thread Id. number
c kmtj : number of vertical layers
c dvsq : (velocity shear re sfc)^2 ((m/s)^2)
c dbloc : local delta buoyancy across interfaces (m/s^2)
c Ritop : numerator of bulk Richardson Number
c =(z-zref)*dbsfc, where dbsfc=delta
c buoyancy with respect to surface ((m/s)^2)
c ustar : surface friction velocity (m/s)
c bo : surface turbulent buoyancy forcing (m^2/s^3)
c bosol : radiative buoyancy forcing (m^2/s^3)
c boplume : haline buoyancy forcing (m^2/s^3)
c coriol : Coriolis parameter (1/s)
INTEGER bi, bj
_RL myTime
integer myIter
integer myThid
integer kmtj(imt)
_RL dvsq (imt,Nr)
_RL dbloc (imt,Nr)
_RL Ritop (imt,Nr)
_RL ustar (imt)
_RL bo (imt)
_RL bosol (imt)
_RL coriol (imt)
integer ikppkey
#ifdef ALLOW_SALT_PLUME
_RL boplume (imt,Nrp1)
_RL SPDepth (imt)
#ifdef SALT_PLUME_SPLIT_BASIN
_RL lon (imt)
_RL lat (imt)
#endif /* SALT_PLUME_SPLIT_BASIN */
#endif /* ALLOW_SALT_PLUME */
c output
c--------
c hbl : boundary layer depth (m)
c bfsfc : Bo+radiation absorbed to d=hbf*hbl (m^2/s^3)
c stable : =1 in stable forcing; =0 unstable
c casea : =1 in case A, =0 in case B
c kbl : -1 of first grid level below hbl
c Rib : Bulk Richardson number
c sigma : normalized depth (d/hbl)
_RL hbl (imt)
_RL bfsfc (imt)
_RL stable (imt)
_RL casea (imt)
integer kbl(imt)
_RL Rib (imt,Nr)
_RL sigma (imt)
#ifdef ALLOW_KPP
c local
c-------
c wm, ws : turbulent velocity scales (m/s)
_RL wm(imt), ws(imt)
_RL worka(imt)
_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
integer i, kl
_RL p5 , eins
parameter ( p5=0.5, eins=1.0 )
_RL minusone
parameter ( minusone=-1.0 )
#ifdef SALT_PLUME_VOLUME
integer km, km1
_RL temp
#endif
#ifdef ALLOW_AUTODIFF_TAMC
integer kkppkey
#endif
#ifdef ALLOW_DIAGNOSTICS
c KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
_RL KPPBFSFC(imt,Nr)
_RL KPPRi(imt,Nr)
#endif /* ALLOW_DIAGNOSTICS */
c find bulk Richardson number at every grid level until > Ricr
c
c note: the reference depth is -epsilon/2.*zgrid(k), but the reference
c u,v,t,s values are simply the surface layer values,
c and not the averaged values from 0 to 2*ref.depth,
c which is necessary for very fine grids(top layer < 2m thickness)
c note: max values when Ricr never satisfied are
c kbl(i)=kmtj(i) and hbl(i)=-zgrid(kmtj(i))
c initialize hbl and kbl to bottomed out values
do i = 1, imt
Rib(i,1) = 0. _d 0
if (kmtj(i).LT.1) kmtj(i) = 1
kbl(i) = kmtj(i)
kl = kbl(i)
hbl(i) = -zgrid(kl)
end
do
#ifdef ALLOW_DIAGNOSTICS
do kl = 1, Nr
do i = 1, imt
KPPBFSFC(i,kl) = 0. _d 0
KPPRi(i,kl) = 0. _d 0
enddo
enddo
#endif /* ALLOW_DIAGNOSTICS */
do kl = 2, Nr
#ifdef ALLOW_AUTODIFF_TAMC
kkppkey = (ikppkey-1)*Nr + kl
#endif
c compute bfsfc = sw fraction at hbf * zgrid
do i = 1, imt
worka(i) = zgrid(kl)
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
#endif
call SWFRAC(
I imt, hbf,
U worka,
I myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
#endif
do i = 1, imt
c use caseA as temporary array
casea(i) = -zgrid(kl)
c compute bfsfc= Bo + radiative contribution down to hbf * hbl
bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
end
do
#ifdef ALLOW_SALT_PLUME
c compute bfsfc = plume fraction at hbf * zgrid
IF ( useSALT_PLUME ) THEN
#ifndef SALT_PLUME_VOLUME
do i = 1, imt
worka(i) = zgrid(kl)
enddo
Ccatn: in original way: accumulate all fractions of boplume above zgrid(kl)
call SALT_PLUME_FRAC(
I imt, hbf,SPDepth,
#ifdef SALT_PLUME_SPLIT_BASIN
I lon,lat,
#endif /* SALT_PLUME_SPLIT_BASIN */
U worka,
I myTime, myIter, myThid)
do i = 1, imt
bfsfc(i) = bfsfc(i) + boplume(i,1)*(worka(i))
C km=max(1,kbl(i)-1)
C temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
C bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
enddo
#else /* def SALT_PLUME_VOLUME */
catn: in vol way: need to integrate down to hbl, so first locate
c k level associated with this hbl, then sum up all SPforc[T,S]
DO i = 1, imt
km =max(1,kbl(i)-1)
km1=max(1,kbl(i))
temp = (boplume(i,km)+boplume(i,km1))*p5
bfsfc(i) = bfsfc(i) + temp
ENDDO
#endif /* ndef SALT_PLUME_VOLUME */
ENDIF
#endif /* ALLOW_SALT_PLUME */
#ifdef ALLOW_DIAGNOSTICS
do i = 1, imt
KPPBFSFC(i,kl) = bfsfc(i)
enddo
#endif /* ALLOW_DIAGNOSTICS */
do i = 1, imt
stable(i) = p5 + sign(p5,bfsfc(i))
sigma(i) = stable(i) + (1. - stable(i)) * epsilon
enddo
c-----------------------------------------------------------------------
c compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
c-----------------------------------------------------------------------
call WSCALE (
I sigma, casea, ustar, bfsfc,
O wm, ws, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store ws = comlev1_kpp_k, key = kkppkey, kind=isbyte
#endif
do i = 1, imt
c-----------------------------------------------------------------------
c compute the turbulent shear contribution to Rib
c-----------------------------------------------------------------------
bvsq = p5 *
1 ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl ))+
2 dbloc(i,kl ) / (zgrid(kl )-zgrid(kl+1)))
if (bvsq .eq. 0. _d 0) then
vtsq = 0. _d 0
else
vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
endif
c compute bulk Richardson number at new level
c note: Ritop needs to be zero on land and ocean bottom
c points so that the following if statement gets triggered
c correctly; otherwise, hbl might get set to (big) negative
c values, that might exceed the limit for the "exp" function
c in "SWFRAC"
c
c rg: assignment to double precision variable to avoid overflow
c ph: test for zero nominator
c
tempVar1 = dvsq(i,kl) + vtsq
tempVar2 = max(tempVar1, phepsi)
Rib(i,kl) = Ritop(i,kl) / tempVar2
#ifdef ALLOW_DIAGNOSTICS
KPPRi(i,kl) = Rib(i,kl)
#endif
end
do
end
do
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(KPPRi ,'KPPRi ',0,Nr,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
cph(
cph without this store, there is a recomputation error for
cph rib in adbldepth (probably partial recomputation problem)
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store Rib = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
#endif
cph)
do kl = 2, Nr
do i = 1, imt
if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl
end
do
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store kbl = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
do i = 1, imt
kl = kbl(i)
c linearly interpolate to find hbl where Rib = Ricr
if (kl.gt.1 .and. kl.lt.kmtj(i)) then
tempVar1 = (Rib(i,kl)-Rib(i,kl-1))
hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
1 (Ricr - Rib(i,kl-1)) / tempVar1
endif
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store hbl = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
c-----------------------------------------------------------------------
c find stability and buoyancy forcing for boundary layer
c-----------------------------------------------------------------------
do i = 1, imt
worka(i) = hbl(i)
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
call SWFRAC(
I imt, minusone,
U worka,
I myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
do i = 1, imt
bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
end
do
#ifdef ALLOW_SALT_PLUME
IF ( useSALT_PLUME ) THEN
#ifndef SALT_PLUME_VOLUME
do i = 1, imt
worka(i) = hbl(i)
enddo
call SALT_PLUME_FRAC(
I imt,minusone,SPDepth,
#ifdef SALT_PLUME_SPLIT_BASIN
I lon,lat,
#endif /* SALT_PLUME_SPLIT_BASIN */
U worka,
I myTime, myIter, myThid )
do i = 1, imt
bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
C km=max(1,kbl(i)-1)
C temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
C bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
enddo
#else /* def SALT_PLUME_VOLUME */
DO i = 1, imt
km =max(1,kbl(i)-1)
km1=max(1,kbl(i))
temp = (boplume(i,km)+boplume(i,km1))/2.0
bfsfc(i) = bfsfc(i) + temp
ENDDO
#endif /* ndef SALT_PLUME_VOLUME */
ENDIF
#endif /* ALLOW_SALT_PLUME */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bfsfc = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
c-- ensure bfsfc is never 0
do i = 1, imt
stable(i) = p5 + sign( p5, bfsfc(i) )
bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
end
do
cph(
cph added stable to store list to avoid extensive recomp.
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bfsfc, stable = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
cph)
c-----------------------------------------------------------------------
c check hbl limits for hekman or hmonob
c ph: test for zero nominator
c-----------------------------------------------------------------------
IF ( LimitHblStable ) THEN
do i = 1, imt
if (bfsfc(i) .gt. 0.0) then
hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)
hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
& / vonk / bfsfc(i)
hlimit = stable(i) * min(hekman,hmonob)
& + (stable(i)-1.) * zgrid(Nr)
hbl(i) = min(hbl(i),hlimit)
end
if
end
do
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store hbl = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
do i = 1, imt
hbl(i) = max(hbl(i),minKPPhbl)
kbl(i) = kmtj(i)
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store hbl = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
c-----------------------------------------------------------------------
c find new kbl
c-----------------------------------------------------------------------
do kl = 2, Nr
do i = 1, imt
if ( kbl(i).eq.kmtj(i) .and. (-zgrid(kl)).gt.hbl(i) ) then
kbl(i) = kl
endif
end
do
end
do
c-----------------------------------------------------------------------
c find stability and buoyancy forcing for final hbl values
c-----------------------------------------------------------------------
do i = 1, imt
worka(i) = hbl(i)
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
call SWFRAC(
I imt, minusone,
U worka,
I myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store worka = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
do i = 1, imt
bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
end
do
#ifdef ALLOW_SALT_PLUME
IF ( useSALT_PLUME ) THEN
#ifndef SALT_PLUME_VOLUME
do i = 1, imt
worka(i) = hbl(i)
enddo
call SALT_PLUME_FRAC(
I imt,minusone,SPDepth,
#ifdef SALT_PLUME_SPLIT_BASIN
I lon,lat,
#endif /* SALT_PLUME_SPLIT_BASIN */
U worka,
I myTime, myIter, myThid )
do i = 1, imt
bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
C km=max(1,kbl(i)-1)
C temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
C bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
enddo
#else /* def SALT_PLUME_VOLUME */
DO i = 1, imt
km =max(1,kbl(i)-1)
km1=max(1,kbl(i)-0)
temp = (boplume(i,km)+boplume(i,km1))/2.0
bfsfc(i) = bfsfc(i) + temp
ENDDO
#endif /* ndef SALT_PLUME_VOLUME */
ENDIF
#endif /* ALLOW_SALT_PLUME */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bfsfc = comlev1_kpp
CADJ & , key=ikppkey, kind=isbyte,
CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
#endif
c-- ensures bfsfc is never 0
do i = 1, imt
stable(i) = p5 + sign( p5, bfsfc(i) )
bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
end
do
c-----------------------------------------------------------------------
c determine caseA and caseB
c-----------------------------------------------------------------------
do i = 1, imt
kl = kbl(i)
casea(i) = p5 +
1 sign(p5, -zgrid(kl) - p5*hwide(kl) - hbl(i))
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine WSCALE (
I sigma, hbl, ustar, bfsfc,
O wm, ws,
I myThid )
c compute turbulent velocity scales.
c use a 2D-lookup table for wm and ws as functions of ustar and
c zetahat (=vonk*sigma*hbl*bfsfc).
c
c note: the lookup table is only used for unstable conditions
c (zehat.le.0), in the stable domain wm (=ws) gets computed
c directly.
c
IMPLICIT NONE
#include "SIZE.h"
#include "KPP_PARAMS.h"
c input
c------
c sigma : normalized depth (d/hbl)
c hbl : boundary layer depth (m)
c ustar : surface friction velocity (m/s)
c bfsfc : total surface buoyancy flux (m^2/s^3)
c myThid : thread number for this instance of the routine
integer myThid
_RL sigma(imt)
_RL hbl (imt)
_RL ustar(imt)
_RL bfsfc(imt)
c output
c--------
c wm, ws : turbulent velocity scales at sigma
_RL wm(imt), ws(imt)
#ifdef ALLOW_KPP
c local
c------
c zehat : = zeta * ustar**3
_RL zehat
integer iz, izp1, ju, i, jup1
_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
_RL wbm, was, wbs, u3, tempVar
c-----------------------------------------------------------------------
c use lookup table for zehat < zmax only; otherwise use
c stable formulae
c-----------------------------------------------------------------------
do i = 1, imt
zehat = vonk*sigma(i)*hbl(i)*bfsfc(i)
if (zehat .le. zmax) then
zdiff = zehat - zmin
C For extremely negative buoyancy forcing bfsfc, zehat and hence
C zdiff can become very negative (default value of zmin = 4.e-7) and
C the extrapolation beyond the limit zmin of the lookup table can
C give very bad values and may make the model crash. Here is a
C simple fix (thanks to Dimitry Sidorenko) that effectively replaces
C linear extrapolation with nearest neighbor extrapolation so that
C only the lower limit values of the lookup tables wmt/wst are used.
C Alternatively, one can get rid of the lookup table altogether
C and compute the coefficients online (done in NEMO, for example).
C zdiff = max( 0. _d 0, zehat - zmin )
iz = int( zdiff / deltaz )
iz = min( iz, nni )
iz = max( iz, 0 )
izp1 = iz + 1
udiff = ustar(i) - umin
ju = int( udiff / deltau )
ju = min( ju, nnj )
ju = max( ju, 0 )
jup1 = ju + 1
zfrac = zdiff / deltaz - float(iz)
ufrac = udiff / deltau - float(ju)
fzfrac= 1. - zfrac
wam = fzfrac * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
wbm = fzfrac * wmt(iz,ju ) + zfrac * wmt(izp1,ju )
wm(i) = (1.-ufrac) * wbm + ufrac * wam
was = fzfrac * wst(iz,jup1) + zfrac * wst(izp1,jup1)
wbs = fzfrac * wst(iz,ju ) + zfrac * wst(izp1,ju )
ws(i) = (1.-ufrac) * wbs + ufrac * was
else
u3 = ustar(i) * ustar(i) * ustar(i)
tempVar = u3 + conc1 * zehat
wm(i) = vonk * ustar(i) * u3 / tempVar
ws(i) = wm(i)
endif
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine RI_IWMIX (
I kmtj, shsq, dbloc, dblocSm,
I diffusKzS, diffusKzT,
I ikppkey,
O diffus,
I myThid )
c compute interior viscosity diffusivity coefficients due
c to shear instability (dependent on a local Richardson number),
c to background internal wave activity, and
c to static instability (local Richardson number < 0).
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "KPP_PARAMS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_PARAMS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c input
c kmtj (imt) number of vertical layers on this row
c shsq (imt,Nr) (local velocity shear)^2 ((m/s)^2)
c dbloc (imt,Nr) local delta buoyancy (m/s^2)
c dblocSm(imt,Nr) horizontally smoothed dbloc (m/s^2)
c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s)
c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s)
c myThid :: My Thread Id. number
integer kmtj (imt)
_RL shsq (imt,Nr)
_RL dbloc (imt,Nr)
_RL dblocSm (imt,Nr)
_RL diffusKzS(imt,Nr)
_RL diffusKzT(imt,Nr)
integer ikppkey
integer myThid
c output
c diffus(imt,0:Nrp1,1) vertical viscosivity coefficient (m^2/s)
c diffus(imt,0:Nrp1,2) vertical scalar diffusivity (m^2/s)
c diffus(imt,0:Nrp1,3) vertical temperature diffusivity (m^2/s)
_RL diffus(imt,0:Nrp1,3)
#ifdef ALLOW_KPP
c local variables
c Rig local Richardson number
c fRi, fcon function of Rig
_RL Rig
_RL fRi, fcon
_RL ratio
integer i, ki, kp1
_RL c1, c0
#ifdef ALLOW_KPP_VERTICALLY_SMOOTH
integer mr
CADJ INIT kpp_ri_tape_mr = common, 1
#endif
c constants
c1 = 1. _d 0
c0 = 0. _d 0
c-----------------------------------------------------------------------
c compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)
c use diffus(*,*,1) as temporary storage of Ri to be smoothed
c use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared
c set values at bottom and below to nearest value above bottom
#ifdef ALLOW_AUTODIFF
C break data flow dependence on diffus
diffus(1,1,1) = 0.0
do ki = 1, Nr
do i = 1, imt
diffus(i,ki,1) = 0.
diffus(i,ki,2) = 0.
diffus(i,ki,3) = 0.
enddo
enddo
#endif
do ki = 1, Nr
do i = 1, imt
if (kmtj(i) .LE. 1 ) then
diffus(i,ki,1) = 0.
diffus(i,ki,2) = 0.
elseif (ki .GE. kmtj(i)) then
diffus(i,ki,1) = diffus(i,ki-1,1)
diffus(i,ki,2) = diffus(i,ki-1,2)
else
diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))
& / max( Shsq(i,ki), phepsi )
diffus(i,ki,2) = dbloc(i,ki) / (zgrid(ki)-zgrid(ki+1))
endif
end
do
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store diffus = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
c-----------------------------------------------------------------------
c vertically smooth Ri
#ifdef ALLOW_KPP_VERTICALLY_SMOOTH
do mr = 1, num_v_smooth_Ri
#ifdef ALLOW_AUTODIFF_TAMC
CADJ store diffus(:,:,1) = kpp_ri_tape_mr
CADJ & , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
#endif
call Z121 (
U diffus(1,0,1),
I myThid )
end
do
#endif
c-----------------------------------------------------------------------
c after smoothing loop
do ki = 1, Nr
do i = 1, imt
c evaluate f of Brunt-Vaisala squared for convection, store in fcon
Rig = max ( diffus(i,ki,2) , BVSQcon )
ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )
fcon = c1 - ratio * ratio
fcon = fcon * fcon * fcon
c evaluate f of smooth Ri for shear instability, store in fRi
Rig = max ( diffus(i,ki,1), c0 )
ratio = min ( Rig / Riinfty , c1 )
fRi = c1 - ratio * ratio
fRi = fRi * fRi * fRi
c ----------------------------------------------------------------------
c evaluate diffusivities and viscosity
c mixing due to internal waves, and shear and static instability
kp1 = MIN(ki+1,Nr)
#ifdef EXCLUDE_KPP_SHEAR_MIX
diffus(i,ki,1) = viscArNr(1)
diffus(i,ki,2) = diffusKzS(i,kp1)
diffus(i,ki,3) = diffusKzT(i,kp1)
#else /* EXCLUDE_KPP_SHEAR_MIX */
# ifdef ALLOW_AUTODIFF
if ( inAdMode ) then
diffus(i,ki,1) = viscArNr(1)
diffus(i,ki,2) = diffusKzS(i,kp1)
diffus(i,ki,3) = diffusKzT(i,kp1)
else
else /* ALLOW_AUTODIFF */
if ( .TRUE. ) then
# endif /* ALLOW_AUTODIFF */
diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0
diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0
diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0
endif
#endif /* EXCLUDE_KPP_SHEAR_MIX */
end
do
end
do
c ------------------------------------------------------------------------
c set surface values to 0.0
do i = 1, imt
diffus(i,0,1) = c0
diffus(i,0,2) = c0
diffus(i,0,3) = c0
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine Z121 (
U v,
I myThid )
c Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
c top (0) value is used as a dummy
c bottom (Nrp1) value is set to input value from above.
c Note that it is important to exclude from the smoothing any points
c that are outside the range of the K(Ri) scheme, ie. >0.8, or <0.0.
c Otherwise, there is interference with other physics, especially
c penetrative convection.
IMPLICIT NONE
#include "SIZE.h"
#include "KPP_PARAMS.h"
c input/output
c-------------
c v : 2-D array to be smoothed in Nrp1 direction
c myThid: thread number for this instance of the routine
integer myThid
_RL v(imt,0:Nrp1)
#ifdef ALLOW_KPP
c local
_RL zwork, zflag
_RL KRi_range(1:Nrp1)
integer i, k, km1, kp1
_RL p0 , p25 , p5 , p2
parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
KRi_range(Nrp1) = p0
#ifdef ALLOW_AUTODIFF_TAMC
C-- dummy assignment to end declaration part for TAMC
i = 0
C-- HPF directive to help TAMC
CHPF$ INDEPENDENT
CADJ INIT z121tape = common, Nr
#endif /* ALLOW_AUTODIFF_TAMC */
do i = 1, imt
k = 1
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE v(i,k) = z121tape
#endif
v(i,Nrp1) = v(i,Nr)
do k = 1, Nr
KRi_range(k) = p5 + SIGN(p5,v(i,k))
KRi_range(k) = KRi_range(k) *
& ( p5 + SIGN(p5,(Riinfty-v(i,k))) )
end
do
zwork = KRi_range(1) * v(i,1)
v(i,1) = p2 * v(i,1) +
& KRi_range(1) * KRi_range(2) * v(i,2)
zflag = p2 + KRi_range(1) * KRi_range(2)
v(i,1) = v(i,1) / zflag
do k = 2, Nr
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE v(i,k), zwork = z121tape
#endif
km1 = k - 1
kp1 = k + 1
zflag = v(i,k)
v(i,k) = p2 * v(i,k) +
& KRi_range(k) * KRi_range(kp1) * v(i,kp1) +
& KRi_range(k) * zwork
zwork = KRi_range(k) * zflag
zflag = p2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1))
v(i,k) = v(i,k) / zflag
end
do
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine SMOOTH_HORIZ (
I k, bi, bj,
U fld,
I myThid )
c Apply horizontal smoothing to global _RL 2-D array
IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "KPP_PARAMS.h"
c input
c bi, bj : array indices
c k : vertical index used for masking
c myThid : thread number for this instance of the routine
INTEGER myThid
integer k, bi, bj
c input/output
c fld : 2-D array to be smoothed
_RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
#ifdef ALLOW_KPP
c local
integer i, j, im1, ip1, jm1, jp1
_RL tempVar
_RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
integer imin , imax , jmin , jmax
parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
_RL p0 , p5 , p25 , p125 , p0625
parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
DO j = jmin, jmax
jm1 = j-1
jp1 = j+1
DO i = imin, imax
im1 = i-1
ip1 = i+1
tempVar =
& p25 * maskC(i ,j ,k,bi,bj) +
& p125 * ( maskC(im1,j ,k,bi,bj) +
& maskC(ip1,j ,k,bi,bj) +
& maskC(i ,jm1,k,bi,bj) +
& maskC(i ,jp1,k,bi,bj) ) +
& p0625 * ( maskC(im1,jm1,k,bi,bj) +
& maskC(im1,jp1,k,bi,bj) +
& maskC(ip1,jm1,k,bi,bj) +
& maskC(ip1,jp1,k,bi,bj) )
IF ( tempVar .GE. p25 ) THEN
fld_tmp(i,j) = (
& p25 * fld(i ,j )*maskC(i ,j ,k,bi,bj) +
& p125 *(fld(im1,j )*maskC(im1,j ,k,bi,bj) +
& fld(ip1,j )*maskC(ip1,j ,k,bi,bj) +
& fld(i ,jm1)*maskC(i ,jm1,k,bi,bj) +
& fld(i ,jp1)*maskC(i ,jp1,k,bi,bj))+
& p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
& fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
& fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
& fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
& / tempVar
ELSE
fld_tmp(i,j) = fld(i,j)
ENDIF
ENDDO
ENDDO
c transfer smoothed field to output array
DO j = jmin, jmax
DO i = imin, imax
fld(i,j) = fld_tmp(i,j)
ENDDO
ENDDO
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine BLMIX (
I ustar, bfsfc, hbl, stable, casea, diffus, kbl
O , dkm1, blmc, ghat, sigma, ikppkey
I , myThid )
c mixing coefficients within boundary layer depend on surface
c forcing and the magnitude and gradient of interior mixing below
c the boundary layer ("matching").
c
c caution: if mixing bottoms out at hbl = -zgrid(Nr) then
c fictitious layer at Nrp1 is needed with small but finite width
c hwide(Nrp1) (eg. epsln = 1.e-20).
c
IMPLICIT NONE
#include "SIZE.h"
#include "KPP_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c input
c ustar (imt) surface friction velocity (m/s)
c bfsfc (imt) surface buoyancy forcing (m^2/s^3)
c hbl (imt) boundary layer depth (m)
c stable(imt) = 1 in stable forcing
c casea (imt) = 1 in case A
c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s)
c kbl (imt) -1 of first grid level below hbl
c myThid thread number for this instance of the routine
integer myThid
_RL ustar (imt)
_RL bfsfc (imt)
_RL hbl (imt)
_RL stable(imt)
_RL casea (imt)
_RL diffus(imt,0:Nrp1,mdiff)
integer kbl(imt)
c output
c dkm1 (imt,mdiff) boundary layer difs at kbl-1 level
c blmc (imt,Nr,mdiff) boundary layer mixing coefficients (m^2/s)
c ghat (imt,Nr) nonlocal scalar transport
c sigma(imt) normalized depth (d / hbl)
_RL dkm1 (imt,mdiff)
_RL blmc (imt,Nr,mdiff)
_RL ghat (imt,Nr)
_RL sigma(imt)
integer ikppkey
#ifdef ALLOW_KPP
c local
c gat1*(imt) shape function at sigma = 1
c dat1*(imt) derivative of shape function at sigma = 1
c ws(imt), wm(imt) turbulent velocity scales (m/s)
_RL gat1m(imt), gat1s(imt), gat1t(imt)
_RL dat1m(imt), dat1s(imt), dat1t(imt)
_RL ws(imt), wm(imt)
integer i, kn, ki, kl
_RL R, dvdzup, dvdzdn, viscp
_RL difsp, diftp, visch, difsh, difth
_RL f1, sig, a1, a2, a3, delhat
_RL Gm, Gs, Gt
_RL tempVar
_RL p0 , eins
parameter (p0=0.0, eins=1.0)
#ifdef ALLOW_AUTODIFF_TAMC
integer kkppkey
#endif
c-----------------------------------------------------------------------
c compute velocity scales at hbl
c-----------------------------------------------------------------------
do i = 1, imt
sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
call WSCALE (
I sigma, hbl, ustar, bfsfc,
O wm, ws, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
do i = 1, imt
wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
end
do
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
do i = 1, imt
kn = int(caseA(i)+phepsi) *(kbl(i) -1) +
$ (1 - int(caseA(i)+phepsi)) * kbl(i)
c-----------------------------------------------------------------------
c find the interior viscosities and derivatives at hbl(i)
c-----------------------------------------------------------------------
delhat = 0.5*hwide(kn) - zgrid(kn) - hbl(i)
R = 1.0 - delhat / hwide(kn)
dvdzup = (diffus(i,kn-1,1) - diffus(i,kn ,1)) / hwide(kn)
dvdzdn = (diffus(i,kn ,1) - diffus(i,kn+1,1)) / hwide(kn+1)
viscp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1 R * (dvdzdn + abs(dvdzdn)) )
dvdzup = (diffus(i,kn-1,2) - diffus(i,kn ,2)) / hwide(kn)
dvdzdn = (diffus(i,kn ,2) - diffus(i,kn+1,2)) / hwide(kn+1)
difsp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1 R * (dvdzdn + abs(dvdzdn)) )
dvdzup = (diffus(i,kn-1,3) - diffus(i,kn ,3)) / hwide(kn)
dvdzdn = (diffus(i,kn ,3) - diffus(i,kn+1,3)) / hwide(kn+1)
diftp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1 R * (dvdzdn + abs(dvdzdn)) )
visch = diffus(i,kn,1) + viscp * delhat
difsh = diffus(i,kn,2) + difsp * delhat
difth = diffus(i,kn,3) + diftp * delhat
f1 = stable(i) * conc1 * bfsfc(i) /
& max(ustar(i)**4,phepsi)
gat1m(i) = visch / hbl(i) / wm(i)
dat1m(i) = -viscp / wm(i) + f1 * visch
gat1s(i) = difsh / hbl(i) / ws(i)
dat1s(i) = -difsp / ws(i) + f1 * difsh
gat1t(i) = difth / hbl(i) / ws(i)
dat1t(i) = -diftp / ws(i) + f1 * difth
end
do
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE gat1m = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE gat1s = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE gat1t = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
do i = 1, imt
dat1m(i) = min(dat1m(i),p0)
dat1s(i) = min(dat1s(i),p0)
dat1t(i) = min(dat1t(i),p0)
end
do
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
do ki = 1, Nr
#ifdef ALLOW_AUTODIFF_TAMC
kkppkey = (ikppkey-1)*Nr + ki
#endif
c-----------------------------------------------------------------------
c compute turbulent velocity scales on the interfaces
c-----------------------------------------------------------------------
do i = 1, imt
sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
end
do
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE wm = comlev1_kpp_k, key = kkppkey
CADJ STORE ws = comlev1_kpp_k, key = kkppkey
#endif
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
#endif
call WSCALE (
I sigma, hbl, ustar, bfsfc,
O wm, ws, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wm = comlev1_kpp_k, key = kkppkey
CADJ STORE ws = comlev1_kpp_k, key = kkppkey
#endif
c-----------------------------------------------------------------------
c compute the dimensionless shape functions at the interfaces
c-----------------------------------------------------------------------
do i = 1, imt
sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
a1 = sig - 2.
a2 = 3. - 2. * sig
a3 = sig - 1.
Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
c-----------------------------------------------------------------------
c compute boundary layer diffusivities at the interfaces
c-----------------------------------------------------------------------
blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
c-----------------------------------------------------------------------
c nonlocal transport term = ghat * o
c-----------------------------------------------------------------------
tempVar = ws(i) * hbl(i)
ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)
end
do
end
do
c-----------------------------------------------------------------------
c find diffusivities at kbl-1 grid level
c-----------------------------------------------------------------------
do i = 1, imt
kl = kbl(i)
sig = -zgrid(kl-1) / hbl(i)
sigma(i) = stable(i) * sig
& + (1. - stable(i)) * min(sig,epsilon)
end
do
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
call WSCALE (
I sigma, hbl, ustar, bfsfc,
O wm, ws, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
#endif
do i = 1, imt
kl = kbl(i)
sig = -zgrid(kl-1) / hbl(i)
a1 = sig - 2.
a2 = 3. - 2. * sig
a3 = sig - 1.
Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
dkm1(i,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
dkm1(i,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
dkm1(i,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
subroutine ENHANCE (
I dkm1, hbl, kbl, diffus, casea
U , ghat
O , blmc
& , myThid )
c enhance the diffusivity at the kbl-.5 interface
IMPLICIT NONE
#include "SIZE.h"
#include "KPP_PARAMS.h"
c input
c dkm1(imt,mdiff) bl diffusivity at kbl-1 grid level
c hbl(imt) boundary layer depth (m)
c kbl(imt) grid above hbl
c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s)
c casea(imt) = 1 in caseA, = 0 in case B
c myThid thread number for this instance of the routine
integer myThid
_RL dkm1 (imt,mdiff)
_RL hbl (imt)
integer kbl (imt)
_RL diffus(imt,0:Nrp1,mdiff)
_RL casea (imt)
c input/output
c nonlocal transport, modified ghat at kbl(i)-1 interface (s/m**2)
_RL ghat (imt,Nr)
c output
c enhanced bound. layer mixing coeff.
_RL blmc (imt,Nr,mdiff)
#ifdef ALLOW_KPP
c local
c fraction hbl lies beteen zgrid neighbors
_RL delta
integer ki, i, md
_RL dkmp5, dstar
do i = 1, imt
ki = kbl(i)-1
if ((ki .ge. 1) .and. (ki .lt. Nr)) then
delta = (hbl(i) + zgrid(ki)) / (zgrid(ki) - zgrid(ki+1))
do md = 1, mdiff
dkmp5 = casea(i) * diffus(i,ki,md) +
1 (1.- casea(i)) * blmc (i,ki,md)
dstar = (1.- delta)**2 * dkm1(i,md)
& + delta**2 * dkmp5
blmc(i,ki,md) = (1.- delta)*diffus(i,ki,md)
& + delta*dstar
end
do
ghat(i,ki) = (1.- casea(i)) * ghat(i,ki)
endif
end
do
#endif /* ALLOW_KPP */
return
end
c*************************************************************************
SUBROUTINE STATEKPP (
O RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
I ikppkey, bi, bj, myThid )
c
c-----------------------------------------------------------------------
c "statekpp" computes all necessary input arrays
c for the kpp mixing scheme
c
c input:
c bi, bj = array indices on which to apply calculations
c
c output:
c rho1 = potential density of surface layer (kg/m^3)
c dbloc = local buoyancy gradient at Nr interfaces
c g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ] (m/s^2)
c dbsfc = buoyancy difference with respect to the surface
c g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ] (m/s^2)
c ttalpha= thermal expansion coefficient without 1/rho factor
c d(rho) / d(potential temperature) (kg/m^3/C)
c ssbeta = salt expansion coefficient without 1/rho factor
c d(rho) / d(salinity) (kg/m^3/PSU)
c
c see also subroutines find_rho.F find_alpha.F find_beta.F
c
c written by: jan morzel, feb. 10, 1995 (converted from "sigma" version)
c modified by: d. menemenlis, june 1998 : for use with MIT GCM UV
c
c-----------------------------------------------------------------------
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "KPP_PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c-------------- Routine arguments -----------------------------------------
INTEGER bi, bj, myThid
_RL RHO1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
_RL DBLOC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
_RL DBSFC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
_RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
_RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
#ifdef ALLOW_KPP
c--------------------------------------------------------------------------
c
c local arrays:
c
c rhok - density of t(k ) & s(k ) at depth k
c rhokm1 - density of t(k-1) & s(k-1) at depth k
c rho1k - density of t(1 ) & s(1 ) at depth k
c work1,2,3 - work arrays for holding horizontal slabs
_RL RHOK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL RHO1K (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER I, J, K
INTEGER ikppkey, kkppkey
c calculate density, alpha, beta in surface layer, and set dbsfc to zero
kkppkey = (ikppkey-1)*Nr + 1
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
CALL FIND_RHO_2D(
I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
I theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
O WORK1,
I 1, bi, bj, myThid )
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
call FIND_ALPHA(
I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
O WORK2, myThid )
call FIND_BETA(
I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
O WORK3, myThid )
DO J = 1-OLy, sNy+OLy
DO I = 1-OLx, sNx+OLx
RHO1(I,J) = WORK1(I,J) + rhoConst
TTALPHA(I,J,1) = WORK2(I,J)
SSBETA(I,J,1) = WORK3(I,J)
DBSFC(I,J,1) = 0.
END
DO
END
DO
c calculate alpha, beta, and gradients in interior layers
CHPF$ INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
DO K = 2, Nr
kkppkey = (ikppkey-1)*Nr + k
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
CALL FIND_RHO_2D(
I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O RHOK,
I k, bi, bj, myThid )
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
CALL FIND_RHO_2D(
I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
I theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
O RHOKM1,
I k-1, bi, bj, myThid )
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
CALL FIND_RHO_2D(
I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
I theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
O RHO1K,
I 1, bi, bj, myThid )
#ifdef KPP_AUTODIFF_MORE_STORE
CADJ STORE rhok (:,:) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE rhokm1(:,:) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
CADJ STORE rho1k (:,:) = comlev1_kpp_k,
CADJ & key=kkppkey, kind=isbyte
#endif /* KPP_AUTODIFF_MORE_STORE */
call FIND_ALPHA(
I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
O WORK1, myThid )
call FIND_BETA(
I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
O WORK2, myThid )
DO J = 1-OLy, sNy+OLy
DO I = 1-OLx, sNx+OLx
TTALPHA(I,J,K) = WORK1 (I,J)
SSBETA(I,J,K) = WORK2 (I,J)
DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
& (RHOK(I,J) + rhoConst)
DBSFC(I,J,K) = gravity * (RHOK(I,J) - RHO1K (I,J)) /
& (RHOK(I,J) + rhoConst)
END
DO
END
DO
END
DO
c compute arrays for K = Nrp1
DO J = 1-OLy, sNy+OLy
DO I = 1-OLx, sNx+OLx
TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
SSBETA(I,J,Nrp1) = SSBETA(I,J,Nr)
DBLOC(I,J,Nr) = 0.
END
DO
END
DO
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_KPP */
RETURN
END
c*************************************************************************
SUBROUTINE KPP_DOUBLEDIFF (
I TTALPHA, SSBETA,
U kappaRT,
U kappaRS,
I ikppkey, imin, imax, jmin, jmax, bi, bj, myThid )
c
c-----------------------------------------------------------------------
c "KPP_DOUBLEDIFF" adds the double diffusive contributions
C as Rrho-dependent parameterizations to kappaRT and kappaRS
c
c input:
c bi, bj = array indices on which to apply calculations
c imin, imax, jmin, jmax = array boundaries
c ikppkey = key for TAMC/TAF automatic differentiation
c myThid = thread id
c
c ttalpha= thermal expansion coefficient without 1/rho factor
c d(rho) / d(potential temperature) (kg/m^3/C)
c ssbeta = salt expansion coefficient without 1/rho factor
c d(rho) / d(salinity) (kg/m^3/PSU)
c output: updated
c kappaRT/S :: background diffusivities for temperature and salinity
c
c written by: martin losch, sept. 15, 2009
c
c-----------------------------------------------------------------------
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "KPP_PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
c-------------- Routine arguments -----------------------------------------
INTEGER ikppkey, imin, imax, jmin, jmax, bi, bj, myThid
_RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
_RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
_RL KappaRT( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
_RL KappaRS( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
#ifdef ALLOW_KPP
C--------------------------------------------------------------------------
C
C local variables
C I,J,K :: loop indices
C kkppkey :: key for TAMC/TAF automatic differentiation
C
INTEGER I, J, K
INTEGER kkppkey
C alphaDT :: d\rho/d\theta * d\theta
C betaDS :: d\rho/dsalt * dsalt
C Rrho :: "density ratio" R_{\rho} = \alpha dT/dz / \beta dS/dz
C nuddt/s :: double diffusive diffusivities
C numol :: molecular diffusivity
C rFac :: abbreviation for 1/(R_{\rho0}-1)
_RL alphaDT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
_RL betaDS ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
_RL nuddt ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
_RL nudds ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
_RL Rrho
_RL numol, rFac, nutmp
INTEGER Km1
C set some constants here
numol = 1.5 _d -06
rFac = 1. _d 0 / (Rrho0 - 1. _d 0 )
C
kkppkey = (ikppkey-1)*Nr + 1
CML#ifdef KPP_AUTODIFF_MORE_STORE
CMLCADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
CMLCADJ & key=kkppkey, kind=isbyte
CMLCADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
CMLCADJ & key=kkppkey, kind=isbyte
CML#endif /* KPP_AUTODIFF_MORE_STORE */
DO K = 1, Nr
Km1 = MAX(K-1,1)
DO J = 1-OLy, sNy+OLy
DO I = 1-OLx, sNx+OLx
alphaDT(I,J) = ( theta(I,J,Km1,bi,bj)-theta(I,J,K,bi,bj) )
& * 0.5 _d 0 * ABS( TTALPHA(I,J,Km1) + TTALPHA(I,J,K) )
betaDS(I,J) = ( salt(I,J,Km1,bi,bj)-salt(I,J,K,bi,bj) )
& * 0.5 _d 0 * ( SSBETA(I,J,Km1) + SSBETA(I,J,K) )
nuddt(I,J) = 0. _d 0
nudds(I,J) = 0. _d 0
ENDDO
ENDDO
IF ( K .GT. 1 ) THEN
DO J = jMin, jMax
DO I = iMin, iMax
Rrho = 0. _d 0
C Now we have many different cases
C a. alphaDT > 0 and betaDS > 0 => salt fingering
C (salinity destabilizes)
IF ( alphaDT(I,J) .GT. betaDS(I,J)
& .AND. betaDS(I,J) .GT. 0. _d 0 ) THEN
Rrho = MIN( alphaDT(I,J)/betaDS(I,J), Rrho0 )
C Large et al. 1994, eq. 31a
C nudds(I,J) = dsfmax * ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )**3
nutmp = ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )
nudds(I,J) = dsfmax * nutmp * nutmp * nutmp
C Large et al. 1994, eq. 31c
nuddt(I,J) = 0.7 _d 0 * nudds(I,J)
ELSEIF ( alphaDT(I,J) .LT. 0. _d 0
& .AND. betaDS(I,J) .LT. 0. _d 0
& .AND.alphaDT(I,J) .GT. betaDS(I,J) ) THEN
C b. alphaDT < 0 and betaDS < 0 => semi-convection, diffusive convection
C (temperature destabilizes)
C for Rrho >= 1 the water column is statically unstable and we never
C reach this point
Rrho = alphaDT(I,J)/betaDS(I,J)
C Large et al. 1994, eq. 32
nuddt(I,J) = numol * 0.909 _d 0
& * exp ( 4.6 _d 0 * exp (
& - 5.4 _d 0 * ( 1. _d 0/Rrho - 1. _d 0 ) ) )
CMLC or
CMLC Large et al. 1994, eq. 33
CML nuddt(I,J) = numol * 8.7 _d 0 * Rrho**1.1
C Large et al. 1994, eqs. 34
nudds(I,J) = nuddt(I,J) * MAX( 0.15 _d 0 * Rrho,
& 1.85 _d 0 * Rrho - 0.85 _d 0 )
ELSE
C Do nothing, because in this case the water colume is unstable
C => double diffusive processes are negligible and mixing due
C to shear instability will dominate
ENDIF
ENDDO
ENDDO
C ENDIF ( K .GT. 1 )
ENDIF
C
DO J = 1-OLy, sNy+OLy
DO I = 1-OLx, sNx+OLx
kappaRT(I,J,K) = kappaRT(I,J,K) + nuddt(I,J)
kappaRS(I,J,K) = kappaRS(I,J,K) + nudds(I,J)
ENDDO
ENDDO
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(nuddt,'KPPnuddt',k,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(nudds,'KPPnudds',k,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C end of K-loop
ENDDO
#endif /* ALLOW_KPP */
RETURN
END