Skip to content
Snippets Groups Projects
Commit 419fe261 authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

diagnostic module implementation update

parent a3aaac62
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_common_phys.f90
+ 1
− 0
View file @ 419fe261
module obl_common_phys module obl_common_phys
!< @brief common physical parameters and constants !< @brief common physical parameters and constants
! --------------------------------------------------------------------------------
implicit none implicit none
private private
......
obl_diag.f90
+ 60
− 74
View file @ 419fe261
...@@ -15,128 +15,114 @@ module obl_diag ...@@ -15,128 +15,114 @@ module obl_diag
! public interface ! public interface
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
public :: get_mld, get_lab_mld public :: get_mld, get_mld_ref
public :: get_eld, get_eld_ref
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
contains contains
subroutine get_mld(mld, N2, dz, nz, z) ! --------------------------------------------------------------------------------
subroutine get_mld(mld, N2, dz, cz)
!< @brief calculate mixed layer depth !< @brief calculate mixed layer depth
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
real, intent(out) :: mld !< mixed layer depth, [m] real, intent(out) :: mld !< mixed layer depth, [m]
integer, intent(in) :: cz !< number of z-steps
real, intent(in), dimension(cz) :: N2 !< square of buoyancy frequency, [s**(-2)]
real, intent(in), dimension(cz) :: dz !< grid steps [m]
integer, intent(in) :: nz !< number of z-steps integer :: k_max ! max(N2) grid cell index
real, intent(in), dimension(nz) :: N2 !< square of buoyancy frequency, [s**(-2)] real :: N2_max
real, intent(in), dimension(nz) :: dz !< depth(z) step [m]
real, intent(in) :: z !< depth, [m]
integer :: max_k !< index of cell with maximum N2 integer :: k
real :: max_N2, mld_tmp !< temporary
real :: mld_min, mld_max !< mld limits
integer :: k !< counter
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
mld_min = 0.5 * dz(nz)
mld_max = z
max_N2 = N2(nz) N2_max = N2(cz)
max_k = 1 k_max = cz
mld = 0.0
mld_tmp = 0.0
do k=nz - 1, 1, -1 do k = cz - 1, 1, -1
if (N2(k) > max_N2) then if (N2(k) > N2_max) then
max_N2 = N2(k) N2_max = N2(k)
max_k = k k_max = k
end if end if
end do end do
do k=nz, max_k + 1, -1 mld = 0.0
mld_tmp = mld_tmp + dz(k) do k = cz, k_max + 1, -1
mld = mld + dz(k)
enddo enddo
mld_tmp = mld_tmp + 0.5 * dz(max_k) mld = mld + 0.5 * dz(k_max)
mld = mld_tmp
if (mld < mld_min) mld = mld_min
if (mld > mld_max) mld = mld_max
end subroutine get_mld end subroutine get_mld
subroutine get_lab_mld(lab_mld, u_dynH, N2_0, time, z) ! --------------------------------------------------------------------------------
subroutine get_mld_ref(mld_ref, u_dynH, N2_0, time, max_depth)
!< @brief calculate theoretical mixed layer depth (Price, 1970) !< @brief calculate theoretical mixed layer depth (Price, 1970)
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
real, intent(out) :: mld_ref !< theoretical mixed layer depth, [m]
real, intent(out) :: lab_mld !< theoretical mixed layer depth, [m] real, intent(in) :: u_dynH !< surface dynamic velocity, [(m/s)]
real, intent(in) :: u_dynH !< u_dyn, [(m/s)]
real, intent(in) :: N2_0 !< square of buoyancy frequency at t=0, [s**(-2)] real, intent(in) :: N2_0 !< square of buoyancy frequency at t=0, [s**(-2)]
real, intent(in) :: time !< time, [s] real, intent(in) :: time !< time, [s]
real, intent(in) :: z !<depth, [m] real, intent(in) :: max_depth !< max depth, [m]
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
lab_mld = 1.05 * u_dynH * sqrt(time) / (N2_0**(1.0/4.0)) ! to fix N2
if (lab_mld > z) lab_mld = z mld_ref = 1.05 * u_dynH * sqrt(time) / (N2_0**(1.0/4.0))
if (mld_ref > max_depth) mld_ref = max_depth
end subroutine get_lab_mld end subroutine get_mld_ref
subroutine get_eld(eld, TKE_B, dz, nz, z) ! --------------------------------------------------------------------------------
subroutine get_eld(eld, B, dz, cz)
!< @brief calculate entrainment layer depth !< @brief calculate entrainment layer depth
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
real, intent(out) :: eld !< mixed layer depth, [m] real, intent(out) :: eld !< entrainment layer depth, [m]
integer, intent(in) :: cz !< number of z-steps
real, intent(in), dimension(cz) :: B !< square of buoyancy frequency, [s**(-2)]
real, intent(in), dimension(cz) :: dz !< depth(z) step [m]
integer, intent(in) :: nz !< number of z-steps integer :: k_min ! min(B) grid cell index
real, intent(in), dimension(nz) :: TKE_B !< square of buoyancy frequency, [s**(-2)] real :: B_min
real, intent(in), dimension(nz) :: dz !< depth(z) step [m]
real, intent(in) :: z !< depth, [m]
integer :: min_k !< index of cell with minimum TKE_B integer :: k
real :: min_TKE_B, eld_tmp !< temporary
real :: eld_min, eld_max !< mld limits
integer :: k !< counter
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
eld_min = 0.5 * dz(nz)
eld_max = z
min_TKE_B = TKE_B(nz) B_min = B(cz)
min_k = 1 k_min = cz
eld = 0.0
eld_tmp = 0.0
do k=nz - 1, 1, -1 do k = cz - 1, 1, -1
if (TKE_B(k) < min_TKE_B) then if (B(k) < B_min) then
min_TKE_B = TKE_B(k) B_min = B(k)
min_k = k k_min = k
end if end if
end do end do
do k=nz, min_k + 1, -1 eld = 0.0
eld_tmp = eld_tmp + dz(k) do k = cz, k_min + 1, -1
eld = eld + dz(k)
enddo enddo
eld = eld + 0.5 * dz(k_min)
eld_tmp = eld_tmp + 0.5 * dz(min_k)
eld = eld_tmp
if (eld < eld_min) eld = eld_min
if (eld > eld_max) eld = eld_max
end subroutine get_eld end subroutine get_eld
subroutine get_lab_eld(lab_eld, N2_0, TKE_B_0, time, z) ! --------------------------------------------------------------------------------
subroutine get_eld_ref(eld_ref, N2_0, Bsurf, time, max_depth)
!< @brief calculate theoretical entrainment layer depth (Turner, 1972) !< @brief calculate theoretical entrainment layer depth (Turner, 1972)
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
real, intent(out) :: lab_eld !< theoretical entrainment layer depth, [m] real, intent(out) :: eld_ref !< theoretical entrainment layer depth, [m]
real, intent(in) :: N2_0 !< square of buoyancy frequency at t=0, [s**(-2)] real, intent(in) :: N2_0 !< square of buoyancy frequency at t=0, [s**(-2)]
real, intent(in) :: TKE_B_0 !<buoyancy production of TKE at t=0, [m**2 / s**3] real, intent(in) :: Bsurf !< buoyancy production of TKE at surface and at t=0, [m**2 / s**3]
real, intent(in) :: time !< time, [s] real, intent(in) :: time !< time, [s]
real, intent(in) :: z !<depth, [m] real, intent(in) :: max_depth !< max depth, [m]
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
lab_eld = sqrt ((2.0 * TKE_B_0 * time) / N2_0)
if (lab_eld > z) lab_eld = z eld_ref = sqrt((2.0 * Bsurf * time) / N2_0)
if (eld_ref > max_depth) eld_ref = max_depth
end subroutine get_lab_eld end subroutine get_eld_ref
endmodule endmodule
\ No newline at end of file
obl_k_epsilon.f90
+ 3
− 3
View file @ 419fe261
...@@ -37,9 +37,6 @@ module obl_k_epsilon ...@@ -37,9 +37,6 @@ module obl_k_epsilon
type (gridDataType), intent(in) :: grid type (gridDataType), intent(in) :: grid
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
call get_N2(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, grid) call get_N2(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, grid)
call get_S2(S2, U, V, & call get_S2(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid) bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
...@@ -58,6 +55,9 @@ module obl_k_epsilon ...@@ -58,6 +55,9 @@ module obl_k_epsilon
real, intent(in) :: dt real, intent(in) :: dt
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
call TKE_bc(TKE, & call TKE_bc(TKE, &
bc%U_dyn0, bc%U_dynH, param, grid%cz) bc%U_dyn0, bc%U_dynH, param, grid%cz)
call EPS_bc(EPS, & call EPS_bc(EPS, &
......
obl_main.f90
+ 4
− 4
View file @ 419fe261
...@@ -384,8 +384,8 @@ program obl_main ...@@ -384,8 +384,8 @@ program obl_main
!> advance screen output !> advance screen output
if (mod(i, nscreen) == 0) then if (mod(i, nscreen) == 0) then
call get_mld(mld, N2, grid%dz, grid%cz, grid%height) call get_mld(mld, N2, grid%dz, grid%cz)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, grid%height) call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
write(*, '(a,g0)') ' mld = ', mld write(*, '(a,g0)') ' mld = ', mld
write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref
...@@ -398,8 +398,8 @@ program obl_main ...@@ -398,8 +398,8 @@ program obl_main
if (mod(i, noutput) == 0) then if (mod(i, noutput) == 0) then
call push_state_vec(grid%cz) call push_state_vec(grid%cz)
call get_mld(mld, N2, grid%dz, grid%cz, grid%height) call get_mld(mld, N2, grid%dz, grid%cz)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, grid%height) call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
call push_value_to_tseries(output_mld, mld) call push_value_to_tseries(output_mld, mld)
call push_value_to_tseries(output_mld_ref, lab_mld) call push_value_to_tseries(output_mld_ref, lab_mld)
......
obl_pph_dyn.f90
+ 1
− 1
View file @ 419fe261
...@@ -56,7 +56,7 @@ module obl_pph_dyn ...@@ -56,7 +56,7 @@ module obl_pph_dyn
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz) call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz) call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
call get_mld(mld, N2, grid%dz, grid%cz, grid%height) call get_mld(mld, N2, grid%dz, grid%cz)
call get_Km_plus(Km, Ri_grad, & call get_Km_plus(Km, Ri_grad, &
bc%U_dynH, mld, param, grid%dz, grid%height, grid%cz) bc%U_dynH, mld, param, grid%dz, grid%height, grid%cz)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment