unstable in p-ph dyn

obl_pph.f90

@@ -42,7 +42,7 @@ module obl_pph
         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 :: Kh_unstable = 100.0 !< unstable eddy diffusivity, [m**2 / s], *:
 
         real :: kappa = 0.4         !< von Karman constant, [n/d]
         real :: PrT0 = 1.0          !< neutral Prandtl number, [n/d]

obl_pph_dyn.f90

@@ -39,6 +39,9 @@ module obl_pph_dyn
         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], *:
+
         real :: C_l = 4.0 / 9.0     !< length scale tuning constant, [n/d], *: = (4.0 / pi^2) in Obukhov length scale
         
         real :: kappa = 0.4         !< von Karman constant
@@ -150,7 +153,11 @@ module obl_pph_dyn
         Km_0 = param%C_l * U_dynH * param%kappa * mld / 2.0
 
         do k = 1, cz
+            if (Ri_grad(k) >= 0.0) then
                 Km(k) = 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
 
@@ -190,7 +197,11 @@ module obl_pph_dyn
         Kh_0 = param%C_l * U_dynH * param%kappa * mld / 2.0
 
         do k = 1, cz
+            if (Ri_grad(k) >= 0.0) then
                 Kh(k) = 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
 
@@ -241,6 +252,22 @@ call c_config_is_varname(trim(tag)//".alpha"//C_NULL_CHAR, status)
                 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)//".C_l"//C_NULL_CHAR, status)
         if (status /= 0) then
             call c_config_get_float(trim(tag)//".C_l"//C_NULL_CHAR, param%C_l, status)