diff --git a/obl_pph.f90 b/obl_pph.f90 index 74c3d75ef97bca832bdcb5a10d1c98b4acbe5702..ceebd3de6c316c919631464b95e413f331461b02 100644 --- a/obl_pph.f90 +++ b/obl_pph.f90 @@ -41,10 +41,13 @@ module obl_pph real :: Km_b = 0.0001 !< background eddy viscosity, [m**2 / s], *: INMCM: 5 * 10**(-6) real :: Kh_b = 0.00001 !< background eddy diffusivity, [m**2 / s] + real :: Km_unstable = 100.0 !< unstable eddy viscosity, [m**2 / s] *: + real :: Kh_unstable = 100.0 !< unstable eddy diffusivity, [m**2 / s], *: INMCM: 5.0 * 0.01 + real :: kappa = 0.4 !< von Karman constant, [n/d] real :: PrT0 = 1.0 !< neutral Prandtl number, [n/d] - real :: c_S2_min = 1e-4 !< min shear frequency, [(1/s)**2] + real :: c_S2_min = 1e-5 !< min shear frequency, [(1/s)**2] end type contains @@ -129,7 +132,11 @@ module obl_pph ! -------------------------------------------------------------------------------- do k = 1, n - Km(k) = param%Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b + if (Ri_grad(k) >= 0.0) then + Km(k) = param%Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b + else + Km(k) = param%Km_unstable + endif end do end subroutine @@ -160,7 +167,11 @@ module obl_pph ! -------------------------------------------------------------------------------- do k = 1, n - Kh(k) = param%Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b + if (Ri_grad(k) >= 0.0) then + Kh(k) = param%Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b + else + Kh(k) = param%Kh_unstable + endif end do end subroutine @@ -227,6 +238,22 @@ module obl_pph return end if endif + call c_config_is_varname(trim(tag)//".Km_unstable"//C_NULL_CHAR, status) + if (status /= 0) then + call c_config_get_float(trim(tag)//".Km_unstable"//C_NULL_CHAR, param%Km_unstable, status) + if (status == 0) then + ierr = 1 ! signal ERROR + return + end if + endif + call c_config_is_varname(trim(tag)//".Kh_unstable"//C_NULL_CHAR, status) + if (status /= 0) then + call c_config_get_float(trim(tag)//".Kh_unstable"//C_NULL_CHAR, param%Kh_unstable, status) + if (status == 0) then + ierr = 1 ! signal ERROR + return + end if + endif call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status) if (status /= 0) then call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)