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]