From d734536c0f4b3e5213bb9088694725421b940107 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 17 Aug 2021 08:04:41 -0600 Subject: [PATCH] Refactor the way Langmuir entrainment and enhancement factor are computed 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. --- .../vertical/MOM_CVMix_KPP.F90 | 251 +++++++----------- 1 file changed, 98 insertions(+), 153 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index cdd205fabf..0abc5d1b2a 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -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) @@ -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 @@ -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 , & @@ -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, & @@ -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)) @@ -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 @@ -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 @@ -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]