Km, Kh unstable added for Pacanowski-Philander

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
+            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
+            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)