sample_monotonic_phi_hat_at_z Subroutine

  subroutine sample_monotonic_phi_hat_at_z(p, branch, phi0_hat, phi_m_hat, n_swe_inf_hat, z_m, phi_hat, side, z_m_hat)
    type(zhao_params_type), intent(in) :: p
    character(len=1), intent(in) :: branch
    real(dp), intent(in) :: phi0_hat, phi_m_hat, n_swe_inf_hat, z_m
    real(dp), intent(out) :: phi_hat, z_m_hat
    character(len=16), intent(out) :: side

    integer :: i, ngrid
    real(dp) :: z_hat_target, x, phi_end_hat, dphi, integral
    real(dp), allocatable :: phi_nodes(:), rho_nodes(:), e2_nodes(:), z_nodes(:), cumulative(:)

    z_hat_target = z_m/p%lambda_d_phe_ref_m
    z_m_hat = 0.0d0
    side = 'monotonic'
    ngrid = max(2000, zhao_profile_grid)
    allocate (phi_nodes(ngrid), rho_nodes(ngrid), e2_nodes(ngrid), z_nodes(ngrid))

    select case (branch)
    case ('B')
      phi_end_hat = abs(zhao_profile_phi_tol_hat)
      do i = 1, ngrid
        x = real(i - 1, dp)/real(ngrid - 1, dp)
        phi_nodes(i) = phi0_hat + (phi_end_hat - phi0_hat)*(2.0d0*x - x*x)
      end do
      do i = 1, ngrid
        call evaluate_zhao_rho_hat(p, branch, side, phi_nodes(i), phi0_hat, phi_m_hat, n_swe_inf_hat, rho_nodes(i))
      end do
      e2_nodes(ngrid) = 0.0d0
      integral = 0.0d0
      do i = ngrid - 1, 1, -1
        dphi = phi_nodes(i + 1) - phi_nodes(i)
        integral = integral + 0.5d0*(rho_nodes(i) + rho_nodes(i + 1))*dphi
        e2_nodes(i) = max(2.0d0*integral, 0.0d0)
      end do
    case ('C')
      phi_end_hat = -abs(zhao_profile_phi_tol_hat)
      allocate (cumulative(ngrid))
      do i = 1, ngrid
        x = real(i - 1, dp)/real(ngrid - 1, dp)
        phi_nodes(i) = phi0_hat + (phi_end_hat - phi0_hat)*(2.0d0*x - x*x)
      end do
      do i = 1, ngrid
        call evaluate_zhao_rho_hat(p, branch, side, phi_nodes(i), phi0_hat, phi_m_hat, n_swe_inf_hat, rho_nodes(i))
      end do
      cumulative(1) = 0.0d0
      do i = 2, ngrid
        dphi = phi_nodes(i) - phi_nodes(i - 1)
        cumulative(i) = cumulative(i - 1) + 0.5d0*(rho_nodes(i - 1) + rho_nodes(i))*dphi
      end do
      do i = 1, ngrid
        e2_nodes(i) = max(2.0d0*(cumulative(ngrid) - cumulative(i)), 0.0d0)
      end do
      deallocate (cumulative)
    case default
      error stop 'Unknown monotonic Zhao branch.'
    end select

    z_nodes(1) = 0.0d0
    do i = 2, ngrid
      dphi = abs(phi_nodes(i) - phi_nodes(i - 1))
      z_nodes(i) = z_nodes(i - 1) + dphi/sqrt(max(0.5d0*(e2_nodes(i - 1) + e2_nodes(i)), 1.0d-14))
    end do

    if (z_hat_target <= z_nodes(ngrid)) then
      phi_hat = interpolate_profile_value(z_nodes, phi_nodes, z_hat_target)
    else
      phi_hat = 0.0d0
    end if
  end subroutine sample_monotonic_phi_hat_at_z