! ============================================================================= !> This submodule defines an equilibrium containing magnetothermal instabilities !! in a cylindrical geometry. The geometry can be overridden in the parfile. !! !! This equilibrium is taken from !! _Van der Linden, R. A. M., Goossens, M., & Hood, A. W. (1992). !! The relevance of the ballooning approximation for magnetic, thermal, !! and coalesced magnetothermal instabilities. Solar physics, 140(2), 317-342_. !! @note Default values are given by !! !! - k2 = 0 !! - k3 = 1 !! - cte_T0 = 1 : used to set the temperature. !! - cooling_curve = 'rosner' !! - parallel thermal conduction, no perpendicular conduction !! !! and normalisations given by !! !! - unit_temperature = 2.6e6 K !! - unit_magneticfield = 10 Gauss !! - unit_length = 1e8 cm !! !! and can all be changed in the parfile. @endnote !! @note The default setup handles _CASE I_ in the original paper. For _CASE II_, !! you can set k2 = 10. To reproduce the thermal continuum plot, !! set unit_length = 2.44e9 cm in _CASE I_. @endnote submodule(mod_equilibrium) smod_equil_magnetothermal_instabilities use mod_equilibrium_params, only: cte_T0 implicit none contains !> Sets the equilibrium module procedure magnetothermal_instability_eq if (settings%equilibrium%use_defaults) then ! LCOV_EXCL_START call settings%grid%set_geometry("cylindrical") call settings%grid%set_grid_boundaries(0.0_dp, 1.0_dp) call settings%physics%enable_cooling(cooling_curve="rosner") call settings%physics%enable_heating(force_thermal_balance=.true.) call settings%physics%enable_parallel_conduction() call settings%units%set_units_from_temperature( & unit_temperature=2.6e6_dp, & unit_magneticfield=10.0_dp, & unit_length=1.00e8_dp, & mean_molecular_weight=1.0_dp & ! this work assumes pure proton plasma ) cte_T0 = 1.0_dp k2 = 0.0_dp k3 = 1.0_dp end if ! LCOV_EXCL_STOP call background%set_density_funcs(rho0_func=rho0, drho0_func=drho0) call background%set_temperature_funcs(T0_func=T0) call background%set_magnetic_2_funcs(B02_func=B02, dB02_func=dB02) end procedure magnetothermal_instability_eq real(dp) function p0(r) real(dp), intent(in) :: r p0 = 1.0_dp / (2.0_dp * (1.0_dp + r**2)**2) end function p0 real(dp) function rho0(r) real(dp), intent(in) :: r rho0 = p0(r) / cte_T0 end function rho0 real(dp) function drho0(r) real(dp), intent(in) :: r drho0 = -2.0_dp * r / (cte_T0 * (r**2 + 1.0_dp)**3) end function drho0 real(dp) function T0() T0 = cte_T0 end function T0 real(dp) function B02(r) real(dp), intent(in) :: r B02 = r / (1.0_dp + r**2) end function B02 real(dp) function dB02(r) real(dp), intent(in) :: r dB02 = (1.0_dp - r**2) / (r**4 + 2.0_dp * r**2 + 1.0_dp) end function dB02 end submodule smod_equil_magnetothermal_instabilities