porous_open_width Function

   pure function porous_open_width(d_min, d_max, d_avg, eta) result(w)
      !! Open WIDTH fraction of a face at interface height `eta`
      !! (dimensionless, in `[0, 1]`).  Zero when the interface is at or
      !! below the deepest along-face point, one when it is above the
      !! shallowest.  A degenerate face (`d_max <= d_min`, i.e. a flat
      !! along-face seafloor) reduces to the binary open/closed step, so
      !! the fit never divides by zero.
      real(wp), intent(in) :: d_min
         !! Deepest along-face topographic height (m, positive up).
      real(wp), intent(in) :: d_max
         !! Shallowest along-face topographic height (m, positive up).
      real(wp), intent(in) :: d_avg
         !! Mean along-face topographic height (m, positive up).
      real(wp), intent(in) :: eta
         !! Interface height at the face (m, positive up).
      real(wp) :: w
      real(wp) :: m, a, zeta, drange

      if (eta <= d_min) then
         w = 0.0_wp
      else if (eta > d_max) then
         w = 1.0_wp
      else if (d_avg <= d_min) then
         ! Degenerate m = 0: the profile collapses to a step at d_min, so
         ! everything above it is open.  Guarded explicitly because
         ! `a = (1-m)/m` would divide by zero here.
         w = 1.0_wp
      else if (d_avg >= d_max) then
         ! Degenerate m = 1: a step at d_max — nothing below it is open.
         w = 0.0_wp
      else
         drange = d_max - d_min
         m = (d_avg - d_min)/drange
         zeta = (eta - d_min)/drange
         if (m < 0.5_wp) then
            ! a = (1-m)/m, so 1/a = m/(1-m) = (d_avg-d_min)/(d_max-d_avg).
            a = (1.0_wp - m)/m
            w = zeta**(1.0_wp/a)
         else if (m > 0.5_wp) then
            a = (1.0_wp - m)/m
            w = 1.0_wp - (1.0_wp - zeta)**a
         else
            w = zeta
         end if
      end if
   end function porous_open_width