diff --git a/Turbulence_PrognosticTKE.pyx b/Turbulence_PrognosticTKE.pyx index e0201e19..8493e3df 100755 --- a/Turbulence_PrognosticTKE.pyx +++ b/Turbulence_PrognosticTKE.pyx @@ -1026,13 +1026,23 @@ cdef class EDMF_PrognosticTKE(ParameterizationBase): cpdef compute_horizontal_eddy_diffusivities(self, GridMeanVariables GMV): cdef: Py_ssize_t i, k - double l, R_up + double R_up,wu_half, we_half, a, velocity_scale + double [:] ae = np.subtract(np.ones((self.Gr.nzg,),dtype=np.double, order='c'),self.UpdVar.Area.bulkvalues) + double l[2] with nogil: for k in xrange(self.Gr.gw, self.Gr.nzg-self.Gr.gw): for i in xrange(self.n_updrafts): if self.UpdVar.Area.values[i,k]>0.0: - self.horizontal_KM[i,k] = self.UpdVar.Area.values[i,k]*self.turbulent_entrainment_factor*sqrt(fmax(GMV.TKE.values[k],0.0))*self.pressure_plume_spacing[i] + wu_half = interp2pt(self.UpdVar.W.values[i,k], self.UpdVar.W.values[i,k-1]) + we_half = interp2pt(self.EnvVar.W.values[k], self.EnvVar.W.values[k-1]) + a = self.UpdVar.Area.values[i,k] + with gil: + l[0] = sqrt(fmax(ae[k]*self.EnvVar.TKE.values[k],0.0)) + l[1] = sqrt(fmax(a*ae[k]*(wu_half-we_half)**2.0,0.0)) + velocity_scale = lamb_smooth_minimum(l, 0.1, 0.05) + self.horizontal_KM[i,k] = self.UpdVar.Area.values[i,k]*self.turbulent_entrainment_factor \ + *velocity_scale*self.pressure_plume_spacing[i] self.horizontal_KH[i,k] = self.horizontal_KM[i,k] / self.prandtl_nvec[k] else: self.horizontal_KM[i,k] = 0.0 @@ -1268,17 +1278,11 @@ cdef class EDMF_PrognosticTKE(ParameterizationBase): wu_half = interp2pt(self.UpdVar.W.values[i,k], self.UpdVar.W.values[i,k-1]) we_half = interp2pt(self.EnvVar.W.values[k], self.EnvVar.W.values[k-1]) if a*wu_half > 0.0: - ed_mf_ratio = fabs(self.EnvVar.TKE.buoy[k])/(fabs(a*(1.0-a)* - (wu_half-we_half)*(self.UpdVar.B.values[i,k] - self.EnvVar.B.values[k]))+1e-8) - self.turb_entr_H[i,k] = (2.0/R_up**2.0)*self.Ref.rho0_half[k] * a * self.horizontal_KH[i,k] * \ - (self.EnvVar.H.values[k] - self.UpdVar.H.values[i,k]) * \ - (1.0/(1.0+exp(self.entrainment_ed_mf_sigma*(ed_mf_ratio-1.0)))) + (self.EnvVar.H.values[k] - self.UpdVar.H.values[i,k]) self.turb_entr_QT[i,k] = (2.0/R_up**2.0)*self.Ref.rho0_half[k]* a * self.horizontal_KH[i,k] * \ - (self.EnvVar.QT.values[k] - self.UpdVar.QT.values[i,k]) * \ - (1.0/(1.0+exp(self.entrainment_ed_mf_sigma*(ed_mf_ratio-1.0)))) - self.frac_turb_entr[i,k] = (2.0/R_up**2.0) * self.horizontal_KH[i,k] * \ - (1.0/(1.0+exp(self.entrainment_ed_mf_sigma*(ed_mf_ratio-1.0))))/ wu_half/a + (self.EnvVar.QT.values[k] - self.UpdVar.QT.values[i,k]) + self.frac_turb_entr[i,k] = (2.0/R_up**2.0) * self.horizontal_KH[i,k] / wu_half/a else: self.turb_entr_H[i,k] = 0.0 self.turb_entr_QT[i,k] = 0.0 @@ -1289,14 +1293,9 @@ cdef class EDMF_PrognosticTKE(ParameterizationBase): b_env_full = interp2pt(self.EnvVar.B.values[k], self.EnvVar.B.values[k-1]) env_tke_full = interp2pt(self.EnvVar.TKE.buoy[k], self.EnvVar.TKE.buoy[k-1]) - ed_mf_ratio = fabs(env_tke_full)/(fabs(a_full*(1.0-a_full)* - (self.UpdVar.W.values[i,k]-self.EnvVar.W.values[k])*(b_upd_full - b_env_full))+1e-8) - self.turb_entr_W[i,k] = (2.0/R_up_full**2.0)*self.Ref.rho0[k] * a_full * K_full * \ - (self.EnvVar.W.values[k]-self.UpdVar.W.values[i,k]) * \ - (1.0/(1.0+exp(self.entrainment_ed_mf_sigma*(ed_mf_ratio-1.0)))) - self.frac_turb_entr_full[i,k] = (2.0/R_up_full**2.0) * K_full / self.UpdVar.W.values[i,k] * \ - (1.0/(1.0+exp(self.entrainment_ed_mf_sigma*(ed_mf_ratio-1.0))))/a_full + (self.EnvVar.W.values[k]-self.UpdVar.W.values[i,k]) + self.frac_turb_entr_full[i,k] = (2.0/R_up_full**2.0) * K_full / self.UpdVar.W.values[i,k] / a_full else: self.turb_entr_W[i,k] = 0.0 diff --git a/generate_paramlist.py b/generate_paramlist.py index 9c3b2088..82e31b25 100644 --- a/generate_paramlist.py +++ b/generate_paramlist.py @@ -45,9 +45,9 @@ def main(): paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['static_stab_coeff'] = 0.4 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['lambda_stab'] = 0.9 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['max_area'] = 0.9 - paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['entrainment_factor'] = 0.11 - paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['detrainment_factor'] = 0.5 - paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['turbulent_entrainment_factor'] = 0.06 + paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['entrainment_factor'] = 0.13 + paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['detrainment_factor'] = 0.51 + paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['turbulent_entrainment_factor'] = 0.08 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['entrainment_ed_mf_sigma'] = 50.0 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['entrainment_smin_tke_coeff'] = 0.3 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['updraft_mixing_frac'] = 0.25 @@ -62,11 +62,10 @@ def main(): # TODO: merge the tan18 buoyancy forluma into normalmode formula -> simply set buoy_coeff1 as 1./3. and buoy_coeff2 as 0. paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_buoy_coeff'] = 1.0/3.0 - - paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_buoy_coeff1'] = 0.22 + paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_buoy_coeff1'] = 0.12 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_buoy_coeff2'] = 0.0 paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_adv_coeff'] = 0.1 - paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_drag_coeff'] = 9.0 + paramlist_defaults['turbulence']['EDMF_PrognosticTKE']['pressure_normalmode_drag_coeff'] = 10.0 if case_name == 'Soares': paramlist = Soares(paramlist_defaults) diff --git a/turbulence_functions.pyx b/turbulence_functions.pyx index af860db4..a3bd93c3 100644 --- a/turbulence_functions.pyx +++ b/turbulence_functions.pyx @@ -104,7 +104,7 @@ cdef entr_struct entr_detr_env_moisture_deficit(entr_in_struct entr_in) nogil: with gil: l[0] = entr_in.tke_coef*fabs(db/sqrt(entr_in.tke+1e-8)) l[1] = fabs(db/dw) - inv_timescale = lamb_smooth_minimum(l, 0.0001, 1e-3) + inv_timescale = lamb_smooth_minimum(l, 0.1, 0.0005) _ret.entr_sc = inv_timescale/dw*(entr_in.c_ent*logistic_e + c_det*moisture_deficit_e) _ret.detr_sc = inv_timescale/dw*(entr_in.c_ent*logistic_d + c_det*moisture_deficit_d)