add legacy_neutral_mld subroutine

obl_legacy.f90

@@ -106,6 +106,34 @@ module obl_legacy
     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(:,:)