zhao_residuals_type_a Subroutine

public subroutine zhao_residuals_type_a(p, x, f)

Arguments

Type IntentOptional Attributes Name
type(zhao_params_type), intent(in) :: p
real(kind=dp), intent(in) :: x(:)
real(kind=dp), intent(out) :: f(:)

Calls

proc~~zhao_residuals_type_a~~CallsGraph proc~zhao_residuals_type_a zhao_residuals_type_a proc~swe_free_current_term swe_free_current_term proc~zhao_residuals_type_a->proc~swe_free_current_term proc~type_a_e2_sum_at_infinity type_a_e2_sum_at_infinity proc~zhao_residuals_type_a->proc~type_a_e2_sum_at_infinity

Source Code

  subroutine zhao_residuals_type_a(p, x, f)
    type(zhao_params_type), intent(in) :: p
    real(dp), intent(in) :: x(:)
    real(dp), intent(out) :: f(:)

    real(dp) :: phi0_v, phi_m_v, n_swe_inf_m3, a_swe, a_phe, ion_term

    phi0_v = x(1)
    phi_m_v = x(2)
    n_swe_inf_m3 = x(3)
    if (phi0_v <= 0.0d0 .or. phi_m_v >= 0.0d0 .or. phi_m_v >= phi0_v .or. n_swe_inf_m3 <= 0.0d0) then
      f = 1.0d6
      return
    end if

    a_swe = sqrt(max(0.0d0, -phi_m_v/p%t_swe_ev)) - p%u
    a_phe = sqrt(max(0.0d0, -phi_m_v/p%t_phe_ev))
    ion_term = p%n_swi_inf_m3*sqrt(2.0d0*pi*p%t_swe_ev/p%t_phe_ev*p%m_e_kg/p%m_i_kg)*p%mach

    f(1) = 0.5d0*n_swe_inf_m3*(1.0d0 + 2.0d0*erf(p%u) + erf(a_swe)) + &
           0.5d0*p%n_phe0_m3*exp(-phi0_v/p%t_phe_ev)*(1.0d0 - erf(a_phe)) - p%n_swi_inf_m3
    f(2) = p%n_phe0_m3*exp((phi_m_v - phi0_v)/p%t_phe_ev) - swe_free_current_term(p, n_swe_inf_m3, a_swe) + ion_term
    f(3) = type_a_e2_sum_at_infinity(p, phi0_v, phi_m_v, n_swe_inf_m3)
  end subroutine zhao_residuals_type_a