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_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