diff --git a/obl_inmom.f90 b/obl_inmom.f90
new file mode 100644
index 0000000000000000000000000000000000000000..beb63413bc3da3d1f5b92520ba10d2eb74d0d5da
--- /dev/null
+++ b/obl_inmom.f90
@@ -0,0 +1,161 @@
+module obl_inmom
+  implicit none
+  integer, parameter :: border_shift = 1
+  private :: border_shift
+
+  contains
+
+  subroutine init_vermix(richnum_mode, kh_km_mode, uu, vv, lu, dx, dy, dxh, dyh, hhu, hhv, hhq, zw, g, rh0, &
+    rit, den, rlh, taux, tauy, aice0, tt, ss, anzt, anzu, anumaxt, anumaxu, anubgrt, anubgru)
+    use obl_legacy
+    use obl_pph
+    use obl_pph_dyn
+    integer, intent(in) :: richnum_mode, kh_km_mode
+    real, intent(inout) :: anzt(:,:,:) ! kh
+    real, intent(inout) :: anzu(:,:,:) ! km
+    real, intent(inout) :: rit(:,:,:)
+    real, intent(inout) :: den(:,:,:)
+    real, intent(in) :: uu(:,:,:), vv(:,:,:)
+    real, intent(in) :: tt(:,:,:)
+    real, intent(in) :: ss(:,:,:)
+    real, intent(in) :: taux(:,:), tauy(:,:)
+    real, intent(in) :: lu(:,:)
+    real, intent(in) :: dx(:,:), dy(:,:), dxh(:,:), dyh(:,:)
+    real, intent(in) :: hhu(:,:), hhv(:,:), hhq(:,:)
+    real, intent(in) :: aice0(:,:)
+    real, intent(in) :: rlh(:,:)
+    real, intent(in) :: zw(:)
+    real, intent(in) :: g, rh0
+    real, intent(in) :: anumaxt ! kh_0
+    real, intent(in) :: anumaxu ! km_0
+    real, intent(in) :: anubgrt ! kh_b0
+    real, intent(in) :: anubgru ! km_b0
+    ! Local variables
+    real, allocatable :: u2(:,:,:), v2(:,:,:)
+    real, allocatable :: s2(:,:,:)
+    real, allocatable :: n2(:,:,:)
+    real, allocatable :: kh_b(:)
+    real, allocatable :: km_b(:)
+    real, allocatable :: kh(:,:,:)
+    real, allocatable :: km(:,:,:)
+    real :: kh_str, km_str
+    real :: kh_undimdepth, km_undimdepth
+    real :: kh_unstable
+    real :: kh_b0, km_b0
+    real :: kh_0, km_0
+    integer :: nx, ny, nz
+    integer :: i, j
+    ! Modules types
+    type(pphParamType) :: pphParams
+    type(gridDataType) :: inmomGrid
+    
+    allocate(u2(size(uu, 1), size(uu, 2), size(uu, 3)))
+    allocate(v2(size(vv, 1), size(vv, 2), size(vv, 3)))
+    allocate(s2(size(u2, 1), size(u2, 2), size(u2, 3) - 1))
+    allocate(n2(size(den, 1), size(den, 2), size(den, 3) - 1))
+    allocate(kh(size(anzt, 1), size(anzt, 2), size(anzt, 3)))
+    allocate(km(size(anzu, 1), size(anzu, 2), size(anzu, 3)))
+    allocate(kh_b(size(kh, 3)))
+    allocate(km_b(size(km, 3)))
+
+    kh = anzt
+    km = anzu
+    kh_0 = anumaxt
+    km_0 = anumaxu
+    kh_b0 = anubgrt
+    km_b0 = anubgru
+
+    nx = size(rit, 1)
+    ny = size(rit, 2)
+    nz = size(rit, 3)
+
+    ! print *, "hello 1"
+    den = legacy_denp(tt, ss + 35.0, 0.0)
+
+    if (richnum_mode == 1) then
+      ! print *, "hello rit calc"
+      call legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
+      call legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
+      call legacy_s2(u2, v2, lu, s2)
+      call legacy_n2(den, hhq, zw, g, lu, n2)
+      call legacy_rit(n2, s2, border_shift, lu, rit(:,:,2:size(rit, 3)))
+      call legacy_rit_top(rlh, taux, tauy, border_shift, lu, rit(:,:,1))
+    end if
+    ! print *, "hello aft rit"
+    call sync_xy_border(rit)
+
+    do j = 1, ny
+      do i = 1, nx
+        if (lu(i, j) > lu_min) then
+          if (kh_km_mode == 1) then
+            call legacy_str("kh", taux(i,j), tauy(i,j), rh0, kh_str)
+            call legacy_undimdepth("kh", kh_str, hhq(i,j), aice0(i,j), kh_undimdepth)
+            call legacy_kh_b(zw(:), hhq(i,j), kh_unstable, kh_undimdepth, kh_b0, kh_b)
+            call legacy_kh_unstable(tt(i,j,1), kh_unstable)
+            call legacy_kh(rit(i,j,:), kh_0, kh_b(:), kh_unstable, kh(i,j,:))
+
+            call legacy_str("km", taux(i,j), tauy(i,j), rh0, km_str)
+            call legacy_undimdepth("km", km_str, hhq(i,j), aice0(i,j), km_undimdepth)
+            call legacy_km_b(zw, km_unstable, km_undimdepth, km_b0, km_b)
+            call legacy_km(rit(i,j,:), km_0, km_b(:), km(i,j,:))
+          else if (kh_km_mode == 2) then
+            ! call get_Kh(kh(i,j,:), rit(i,j,:), nz)
+            ! call get_Km(km(i,j,:), rit(i,j,:), nz)
+            call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), pphParams, inmomGrid)
+            call get_eddy_viscosity(km(i,j,:), rit(i,j,:), pphParams, inmomGrid)
+          else if (kh_km_mode == 3) then
+            ! u_dynH = 
+            ! mld = 
+            call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
+            call get_eddy_viscosity(km(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
+          end if
+        end if
+      end do
+    end do
+
+    if (kh_km_mode == 2) then
+      kh = kh * 10000.0
+      km = km * 10000.0
+    end if
+    
+    ! print *, "Kh first:", kh(3,3,1:4)
+    ! print *, "Km first:", km(3,3,1:4)
+
+    ! print *, "Kh last:", kh(3,3,size(kh,3)-2:size(kh,3))
+    ! print *, "Km last:", km(3,3,size(km,3)-2:size(km,3))
+    anzt = kh
+    anzu = km
+  end subroutine init_vermix
+
+  subroutine sync_xy_border(array)
+    implicit none
+    real, intent(inout) :: array(:,:,:)
+    integer :: nx, ny, nz, k
+  
+    ! Determine array dimensions
+    nx = size(array, 1)
+    ny = size(array, 2)
+    nz = size(array, 3)
+
+    ! Update boundary points along x (first_x and end_x boundaries)
+    do k = 1, nz
+      array(1,2:ny-1,k) = array(2,2:ny-1,k)      ! first_x boundary
+      array(nx,2:ny-1,k) = array(nx-1,2:ny-1,k)  ! end_x boundary
+    end do
+  
+    ! Update boundary points along y (first_y and end_y boundaries)
+    do k = 1, nz
+      array(2:nx-1,1,k) = array(2:nx-1,2,k)      ! first_y boundary
+      array(2:nx-1,ny,k) = array(2:nx-1,ny-1,k)  ! end_y boundary
+    end do
+  
+    ! Update corner points
+    do k = 1, nz
+      array(1,1,k) = (array(2,1,k) + array(1,2,k)) / 2.0
+      array(1,ny,k) = (array(2,ny,k) + array(1,ny-1,k)) / 2.0
+      array(nx,1,k) = (array(nx-1,1,k) + array(nx,2,k)) / 2.0
+      array(nx,ny,k) = (array(nx-1,ny,k) + array(nx,ny-1,k)) / 2.0
+    end do
+  
+  end subroutine sync_xy_border
+end module obl_inmom
\ No newline at end of file
diff --git a/obl_legacy.f90 b/obl_legacy.f90
new file mode 100644
index 0000000000000000000000000000000000000000..88334706a175bd029567ad66973e49720672a225
--- /dev/null
+++ b/obl_legacy.f90
@@ -0,0 +1,465 @@
+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
\ No newline at end of file