pacanowski-philander dynamic closure update

obl_diag.f90

@@ -2,7 +2,11 @@ module obl_diag
     !< @brief diagnostics module
     ! --------------------------------------------------------------------------------
 
+    ! TO DO:
+    !   -- add mld def. using momentum flux
+
     ! modules used
+    ! --------------------------------------------------------------------------------
     use obl_grid
     use obl_math
     use obl_common_phys
@@ -10,6 +14,7 @@ module obl_diag
     use obl_turb_common
 
     ! directives list
+    ! --------------------------------------------------------------------------------
     implicit none 
     private
 

obl_pph_dyn.f90

@@ -20,13 +20,13 @@ module obl_pph_dyn
     public :: advance_turbulence_closure
     public :: define_stability_functions
 
-    !public :: get_eddy_viscosity, get_eddy_diffusivity
+    public :: get_eddy_viscosity, get_eddy_diffusivity
     public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
     ! --------------------------------------------------------------------------------
 
     !> @brief Pacanowski-Philander dynamic parameters
     type, public :: pphDynParamType
-        real :: alpha = 5.0         !< constant
+        real :: alpha = 4.0         !< constant, *: default = 5.0
         real :: n = 2.0             !< constant
         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]
@@ -34,6 +34,8 @@ module obl_pph_dyn
         real :: kappa = 0.4         !< von Karman constant
         real :: PrT0 = 1.0          !< turbulent Prandtl number, *: probably redundant
 
+        real :: C_l = 4.0 / 9.0     !< length scale tuning constant, [n/d], *: = (4.0 / pi^2) in Obukhov length scale
+
         real :: c_S2_min = 1e-5     !< min shear frequency, [(1/s)**2]
     end type
 
@@ -77,57 +79,89 @@ module obl_pph_dyn
 
         call get_mld(mld, N2, grid%dz, grid%cz)
 
-        call get_eddy_viscosity_vec(Km, Ri_grad, &
-            bc%U_dynH, mld, param, grid%dz, grid%height, grid%cz)
-        call get_eddy_diffusivity_vec(Kh, Km, param, grid%cz)
+        call get_eddy_viscosity(Km, Ri_grad, &
+            bc%U_dynH, mld, param, grid)
+        call get_eddy_diffusivity(Kh, Ri_grad, &
+            bc%U_dynH, mld, param, grid)
     end subroutine
 
     ! --------------------------------------------------------------------------------
-    subroutine get_eddy_viscosity_vec(Km, Ri_grad, U_dynH, mld, param, dz, max_depth, nz)
+    subroutine get_eddy_viscosity(Km, Ri_grad, U_dynH, mld, param, grid)
         !< @brief calculate eddy viscosity
         ! ----------------------------------------------------------------------------
         type(pphDynParamType), intent(in) :: param
+        type (gridDataType), intent(in) :: grid
+
+        real, intent(out) :: Km(grid%cz)     !< eddy viscosity, [m**2 / s]
+        real, intent(in) :: Ri_grad(grid%cz) !< Richardson gradient number, [n/d] 
 
-        real, intent(out) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(in) :: Ri_grad(nz) !< Richardson gradient number 
-        real, intent(in) :: dz(nz) !< depth(z) step [m] 
-        real, intent(in) :: max_depth !< dept
         real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
         real, intent(in) :: mld         !< mixed layer depth, [m]
-        real :: mld_tmp, mld_max, mld_min !< temporary for mld, [m]
-        real :: Km_0 !< eddy viscosity for momentum, [m**2 / s]
-        integer :: k !< counter
         ! ----------------------------------------------------------------------------
         
-        mld_max = max_depth / 3.1416
-        mld_min = dz(nz)
+        call get_eddy_viscosity_vec(Km, Ri_grad, U_dynH, mld, param, grid%cz)
+    end subroutine
 
-        mld_tmp = mld 
+    subroutine get_eddy_viscosity_vec(Km, Ri_grad, U_dynH, mld, param, cz)
+        !< @brief calculate eddy viscosity
+        ! ----------------------------------------------------------------------------
+        type(pphDynParamType), intent(in) :: param
 
-        if (mld < mld_min) mld_tmp = mld_min
-        if (mld > mld_max) mld_tmp = mld_max
+        integer, intent(in) :: cz       !< grid size 
+        real, intent(out) :: Km(cz)     !< eddy viscosity, [m**2 / s]
+        real, intent(in) :: Ri_grad(cz) !< Richardson gradient number, [n/d] 
 
-        Km_0 = U_dynH * param%kappa * mld_tmp / 2.0
+        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
+        real, intent(in) :: mld         !< mixed layer depth, [m]
+
+        real :: Km_0                    !< neutral eddy viscosity, [m**2 / s]
+        integer :: k
+        ! ----------------------------------------------------------------------------
+        
+        Km_0 = param%C_l * U_dynH * param%kappa * mld / 2.0
 
-        do k = 1, nz
+        do k = 1, cz
             Km(k) = Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
         end do
     end subroutine
 
     ! --------------------------------------------------------------------------------
-    subroutine get_eddy_diffusivity_vec(Kh, Km, param, nz)
-        !< @brief calculate scalar transfer coefficient in in P-Ph+ closure
+    subroutine get_eddy_diffusivity(Kh, Ri_grad, U_dynH, mld, param, grid)
+        !< @brief calculate eddy diffusivity
         ! ----------------------------------------------------------------------------
         type(pphDynParamType), intent(in) :: param
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(in) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
-        real, intent(out) :: Kh(nz) !eddy diffusivity for tracers, [m**2 / s]
-        integer :: k !< counter
+        type (gridDataType), intent(in) :: grid
+
+        real, intent(out) :: Kh(grid%cz)     !< eddy diffusivity, [m**2 / s]
+        real, intent(in) :: Ri_grad(grid%cz) !< Richardson gradient number, [n/d]
+
+        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
+        real, intent(in) :: mld         !< mixed layer depth, [m]
+        ! ----------------------------------------------------------------------------
+
+        call get_eddy_diffusivity_vec(Kh, Ri_grad, U_dynH, mld, param, grid%cz)
+    end subroutine
+
+    subroutine get_eddy_diffusivity_vec(Kh, Ri_grad, U_dynH, mld, param, cz)
+        !< @brief calculate eddy diffusivity
+        ! ----------------------------------------------------------------------------
+        type(pphDynParamType), intent(in) :: param
+
+        integer, intent(in) :: cz       !< grid size 
+        real, intent(out) :: Kh(cz)     !< eddy diffusivity, [m**2 / s]
+        real, intent(in) :: Ri_grad(cz) !< Richardson gradient number, [n/d]
+
+        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
+        real, intent(in) :: mld         !< mixed layer depth, [m]
+
+        real :: Kh_0                    !< neutral eddy diffusivity, [m**2 / s]
+        integer :: k
         ! ----------------------------------------------------------------------------
 
-        do k = 1, nz
-            Kh(k) = Km(k) / param%PrT0
+        Kh_0 = param%C_l * U_dynH * param%kappa * mld / 2.0
+
+        do k = 1, cz
+            Kh(k) = Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b
         end do
     end subroutine
 

obl_state.f90

@@ -9,6 +9,7 @@ module obl_state
     ! --------------------------------------------------------------------------------
 
     ! directives list
+    ! --------------------------------------------------------------------------------
     implicit none
     private