pacanowski-philander module update

obl_main.f90

@@ -42,6 +42,9 @@ program obl_main
 
     type(gridDataType) :: grid
 
+    ! turbulence closure parameters
+    type(pacanowskiParamType) :: param_pacanowski
+
 #ifdef OBL_USE_TSLICE_OUTPUT
     ! time slice
     type (oblTimeSlice) :: output_Theta, output_Salin
@@ -362,8 +365,9 @@ program obl_main
             call get_dyn_velocity(u_dynH, flux_u_surf_res, flux_v_surf_res)
             call get_lab_mld(lab_mld, u_dynH, N2_0, time_current, z)
             if (closure_mode.eq.1) then
-                call get_Km(Km, Ri_grad, nz)
-                call get_Kh(Kh, Ri_grad, nz)
+                call get_eddy_viscosity(Km, Ri_grad, param_pacanowski, nz)
+                call get_eddy_diffusivity(Kh, Ri_grad, param_pacanowski, nz)
+                
                 call get_TKE_generation(P_TKE, Km, S2, nz)
                 call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
             else if (closure_mode.eq.2) then

obl_pacanowski_philander.f90

@@ -7,13 +7,6 @@ module obl_pacanowski_philander
     ! directives list
     implicit none
 
-    real, parameter :: Km_0 = 2.0 * 0.01 !< eddy viscosity for momentum, [m**2 / s]   7.0 * 0.01 
-    real, parameter :: alpha_turb = 5.0 !< constant
-    real, parameter :: n_turb = 2.0 !< constant
-    real, parameter :: Km_b =  0.0001 !< background eddy viscosity for momentum, [m**2 / s] - 5 * 10**(-6) in inmcm 
-    real, parameter :: Kh_0 = 2.0 * 0.01 !< background eddy diffusivity for tracers, [m**2 / s] 5.0 * 0.01 
-    real, parameter :: Kh_b = 0.00001 !< background eddy diffusivity for tracers, [m**2 / s]
-
     !> @brief Pacanowski-Philander parameters
     type, public :: pacanowskiParamType
         real :: Km_0 = 2.0 * 0.01   !< neutral eddy viscosity, [m**2 / s] *: INMCM: 7.0 * 0.01
@@ -26,18 +19,6 @@ module obl_pacanowski_philander
 
     contains
     
-    subroutine get_Km(Km, Ri_grad, nz)
-    !< @brief calculate momentum transfer coefficient in P-Ph closure
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(in) :: Ri_grad(nz) !< Richardson gradient number 
-        real, intent(out) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
-        integer :: k !< counter
-
-        do k = 1, nz
-            Km(k) = Km_0 / (1.0 + alpha_turb * Ri_grad(k))**(n_turb) + Km_b
-        end do
-    end subroutine
-
     subroutine get_eddy_viscosity(Km, Ri_grad, param, nz)
         !< @brief calculate eddy viscosity
         type(pacanowskiParamType), intent(in) :: param
@@ -49,20 +30,23 @@ module obl_pacanowski_philander
         integer :: k
     
         do k = 1, nz
-            Km(k) = Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
+            Km(k) = param%Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
         end do
     end subroutine
 
-    subroutine get_Kh(Kh, Ri_grad, nz)
+
+    subroutine get_eddy_diffusivity(Kh, Ri_grad, param, nz)
         !< @brief calculate scalar transfer coefficient in in P-Ph+ closure
+        type(pacanowskiParamType), intent(in) :: param
+
         integer, intent(in) :: nz                       !< number of z-steps
-        real, intent(in) :: Ri_grad(nz)  !< Richardson gradient number 
-        real, intent(out) :: Kh(nz) !eddy diffusivity for tracers, [m**2 / s]
-        integer :: k !< counter
+        real, dimension(nz), intent(in) :: Ri_grad      !< Richardson gradient number 
+        real, dimension(nz), intent(out) :: Kh          !eddy diffusivity for tracers, [m**2 / s]
+        
+        integer :: k
 
         do k = 1, nz
-            Kh(k) = Kh_0 / (1.0 + alpha_turb * Ri_grad(k))**(n_turb + 1.0) &
-                + Kh_b
+            Kh(k) = param%Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b
         end do
     end subroutine