meke_length_scales Subroutine

   pure subroutine meke_length_scales(nx, ny, cd_scale, cb, ct, min_gamma2, &
                                      cdrag, a_deform, a_rhines, &
                                      a_eady, a_frict, a_grid, &
                                      areaT, idxT, idyT, f_centre, &
                                      depth_tot, meke, &
                                      rd_over_dx, sn_u, sn_v, &
                                      bottom_fac2, barotr_fac2, le)
      !! Fill the structure factors gamma_b^2 (`bottom_fac2`) and gamma_t^2
      !! (`barotr_fac2`) plus the mixing length `le` (Lmix) at each cell
      !! centre.  `Ldeform/Lfrict` drives both gammas; `Lmix` is the harmonic
      !! sum of the alpha-weighted scales.  `beta = |grad f|` from centred
      !! `f_centre` differences scaled by `idxT`/`idyT` (zero when f_centre
      !! is unfilled ⇒ Rhines inert).  SN = 0.25*(sn_u(i)+sn_u(i-1)+
      !! sn_v(j)+sn_v(j-1)) only when aEady>0.
      integer, intent(in) :: nx, ny
      real(wp), intent(in) :: cd_scale, cb, ct, min_gamma2, cdrag
      real(wp), intent(in) :: a_deform, a_rhines, a_eady, a_frict, a_grid
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: idxT(nx, ny)
      real(wp), intent(in) :: idyT(nx, ny)
      real(wp), intent(in) :: f_centre(nx, ny)
      real(wp), intent(in) :: depth_tot(nx, ny)
      real(wp), intent(in) :: meke(nx, ny)
      real(wp), intent(in) :: rd_over_dx(nx, ny)
      real(wp), intent(in) :: sn_u(nx + 1, ny)
      real(wp), intent(in) :: sn_v(nx, ny + 1)
      real(wp), intent(out) :: bottom_fac2(nx, ny)
      real(wp), intent(out) :: barotr_fac2(nx, ny)
      real(wp), intent(out) :: le(nx, ny)
      integer :: i, j
      real(wp) :: lgrid, ldeform, lfrict, ratio, bf2, tf2
      real(wp) :: ueddy, sn, beta, inv_l, rd

      do concurrent(j=1:ny, i=1:nx) &
         local(lgrid, ldeform, lfrict, ratio, bf2, tf2, ueddy, sn, beta, inv_l, rd)
         rd = rd_over_dx(i, j)
         lgrid = sqrt(max(areaT(i, j), 0.0_wp))
         ldeform = lgrid*rd
         lfrict = 0.0_wp
         if (cdrag > 0.0_wp) lfrict = depth_tot(i, j)/cdrag

         ! gamma_b^2 = cd_scale^2 + 1/(1+Cb*Ldeform/Lfrict)^0.8 (floor).
         bf2 = cd_scale*cd_scale
         if (lfrict*cb > 0.0_wp) then
            ratio = ldeform/lfrict
            bf2 = bf2 + 1.0_wp/(1.0_wp + cb*ratio)**0.8_wp
         end if
         bf2 = max(bf2, min_gamma2)
         bottom_fac2(i, j) = bf2

         ! gamma_t^2 = 1/(1+Ct*Ldeform/Lfrict)^0.25 (floor).
         tf2 = 1.0_wp
         if (lfrict*ct > 0.0_wp) then
            ratio = ldeform/lfrict
            tf2 = 1.0_wp/(1.0_wp + ct*ratio)**0.25_wp
         end if
         tf2 = max(tf2, min_gamma2)
         barotr_fac2(i, j) = tf2

         ! Mixing length: harmonic sum of alpha-weighted scales.
         ueddy = sqrt(2.0_wp*max(0.0_wp, tf2*meke(i, j)))
         ! beta = |grad f|; centred f_centre differences scaled to a true
         ! gradient by the cell-centre inverse metrics (df/dx ~ (f_{i+1} -
         ! f_{i-1})*idxT/2).  When f_centre is unfilled (default) beta=0 ⇒
         ! Rhines weight inert (alpha_rhines default 0).
         beta = 0.0_wp
         if (i > 1 .and. i < nx) then
            beta = beta + (0.5_wp*(f_centre(i + 1, j) - f_centre(i - 1, j))*idxT(i, j))**2
         end if
         if (j > 1 .and. j < ny) then
            beta = beta + (0.5_wp*(f_centre(i, j + 1) - f_centre(i, j - 1))*idyT(i, j))**2
         end if
         beta = sqrt(beta)
         sn = 0.0_wp
         if (a_eady > 0.0_wp) sn = 0.25_wp*((sn_u(i, j) + sn_u(i + 1, j)) + &
                                            (sn_v(i, j) + sn_v(i, j + 1)))

         inv_l = meke_inv_lmix(ueddy, sn, beta, areaT(i, j), rd, &
                               depth_tot(i, j), cdrag, &
                               a_deform, a_rhines, a_eady, a_frict, a_grid)
         if (inv_l > 0.0_wp) then
            le(i, j) = 1.0_wp/inv_l
         else
            le(i, j) = 0.0_wp
         end if
      end do
   end subroutine meke_length_scales