module obl_legacy
implicit none
!!! OCEVERMIX LEGACY PARAMETERS
real, parameter :: lu_min = 0.5
! legacy_rit_top
real, parameter :: coriolis_rit_top_min = 2.5e-05
! legacy_rit
real, parameter :: rit_max = 1E+05
real, parameter :: s2_add = 0.0025
! legacy_kh_unstable
real, parameter :: kh_ttmx = 10.0
real, parameter :: kh_ttmn = 5.0
real, parameter :: kh_unstmx = 500.0
real, parameter :: kh_unstmn = 500.0 !50
! legacy_undimdepth & legacy_str
real, parameter :: hm_ice_max = 0.7
real, parameter :: hhq_add = 1.0e-30
real, parameter :: kh_strmin1 = 0.5
real, parameter :: kh_strmax1 = 2.0
real, parameter :: kh_strmin2 = 2.0
real, parameter :: kh_strmax2 = 10.0
real, parameter :: kh_hmmin1 = 300.0
real, parameter :: kh_hmmax1 = 3000.0
real, parameter :: kh_hmmin2 = 3000.0
real, parameter :: kh_hmmax2 = 10000.0
real, parameter :: kh_dimdepth = 1000.0
real, parameter :: km_strmin1 = 0.4
real, parameter :: km_strmax1 = 1.0
real, parameter :: km_strmin2 = 1.0
real, parameter :: km_strmax2 = 8.0
real, parameter :: km_hmmin1 = 300.0
real, parameter :: km_hmmax1 = 3000.0
real, parameter :: km_hmmin2 = 3000.0
real, parameter :: km_hmmax2 = 10000.0
real, parameter :: km_dimdepth = 1000.0 !1000.0
! legacy_kh_b
real, parameter :: kh_amnazmin = 1.0
real, parameter :: kh_amnazmax = 20.0
real, parameter :: kh_zzmin = 100000.0
real, parameter :: kh_zzmax = 150000.0
! legacy_km_b
real, parameter :: km_unstable = 500.0
!!! OCEVERMIX PARAMETERS VISIBILITY
public :: lu_min, km_unstable
private :: coriolis_rit_top_min, rit_max, s2_add, kh_ttmx, kh_ttmn, kh_unstmx, kh_unstmn
private :: hm_ice_max, hhq_add, kh_strmin1, kh_strmax1, kh_strmin2, kh_strmax2
private :: kh_hmmin1, kh_hmmax1, kh_hmmin2, kh_hmmax2, kh_dimdepth
private :: km_strmin1, km_strmax1, km_strmin2, km_strmax2
private :: km_hmmin1, km_hmmax1, km_hmmin2, km_hmmax2, km_dimdepth
private :: kh_amnazmin, kh_amnazmax, kh_zzmin, kh_zzmax
contains
subroutine legacy_rit(n2, s2, border_shift, lu, rit)
real, intent(in) :: n2(:,:,:)
real, intent(in) :: s2(:,:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: rit(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
! Parameters: rit_max, s2_add, lu_min
nx = size(rit, 1)
ny = size(rit, 2)
nz = size(rit, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
rit(i,j,k) = n2(i,j,k) / (s2(i,j,k) + s2_add)
rit(i,j,k) = min(rit(i,j,k), rit_max)
end if
end do
end do
end do
end subroutine legacy_rit
subroutine legacy_rit_top(rlh, tau_u, tau_v, border_shift, lu, rit_top)
real, intent(in) :: rlh(:,:)
real, intent(in) :: tau_u(:,:), tau_v(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: rit_top(:,:)
! Local variables
real :: rlt, coriolis, u_star ! Auxiliary variables
integer :: i, j ! Loop indices
integer :: nx, ny ! Array sizes
! Parameters: coriolis_rit_top_min, lu_min
nx = size(rit_top, 1)
ny = size(rit_top, 2)
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
rlt = (rlh(i,j) + rlh(i-1,j) + rlh(i,j-1) + rlh(i-1,j-1)) / 4.0
coriolis = max(abs(rlt), coriolis_rit_top_min) ! freezing in 10N,S
u_star = sqrt(sqrt(tau_u(i,j)**2 + tau_v(i,j)**2))
rit_top(i,j) = 0.5 * u_star / coriolis
end if
end do
end do
end subroutine legacy_rit_top
subroutine legacy_neutral_mld(rlh, tau_u, tau_v, border_shift, lu, neutral_mld, u_star)
real, intent(in) :: rlh(:,:)
real, intent(in) :: tau_u(:,:), tau_v(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: neutral_mld(:,:)
real, intent(out) :: u_star(:,:)
! Local variables
real :: rlt, coriolis ! Auxiliary variables
integer :: i, j ! Loop indices
integer :: nx, ny ! Array sizes
! Parameters: coriolis_rit_top_min, lu_min
nx = size(neutral_mld, 1)
ny = size(neutral_mld, 2)
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
rlt = (rlh(i,j) + rlh(i-1,j) + rlh(i,j-1) + rlh(i-1,j-1)) / 4.0
coriolis = max(abs(rlt), coriolis_rit_top_min) ! freezing in 10N,S
u_star(:,:) = sqrt(sqrt(tau_u(i,j)**2 + tau_v(i,j)**2))
neutral_mld(i,j) = 0.25 * u_star(i,j) / coriolis
end if
end do
end do
end subroutine legacy_neutral_mld
subroutine legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
real, intent(in) :: uu(:,:,:)
real, intent(in) :: dy(:,:)
real, intent(in) :: dyh(:,:)
real, intent(in) :: hhu(:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: u2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(uu, 1)
ny = size(uu, 2)
nz = size(uu, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
u2(i,j,k) = lu(i,j) * 0.5 / dy(i,j) * &
(uu(i-1,j,k) * hhu(i-1,j) * dyh(i-1,j) + &
uu(i,j,k) * hhu(i,j) * dyh(i,j)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_u2
subroutine legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
real, intent(in) :: vv(:,:,:)
real, intent(in) :: dx(:,:)
real, intent(in) :: dxh(:,:)
real, intent(in) :: hhv(:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: v2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(vv, 1)
ny = size(vv, 2)
nz = size(vv, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
v2(i,j,k) = lu(i,j) * 0.5 / dx(i,j) * &
(vv(i,j-1,k) * hhv(i,j-1) * dxh(i,j-1) + &
vv(i,j,k) * hhv(i,j) * dxh(i,j)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_v2
subroutine legacy_s2(u2, v2, lu, s2)
real, intent(in) :: u2(:,:,:)
real, intent(in) :: v2(:,:,:)
real, intent(in) :: lu(:,:)
real, intent(inout) :: s2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(s2, 1)
ny = size(s2, 2)
nz = size(s2, 3)
do k = 1, nz
do j = 1, ny
do i = 1, nx
if (lu(i,j) > lu_min) then
s2(i,j,k) = (u2(i,j,k+1) - u2(i,j,k))**2 + (v2(i,j,k+1) - v2(i,j,k))**2
end if
end do
end do
end do
end subroutine legacy_s2
subroutine legacy_n2(den, hhq, zw, g, lu, n2)
real, intent(in) :: den(:,:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
real, intent(in) :: zw(:)
real, intent(in) :: g
real, intent(inout) :: n2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(n2, 1)
ny = size(n2, 2)
nz = size(n2, 3)
do k = 1, nz
do j = 1, ny
do i = 1, nx
if (lu(i,j) > lu_min) then
n2(i,j,k) = g * (den(i,j,k+1) - den(i,j,k)) / 2.0 * &
(zw(k+2) - zw(k)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_n2
subroutine legacy_km_b(zw, unstable, undimdepth, km_b0, km_b)
real, intent(in) :: zw(:)
real, intent(in) :: unstable
real, intent(in) :: undimdepth
real, intent(in) :: km_b0
real, intent(inout) :: km_b(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(km_b) - 1
do k = 2, nz
! if (zw(k) <= undimdepth) then
! km_b(k) = unstable
! else
! km_b(k) = km_b0
! end if
km_b(k) = km_b0
end do
end subroutine legacy_km_b
subroutine legacy_km(rit, km_0, km_b, km)
real, intent(in) :: rit(:)
real, intent(in) :: km_0
real, intent(in) :: km_b(:)
real, intent(inout) :: km(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(rit)
do k = 2, nz
if (rit(k) > 0.0) then
km(k) = km_0 / (1.0 + 5.0 * rit(k))**2 + km_b(k)
else
km(k) = km_0 + km_b(k)
end if
end do
end subroutine legacy_km
subroutine legacy_kh(rit, kh_0, kh_b, unstable, kh)
real, intent(in) :: rit(:)
real, intent(in) :: kh_0
real, intent(in) :: unstable
real, intent(in) :: kh_b(:)
real, intent(inout) :: kh(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(rit)
do k = 2, nz
if (rit(k) >= 0.0) then
kh(k) = kh_0 / (1.0 + 5.0 * rit(k)) + kh_b(k)
else
kh(k) = unstable
end if
end do
! For unstable stratification
! do k = 2, nz
! if (rit(k) < 0.0) then
! Mixing is decreased in vicinity of equator in unstable situation
! kh(k) = unstable * (1.0 - 0.7 * (rn / rm) ** 10)
! kh(k) = unstable
! kh(k + 1) = unstable
! ! Convective diffusion as dZ * dZ * Vajsjala-Brenta frequency
! kh(k) = az0 * hhq * hzt(k) * &
! sqrt(abs(den(k) - den(k - 1)) * &
! g * hhq * hzt(k)) + &
! 0.3 * hhq**2 * hzt(k-1) * &
! hzt(k) / (12.0 * 3600.0)
! end if
! end do
end subroutine legacy_kh
subroutine legacy_kh_b(zw, hhq, unstable, undimdepth, kh_b0, kh_b)
real, intent(in) :: zw(:)
real, intent(in) :: hhq
real, intent(in) :: unstable
real, intent(in) :: undimdepth
real, intent(in) :: kh_b0
real, intent(inout) :: kh_b(:)
! Local variables
real :: amnaz, zz
integer :: k ! Loop indices
integer :: nz ! Array sizes
! Parameters: kh_amnazmin, kh_amnazmax, kh_zzmin, kh_zzmax
nz = size(kh_b) - 1
do k = 2, nz
if (zw(k) <= undimdepth) then
kh_b = unstable ! Addition mixing in upper layer for T & S
else
! BACKGROUND COEFFICIENT DEPENDS ON DEPTH
zz = zw(k) * hhq
if (zz <= kh_zzmin) then
amnaz = kh_amnazmin
else if (zz >= kh_zzmax) then
amnaz = kh_amnazmax
else
amnaz = kh_amnazmin + (kh_amnazmax - kh_amnazmin) * &
(zz - kh_zzmin) / (kh_zzmax - kh_zzmin)
end if
kh_b(k) = kh_b0 * amnaz
end if
end do
end subroutine legacy_kh_b
subroutine legacy_undimdepth(mode, str, hhq, aice0, undimdepth)
character(len=2) :: mode
real, intent(in) :: str
real, intent(in) :: hhq
real, intent(in) :: aice0
real, intent(inout) :: undimdepth
! Local variables
real :: strmin1, strmax1, strmin2, strmax2
real :: hmmin1, hmmax1, hmmin2, hmmax2
real :: dimdepth
real :: hm
! Parameters: kh_strmin1, kh_strmax1, kh_strmin2, kh_strmax2
! kh_hmmin1, kh_hmmax1, kh_hmmin2, kh_hmmax2
! kh_dimdepth
! km_strmin1, km_strmax1, km_strmin2, km_strmax2
! km_hmmin1, km_hmmax1, km_hmmin2, km_hmmax2
! km_dimdepth
if (mode == "kh") then
strmin1 = kh_strmin1
strmax1 = kh_strmax1
strmin2 = kh_strmin2
strmax2 = kh_strmax2
hmmin1 = kh_hmmin1
hmmax1 = kh_hmmax1
hmmin2 = kh_hmmin2
hmmax2 = kh_hmmax2
dimdepth = kh_dimdepth
else if (mode == "km") then
strmin1 = km_strmin1
strmax1 = km_strmax1
strmin2 = km_strmin2
strmax2 = km_strmax2
hmmin1 = km_hmmin1
hmmax1 = km_hmmax1
hmmin2 = km_hmmin2
hmmax2 = km_hmmax2
dimdepth = km_dimdepth
else
stop "legacy_undimdepth: ERROR: mode must be 'kh' or 'km'"
end if
if (str <= strmax1) then
hm = hmmin1 + (hmmax1 - hmmin1) * (str - strmin1) / (strmax1 - strmin1)
else
hm = hmmin2 + (hmmax2 - hmmin2) * (str - strmin2) / (strmax2 - strmin2)
end if
if (aice0 < hm_ice_max) then
hm = dimdepth
end if
undimdepth = amax1(dimdepth, hm) / (hhq + hhq_add)
end subroutine legacy_undimdepth
subroutine legacy_str(mode, taux, tauy, rh0, str)
character(len=2) :: mode
real, intent(in) :: taux, tauy
real, intent(in) :: rh0
real, intent(inout) :: str
! Local variables
real :: strmin, strmax
! Parameters: kh_strmin1, kh_strmax2, km_strmin1, km_strmax2
if (mode == "kh") then
strmin = kh_strmin1
strmax = kh_strmax2
else if (mode == "km") then
strmin = km_strmin1
strmax = km_strmax2
else
stop "legacy_str: ERROR: mode must be 'kh' or 'km'"
end if
str = rh0 * sqrt(taux**2 + tauy**2)
str = amax1(str, strmin) ! TODO: check is amax1, amin1 is not outdated
str = amin1(str, strmax)
end subroutine legacy_str
subroutine legacy_kh_unstable(tt, unstable)
real, intent(in) :: tt
real, intent(inout) :: unstable
! Parameters: kh_ttmx, kh_ttmn, kh_unstmx, kh_unstmn
if (tt > kh_ttmx) then
unstable = kh_unstmx
else if (tt < kh_ttmn) then
unstable = kh_unstmn
else
unstable = kh_unstmn + (kh_unstmx - kh_unstmn) * &
(tt - kh_ttmn) / (kh_ttmx - kh_ttmn)
end if
end subroutine legacy_kh_unstable
!> Computes the sea water density deviation from 1.02 [g/cm³].
!!
!! This function calculates the potential density deviation
!! of seawater as a function of potential temperature (**t** in °C),
!! salinity (**s** in parts per thousand), and pressure (**p** in MPa).
!! real, parameter :: prestompa = 1.0e-7 ! conversion factor for pressure [1MPa=1E-7din/cm**2] <init: Inc/0DENP.INC; using in: occont.f: ocforcing>
!!
!! **References:**
!! The calculation based on Bryden et al. (1999): "A new approximation of
!! the equation of state for seawater, suitable for numerical ocean models."
!! https://doi.org/10.1029/1998JC900059
!!
!! **Valid ranges:**
!! - Temperature: -2 < t < 40 °C
!! - Salinity: 0 < s < 42 ppt
!! - Pressure: 0 < p < 100 MPa
!!
!! **Input Variables:**
!! @param[in] t Potential temperature in degrees Celsius.
!! @param[in] s Salinity in parts per thousand (ppt).
!! @param[in] p Pressure in megapascals (MPa).
!! **Return Value:**
!! @return Potential density deviation from 1.02 [g/cm³].
elemental real function legacy_denp(t, s, p) result(denp)
real, intent(in) :: t, s, p
denp = -2.00920601e-02 + &
( 5.07043e-04 * p - 5.43283e-07 * p * p) + &
( 5.10768e-05 - 3.69119e-06 * p + 6.54837e-09 * p * p) * t + &
( 8.05999e-04 - 9.34012e-07 * p + 1.38777e-09 * p * p) * s + &
(-7.40849e-06 + 5.33243e-08 * p - 1.01563e-10 * p * p) * t * t + &
(-3.01036e-06 + 1.75145e-08 * p - 2.34892e-11 * p * p) * t * s + &
( 3.32267e-08 - 3.25887e-10 * p + 4.98612e-13 * p * p) * t * t * t + &
( 3.21931e-08 - 1.65849e-10 * p + 2.17612e-13 * p * p) * t * t * s
end function legacy_denp
end module obl_legacy