Skip to content

Commit

Permalink
Refactor the way Langmuir entrainment and enhancement factor are comp…
Browse files Browse the repository at this point in the history 
…uted and applied: (1) When available, pass the enhancement factor received from WW3. Otherwise, let CVmix compute the enhancement factor. (2) Instead of explicitly multiplying diffusivities and viscosities with the enhancement factor, let CVMix handle enhancement internally.
  • Loading branch information
@alperaltuntas
alperaltuntas committed Aug 17, 2021
1 parent 734e537 commit d734536
Showing 1 changed file with 98 additions and 153 deletions.
251 changes: 98 additions & 153 deletions src/parameterizations/vertical/MOM_CVMix_KPP.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,8 +198,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
# include "version_variable.h"
character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
character(len=20) :: string !< local temporary string
character(len=20) :: langmuir_mixing_opt !< langmuir mixing opt to be passed to CVMix, e.g., LWF16
character(len=20) :: langmuir_entrainment_opt !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16
character(len=20) :: langmuir_mixing_opt = 'NONE' !< langmuir mixing opt to be passed to CVMix, e.g., LWF16
character(len=20) :: langmuir_entrainment_opt = 'NONE' !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16
character(len=20) :: wave_method
logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local
logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function)
Expand Down Expand Up @@ -418,10 +418,16 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
'\t RW16 = Function of Langmuir number based on RW16', &
default='CONSTANT')
select case ( trim(string))
case ("CONSTANT") ; CS%LT_K_METHOD = LT_K_MODE_CONSTANT
case ("VR12") ; CS%LT_K_METHOD = LT_K_MODE_VR12
case ("RW16") ; CS%LT_K_METHOD = LT_K_MODE_RW16
case default ; call MOM_error(FATAL,"KPP_init: "//&
case ("CONSTANT")
CS%LT_K_METHOD = LT_K_MODE_CONSTANT
case ("VR12")
CS%LT_K_METHOD = LT_K_MODE_VR12
langmuir_mixing_opt = 'LWF16'
case ("RW16")
CS%LT_K_METHOD = LT_K_MODE_RW16
langmuir_mixing_opt = 'RWHGK16'
case default
call MOM_error(FATAL,"KPP_init: "//&
"Unrecognized KPP_LT_K_METHOD option: "//trim(string))
end select
if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then
Expand All @@ -444,12 +450,20 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
'\t LF17 = Function of Langmuir number based on LF17', &
default='CONSTANT')
select case ( trim(string))
case ("CONSTANT") ; CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT
case ("VR12") ; CS%LT_VT2_METHOD = LT_VT2_MODE_VR12
case ("RW16") ; CS%LT_VT2_METHOD = LT_VT2_MODE_RW16
case ("LF17") ; CS%LT_VT2_METHOD = LT_VT2_MODE_LF17
case default ; call MOM_error(FATAL,"KPP_init: "//&
"Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
case ("CONSTANT")
CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT
case ("VR12")
CS%LT_VT2_METHOD = LT_VT2_MODE_VR12
langmuir_entrainment_opt = 'LWF16'
case ("RW16")
CS%LT_VT2_METHOD = LT_VT2_MODE_RW16
langmuir_entrainment_opt = 'RWHGK16'
case ("LF17")
CS%LT_VT2_METHOD = LT_VT2_MODE_LF17
langmuir_entrainment_opt = 'LF17'
case default
call MOM_error(FATAL,"KPP_init: "//&
"Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
end select
if (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then
call get_param(paramFile, mdl, "KPP_VT2_ENH_FAC",CS%KPP_VT2_ENH_FAC , &
Expand All @@ -460,17 +474,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)

call closeParameterBlock(paramFile)

! since wave interface control structure is not initialized yet, temporarily read in wave_method from param file
call get_param(paramFile, mdl, "WAVE_METHOD", wave_method, default="EMPTY", do_not_log=.true.)
select case ( trim(wave_method) )
case ("VR12-MA")
langmuir_mixing_opt = 'LWF16'
langmuir_entrainment_opt = 'LWF16'
case default
langmuir_mixing_opt = 'NONE'
langmuir_entrainment_opt = 'NONE'
end select

call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.)

call CVMix_init_kpp( Ri_crit=CS%Ri_crit, &
Expand Down Expand Up @@ -732,62 +735,16 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
Kviscosity(:) = US%Z2_T_to_m2_s * Kv(i,j,:)
endif

if (get_wave_method(Waves) == "VR12-MA") then
call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
iFaceHeight, & ! (in) Height of interfaces [m]
cellHeight, & ! (in) Height of level centers [m]
Kviscosity(:), & ! (in) Original viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
CS%OBLdepth(i,j), & ! (in) OBL depth [m]
CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim]
nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim]
surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
GV%ke, & ! (in) Number of levels to compute coeffs for
GV%ke, & ! (in) Number of levels in array shape
Langmuir_EFactor=lamult(i,j),& ! Langmuir enhancement multiplier
CVMix_kpp_params_user=CS%KPP_params )
else
call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
iFaceHeight, & ! (in) Height of interfaces [m]
cellHeight, & ! (in) Height of level centers [m]
Kviscosity(:), & ! (in) Original viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
CS%OBLdepth(i,j), & ! (in) OBL depth [m]
CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim]
nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim]
surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
GV%ke, & ! (in) Number of levels to compute coeffs for
GV%ke, & ! (in) Number of levels in array shape
CVMix_kpp_params_user=CS%KPP_params )
endif

! safety check, Kviscosity and Kdiffusivity must be >= 0
do k=1, GV%ke+1
if (Kviscosity(k) < 0. .or. Kdiffusivity(k,1) < 0.) then
call MOM_error(FATAL,"KPP_calculate, after CVMix_coeffs_kpp: "// &
"Negative vertical viscosity or diffusivity has been detected. " // &
"This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //&
"You might consider using the default options for these parameters." )
endif
enddo

IF (CS%LT_K_ENHANCEMENT) then
if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then
LangEnhK = CS%KPP_K_ENH_FAC
elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then
! Added minimum value for La_SL, so removed maximum value for LangEnhK.
LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
if (present(lamult)) then
LangEnhK = lamult(i,j)
else
LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
(5.4*CS%La_SL(i,j))**(-4))
endif
elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then
!This maximum value is proposed in Reichl et al., 2016 JPO formula
LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j))
Expand All @@ -796,26 +753,59 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
!call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
LangEnhK = 1.0
endif

! diffusivities don't need to be enhanced below anymore since LangEnhK is applied within CVMix.
! todo: need to double check if the distinction between the two different options of LT_K_SHAPE may need to be
! treated specially.
do k=1,GV%ke
if (CS%LT_K_SHAPE== LT_K_CONSTANT) then
if (CS%id_EnhK > 0) CS%EnhK(i,j,:) = LangEnhK
Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK
Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK
Kviscosity(k) = Kviscosity(k) * LangEnhK
!Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK
!Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK
!Kviscosity(k) = Kviscosity(k) * LangEnhK
elseif (CS%LT_K_SHAPE == LT_K_SCALED) then
sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
SigmaRatio = sigma * (1. - sigma)**2 / 0.148148037
if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio)
Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + &
( LangEnhK - 1.)*sigmaRatio)
Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + &
( LangEnhK - 1.)*sigmaRatio)
Kviscosity(k) = Kviscosity(k) * ( 1. + &
( LangEnhK - 1.)*sigmaRatio)
!Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + &
! ( LangEnhK - 1.)*sigmaRatio)
!Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + &
! ( LangEnhK - 1.)*sigmaRatio)
!Kviscosity(k) = Kviscosity(k) * ( 1. + &
! ( LangEnhK - 1.)*sigmaRatio)
endif
enddo
endif

call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
iFaceHeight, & ! (in) Height of interfaces [m]
cellHeight, & ! (in) Height of level centers [m]
Kviscosity(:), & ! (in) Original viscosity [m2 s-1]
Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
CS%OBLdepth(i,j), & ! (in) OBL depth [m]
CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim]
nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim]
surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
GV%ke, & ! (in) Number of levels to compute coeffs for
GV%ke, & ! (in) Number of levels in array shape
Langmuir_EFactor=LangEnhK,& ! Langmuir enhancement multiplier
CVMix_kpp_params_user=CS%KPP_params )

! safety check, Kviscosity and Kdiffusivity must be >= 0
do k=1, GV%ke+1
if (Kviscosity(k) < 0. .or. Kdiffusivity(k,1) < 0.) then
call MOM_error(FATAL,"KPP_calculate, after CVMix_coeffs_kpp: "// &
"Negative vertical viscosity or diffusivity has been detected. " // &
"This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //&
"You might consider using the default options for these parameters." )
endif
enddo

! Over-write CVMix NLT shape function with one of the following choices.
! The CVMix code has yet to update for thse options, so we compute in MOM6.
! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
Expand Down Expand Up @@ -993,8 +983,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
integer :: kk, ksfc, ktmp

! For Langmuir Calculations
real :: LangEnhW ! Langmuir enhancement for turbulent velocity scale
real, dimension(GV%ke) :: LangEnhVt2 ! Langmuir enhancement for unresolved shear
real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear
real, dimension(GV%ke) :: U_H, V_H
real :: MLD_GUESS, LA
real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir
Expand Down Expand Up @@ -1186,87 +1175,43 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
CVMix_kpp_params_user=CS%KPP_params )

!Compute CVMix VT2
if (get_wave_method(Waves) == "VR12-MA") then
CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( &
zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m]
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
EFactor=lamult(i,j), & ! Langmuir enhancement multiplier
CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
else
CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( &
zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m]
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
endif

!Modify CVMix VT2
! Determine the enhancement factor for unresolved shear
IF (CS%LT_VT2_ENHANCEMENT) then
IF (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then
do k=1,GV%ke
LangEnhVT2(k) = CS%KPP_VT2_ENH_FAC
enddo
LangEnhVT2 = CS%KPP_VT2_ENH_FAC
elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then
!Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
enhvt2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
(5.4*CS%La_SL(i,j))**(-4))
do k=1,GV%ke
LangEnhVT2(k) = enhvt2
enddo
elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_RW16) then
!Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
enhvt2 = 1. + 2.3*CS%La_SL(i,j)**(-0.5)
do k=1,GV%ke
LangEnhVT2(k) = enhvt2
enddo
elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_LF17) then
CS%CS=cvmix_get_kpp_real('c_s',CS%KPP_params)
do k=1,GV%ke
WST = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellHeight(k)))**(1./3.)
LangEnhVT2(k) = sqrt((0.15*WST**3. + 0.17*surfFricVel**3.* &
(1.+0.49*CS%La_SL(i,j)**(-2.))) / &
(0.2*ws_1d(k)**3/(CS%cs*CS%surf_layer_ext*CS%vonKarman**4.)))
enddo
if (present(lamult)) then
LangEnhVT2 = lamult(i,j)
else
LangEnhVT2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
(5.4*CS%La_SL(i,j))**(-4))
endif
else
!This shouldn't be reached.
!call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2")
LangEnhVT2(:) = 1.0
! for other methods (e.g., LT_VT2_MODE_RW16, LT_VT2_MODE_LF17), the enhancement factor is
! computed internally within CVMix using LaSL, bfsfc, and ustar to be passed to CVMix.
LangEnhVT2 = 1.0
endif
else
LangEnhVT2(:) = 1.0
LangEnhVT2 = 1.0
endif

do k=1,GV%ke
CS%Vt2(i,j,k)=CS%Vt2(i,j,k)*LangEnhVT2(k)
if (CS%id_EnhVt2 > 0) CS%EnhVt2(i,j,k)=LangEnhVT2(k)
enddo
surfBuoyFlux = buoy_scale * buoyFlux(i,j,1)

! Calculate Bulk Richardson number from eq (21) of LMD94
if (get_wave_method(Waves) == "VR12-MA") then
BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( &
zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m]
delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2]
Vt_sqr_cntr=CS%Vt2(i,j,:), &
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1]
EFactor=lamult(i,j)) ! Langmuir enhancement factor
else
BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( &
zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m]
delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2]
Vt_sqr_cntr=CS%Vt2(i,j,:), &
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
N_iface=CS%N(i,j,:)) ! Buoyancy frequency [s-1]
endif
BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( &
zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m]
delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2]
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1]
EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim]
LaSL = CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim]
bfsfc = surfBuoyFlux, & ! surface buoyancy flux [m2 s-3]
uStar = uStar(i,j), & ! surface friction velocity [m s-1]
CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters


surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
! h to Monin-Obukov (default is false, ie. not used)

call CVMix_kpp_compute_OBL_depth( &
BulkRi_1d, & ! (in) Bulk Richardson number
iFaceHeight, & ! (in) Height of interfaces [m]
Expand Down

0 comments on commit d734536

Please sign in to comment.