! $Id$ ! ! This modules contains the old routines for SNe-driven ISM simulations. ! It may be necessary to select this for using old simulation data ! !***************** AUTOMATIC CPARAM.INC GENERATION *************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: linterstellar = .true. ! ! MAUX CONTRIBUTION 2 ! COMMUNICATED AUXILIARIES 1 ! !***************************************************************************** module Interstellar ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! include 'interstellar.h' include 'record_types.h' ! type ExplosionSite real :: rho, lnrho, yH, lnTT, TT, ss, ee endtype ! type SNRemnant real :: x, y, z, t ! Time and location real :: EE, MM ! Mass and energy injected real :: rhom ! Local mean density at explosion time real :: radius ! Injection radius real :: t_sedov real :: t_damping real :: heat_energy real :: damping_factor real :: energy_loss integer :: l,m,n ! Grid position integer :: iproc,ipx,ipy,ipz integer :: SN_type integer :: state type (ExplosionSite) :: site endtype ! ! Enumeration of Supernovae types. ! integer, parameter :: SNtype_I = 1, SNtype_II = 2 ! ! Enumeration of Supernovae states. ! integer, parameter :: SNstate_invalid = 0 integer, parameter :: SNstate_waiting = 1 integer, parameter :: SNstate_exploding = 2 integer, parameter :: SNstate_damping = 3 integer, parameter :: SNstate_finished = 4 ! ! Enumeration of Explosion Errors ! integer, parameter :: iEXPLOSION_OK = 0 integer, parameter :: iEXPLOSION_TOO_HOT = 1 integer, parameter :: iEXPLOSION_TOO_RARIFIED = 2 integer, parameter :: iEXPLOSION_TOO_UNEVEN = 3 ! ! 04-sep-09/fred: amended xsi_sedov ! ref Dyson & Williams Ch7 value = (25/3/pi)**(1/5)=1.215440704 for gamma=5/3 ! 2nd ref Ostriker & McKee 1988 Rev.Mod.Phys 60,1 ! Est'd value for similarity variable at shock ! ! real :: xsi_sedov=1.215440704 real :: xsi_sedov=1.15166956 ! ! 'Current' SN Explosion site parameters ! integer, parameter :: mSNR = 120 integer :: nSNR = 0 type (SNRemnant), dimension(mSNR) :: SNRs integer, dimension(mSNR) :: SNR_index integer, parameter :: npreSN = 5 integer, dimension(4,npreSN) :: preSN ! ! integer :: icooling=0 ! integer :: inetheat=0 ! ! Squared distance to the SNe site along the current pencil ! Outward normal vector from SNe site along the current pencil ! real, dimension(nx) :: dr2_SN real, dimension(nx,3) :: outward_normal_SN ! ! Allocate time of next SNI/II and intervals until next ! real :: t_next_SNI=0.0, t_next_SNII=0.0, t_next_mass=0.0 real :: x_next_SNII=0.0, y_next_SNII=0.0 real :: t_interval_SNI=impossible, t_interval_SNII=impossible ! ! normalisation factors for 1-d, 2-d, and 3-d profiles like exp(-r^6) ! ( 1d: 2 int_0^infty exp(-(r/a)^6) dr) / a ! 2d: 2 pi int_0^infty exp(-(r/a)^6) r dr) / a^2 ! 3d: 4 pi int_0^infty exp(-(r/a)^6) r^2 dr) / a^3 ) ! ( cf. 3.128289613 -- from where ?!? ) ! NB: 1d and 2d results just from numerical integration -- calculate ! exact integrals at some point... ! 3-D was 3.71213666 but replaced with Maple result.... ! real, parameter, dimension(3) :: & cnorm_gaussian_SN = (/ 0.8862269255, 3.141592654, 5.568327998 /) real, parameter, dimension(3) :: & cnorm_gaussian2_SN = (/ 0.9064024771, 2.784163999, 3.849760109 /) real, parameter, dimension(3) :: & cnorm_SN = (/ 1.855438667 , 2.805377875 , 3.712218666 /) real, parameter, dimension(3) :: & cnorm_para_SN = (/ 1.33333333, 1.5707963, 1.6755161 /) real, parameter, dimension(3) :: & cnorm_quar_SN = (/ 0., 2.0943951, 0. /) ! ! cp1=1/cp used to convert TT (and ss) into interstellar code units ! (useful, as many conditions conveniently expressed in terms of TT) ! code units based on: ! [length] = 1kpc = 3.09 10^21 cm ! [time] = 1Gyr = 3.15 10^16 s !no, on [u]=1km/s... ! [rho] = = 1.00 10^-24 g/cm^3 ! Lambdaunits converts coolH into interstellar code units. ! real :: unit_Lambda ! ! Minimum resulting central temperature of a SN explosion. ! If this is not reached then consider moving mass to achieve this. ! real, parameter :: TT_SN_min_cgs=1.E6 ! ! 22-jan-10/fred: ! With lSN_velocity kinetic energy lower limit no longer required for shock ! speed. ! 10-aug-10/fred: ! As per joung et al apj653 2005 min temp 1E6 to avoid excess radiative ! energy losses in early stages. ! real :: uu_sedov_max=0. real :: TT_SN_min=impossible real :: TT_cutoff_cgs=100. real :: TT_cutoff=impossible, TT_cutoff1=impossible logical :: lTT_cutoff=.false. ! ! SNe placement limitations (for code stability) ! real, parameter :: rho_SN_min_cgs=1E-28,rho_SN_max_cgs=5E-24 real, parameter :: TT_SN_max_cgs=5E9 real :: rho_SN_min=impossible, TT_SN_max=impossible, rho_SN_max=impossible ! ! SNI per (x,y)-area explosion rate ! double precision, parameter :: SNI_area_rate_cgs=1.330982784D-56 real :: SNI_area_rate=impossible, SNII_area_rate=impossible ! ! SNII rate=5.E-12 mass(H1+HII)/solar_mass ! van den Bergh/Tammann Annu. Rev Astron. Astrophys. 1991 29:363-407 ! SNI rate=4.7E-14/solar_mass + 0.35 x SNII rate ! Mannucci et al A&A 433, 807-814 (2005) ! real, parameter :: SNII_mass_rate_cgs=1.584434515E-19 real, parameter :: SNI_mass_rate_cgs=1.489368444E-21 real :: SNII_mass_rate, SNI_mass_rate logical :: lSN_mass_rate=.false. ! ! Some useful constants ! real, parameter :: kpc_cgs=3.086E+21 ! [cm] real, parameter :: yr_cgs=3.155692E7 ! [s] real, parameter :: solar_mass_cgs=1.989E33! [g] real :: solar_mass=impossible ! ! Scale heights for SNI/II with Gaussian z distributions ! real, parameter :: h_SNI_cgs=1.00295E21, h_SNII_cgs=2.7774E20 real :: h_SNI=impossible, h_SNII=impossible ! ! Self regulating SNII explosion coefficients ! real, parameter :: cloud_rho_cgs=1.67262158E-24, cloud_TT_cgs=4000. real, parameter :: cloud_tau_cgs=2.E7 * yr_cgs, minTT_cgs = 0.75E2 real, parameter :: mass_SN_progenitor_cgs=10.*solar_mass_cgs real, parameter :: frac_converted=0.02, frac_heavy=0.10 ! real, parameter :: tosolarMkpc3=1.483E7 real :: cloud_rho=impossible, cloud_TT=impossible real :: cloud_tau=impossible real :: mass_SN_progenitor=impossible ! ! Total SNe energy ! double precision, parameter :: ampl_SN_cgs=1D51 real :: frac_ecr=0.1, frac_eth=0.9 real :: ampl_SN=impossible, kampl_SN=impossible ! ! SNe composition ! logical :: lSN_eth=.true., lSN_ecr=.true., lSN_mass=.true., & lSN_velocity=.false., lSN_fcr=.false. ! ! Total mass added by a SNe ! real, parameter :: mass_SN_cgs=10.*solar_mass_cgs real :: mass_SN=impossible real :: velocity_SN=impossible ! ! Size of SN insertion site (energy and mass) and shell in mass movement ! real :: sigma_SN, sigma_SN1 real, parameter :: width_SN_cgs=3.086E19 real :: energy_width_ratio=1. real :: mass_width_ratio=1. real :: velocity_width_ratio=1. real :: outer_shell_proportion = 1.2 real :: inner_shell_proportion = 1. real :: width_SN=impossible ! ! Parameters for 'averaged'-SN heating ! real :: r_SNI_yrkpc2=4.E-6, r_SNII_yrkpc2=3.E-5 real :: r_SNI=3.E+4, r_SNII=4.E+3 real :: average_SNI_heating=0., average_SNII_heating=0. ! ! Limit placed of minimum density resulting from cavity creation and ! parameters for thermal_hse(hydrostatic equilibrium) assuming RBNr ! real, parameter :: rho_min_cgs=1.E-34, rho0ts_cgs=3.5E-24, T0hs_cgs=7.088E2 real :: rho0ts=impossible, T0hs=impossible, rho_min=impossible ! ! Cooling timestep limiter coefficient ! (This value 0.08 is overly restrictive. cdt_tauc=0.5 is a better value.) ! real :: cdt_tauc=0.5 ! ! Time of most recent SNII event ! real :: last_SN_t=0. ! ! Time to wait before allowing SN to begin firing ! real :: t_settle=0. ! ! Initial dist'n of explosions ! character (len=labellen), dimension(ninit) :: initinterstellar='nothing' ! ! Number of randomly placed SNe to generate to mix the initial medium ! integer :: initial_SNI=0 ! ! Parameters for UV heating of Wolfire et al. ! real, parameter :: rhoUV_cgs=0.1 real, parameter :: GammaUV_cgs=0.0147 real, parameter :: TUV_cgs=7000., T0UV_cgs=20000., cUV_cgs=5.E-4 real :: GammaUV=impossible, T0UV=impossible, cUV=impossible ! ! 04-jan-10/fred: ! Amended cool dim from 7 to 11 to accomodate WSW dimension. ! Appended null last term to all arrays for RBN and SS cooling ! double precision, dimension(11) :: coolH_cgs real, dimension(11) :: coolT_cgs real, dimension(11) :: coolB, lncoolH, lncoolT integer :: ncool ! ! TT & z-dependent uv-heating profile ! real, dimension(mz) :: heat_z, zrho logical :: lthermal_hse=.false., lheatz_min=.true. ! real :: coolingfunction_scalefactor=1. real :: heatingfunction_scalefactor=1. ! real :: heating_rate = 0.015 real :: heating_rate_code = impossible ! real :: heatcool_shock_cutoff = 0. real :: heatcool_shock_cutoff_rate = 0. real :: heatcool_shock_cutoff_rate1 = 0.0 ! real :: cooltime_despike_factor = 2. ! ! Set .true. to smooth the radiative cooling in cooling time space. ! logical :: lcooltime_despike = .false. logical :: lcooltime_smooth = .false. ! ! Set .true. to smoothly turn off the heating and cooling where the ! shock_profile is > heatcool_shock_cutoff ! logical :: lheatcool_shock_cutoff = .false. ! ! SN type flags ! logical :: lSNI=.true., lSNII=.false., lSNII_gaussian=.true. ! ! Damp central regions of SNRs at early times ! logical :: lSNR_damping = .false. real :: SNR_damping = 0. real :: SNR_damping_time_cgs = 5E4*yr_cgs real :: SNR_damping_rate_cgs = 1E3*yr_cgs real :: SNR_damping_time = impossible real :: SNR_damping_rate = impossible ! ! Cooling & heating flags ! logical :: lsmooth_coolingfunc = .false. logical :: laverage_SN_heating = .false. logical :: lheating_UV = .true. ! ! Remnant location flags ! logical :: lforce_locate_SNI=.false. logical :: uniform_zdist_SNI = .false. logical :: lOB_cluster = .false. ! SN clustering real :: p_OB=0.7 ! Probability that an SN is in a cluster real :: SN_clustering_radius=impossible real :: SN_clustering_radius_cgs=9.258e20 ! cm (300 pc) ! ! Adjust SNR%radius inversely with density ! logical :: lSN_scale_rad=.false. real :: N_mass=60.0 ! ! Requested SNe location (used for test SN) ! real :: center_SN_x = impossible real :: center_SN_y = impossible real :: center_SN_z = impossible ! ! Volume element ! real :: dv ! ! Cooling time diagnostic ! integer :: idiag_taucmin=0 integer :: idiag_Hmax=0 integer :: idiag_Lamm=0 integer :: idiag_nrhom=0 integer :: idiag_rhoLm=0 integer :: idiag_Gamm=0 ! ! Heating function, cooling function and mass movement ! method selection. ! character (len=labellen) :: cooling_select = 'RBN' character (len=labellen) :: heating_select = 'wolfire' character (len=labellen) :: thermal_profile = 'gaussian3' character (len=labellen) :: velocity_profile= 'lineartanh' character (len=labellen) :: mass_profile = 'gaussian3' character (len=labellen) :: mass_movement = 'off' character (len=labellen) :: cavity_profile = 'gaussian3' ! ! Variables required for returning mass to disk given no inflow ! boundary condition used in addmassflux ! real :: addflux_dim1, addrate=1.0 real :: boldmass=0.0 logical :: ladd_massflux = .false. ! ! start parameters ! namelist /interstellar_init_pars/ & initinterstellar, initial_SNI, h_SNI, h_SNII, lSNII, lSNI, & lSN_scale_rad, ampl_SN, kampl_SN, mass_SN, velocity_SN, width_SN, & mass_width_ratio, energy_width_ratio, velocity_width_ratio, & t_next_SNI, t_next_SNII, center_SN_x, center_SN_y, center_SN_z, & lSN_velocity, lSN_eth, lSN_ecr, lSN_fcr, lSN_mass, mass_movement, & uu_sedov_max, frac_ecr, frac_eth, thermal_profile, velocity_profile, & mass_profile, uniform_zdist_SNI, inner_shell_proportion, & outer_shell_proportion, SNR_damping, cooling_select, heating_select, & heating_rate, rho0ts, T0hs, TT_SN_max, rho_SN_min, N_mass, & lSNII_gaussian, rho_SN_max, lthermal_hse, lheatz_min ! ! run parameters ! namelist /interstellar_run_pars/ & ampl_SN, kampl_SN, mass_SN, velocity_SN, t_next_SNI, t_next_SNII, & mass_width_ratio, energy_width_ratio, velocity_width_ratio, & lSN_velocity, lSN_eth, lSN_ecr, lSN_fcr, lSN_mass, width_SN, lSNI, & lSNII, uniform_zdist_SNI, mass_movement, SNI_area_rate, SNII_area_rate, & inner_shell_proportion, outer_shell_proportion, uu_sedov_max, & frac_ecr, frac_eth, thermal_profile,velocity_profile, mass_profile, & h_SNI, h_SNII, TT_SN_min, SNR_damping, uu_sedov_max, lSN_scale_rad, & mass_SN_progenitor, cloud_tau, cdt_tauc, cloud_rho, cloud_TT, & laverage_SN_heating, coolingfunction_scalefactor, lforce_locate_SNI,& lsmooth_coolingfunc, heatingfunction_scalefactor, t_settle, & center_SN_x, center_SN_y, center_SN_z, rho_SN_min, TT_SN_max, & lheating_UV, cooling_select, heating_select, heating_rate, & lcooltime_smooth, lcooltime_despike, cooltime_despike_factor, & heatcool_shock_cutoff, heatcool_shock_cutoff_rate, ladd_massflux, & lTT_cutoff, TT_cutoff, N_mass, addrate, T0hs, rho0ts, & lSNII_gaussian, rho_SN_max, lSN_mass_rate, lthermal_hse, lheatz_min, & p_OB, SN_clustering_radius, lOB_cluster ! contains ! !*********************************************************************** subroutine register_interstellar() ! ! 19-nov-02/tony: coded ! use FArrayManager ! call farray_register_auxiliary('netcool',inetheat,communicated=.true.) call farray_register_auxiliary('cooling',icooling) ! ! identify version number ! if (lroot) call svn_id( & "$Id$") ! ! Invalidate all SNRs ! nSNR=0 SNRs(:)%state=SNstate_invalid ! ! Writing files for use with IDL ! if (naux+naux_com < maux+maux_com) aux_var(aux_count)=',netcool $' if (naux+naux_com == maux+maux_com) aux_var(aux_count)=',netcool' aux_count=aux_count+1 if (lroot) write(15,*) 'netcool = fltarr(mx,my,mz)*one' ! endsubroutine register_interstellar !*********************************************************************** subroutine initialize_interstellar(f) ! ! Perform any post-parameter-read initialization eg. set derived ! parameters ! ! 24-nov-02/tony: coded ! ! read parameters from seed.dat and interstellar.dat ! use General, only: random_seed_wrapper use Mpicomm, only: stop_it use EquationOfState, only: getmu ! real, dimension (mx,my,mz,mfarray) :: f ! real :: mu ! f(:,:,:,icooling)=0.0 f(:,:,:,inetheat)=0.0 ! if (lroot) print*,'initialize_interstellar: t_next_SNI',t_next_SNI ! if (lroot.and.uniform_zdist_SNI) then print*,'initialize_interstellar: using UNIFORM z-distribution of SNI' endif ! dv=1. if (nxgrid/=1) dv=dv*dx if (nygrid/=1) dv=dv*dy if (nzgrid/=1) dv=dv*dz ! call getmu(f,mu) if (unit_system=='cgs') then ! ! this Lambda as such enters as n^2*Lambda(T) on the rhs of the ! energy equation per unit volume ! unit_Lambda = unit_velocity**2 / unit_density / unit_time elseif (unit_system=='SI') then call stop_it('initialize_interstellar: SI unit conversions not implemented') endif if (lroot) print*,'initialize_interstellar: unit_Lambda',unit_Lambda ! ! Mara: Initialize cooling parameters according to selection ! Default selection 'RBN' Rosen & Bregman (1993) ! Alternative selection 'SS' Sanchez-Salcedo et al. (2002) ! Turn off cooling: cooling_select='off' ! cooling_select in interstellar_init_pars added ! if (cooling_select == 'RBN') then if (lroot) print*,'initialize_interstellar: default RBN cooling fct' coolT_cgs = (/ 100.0, & 2000.0, & 8000.0, & 1.0E5, & 4.0E7, & 1.0E9, & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0) /) coolH_cgs = (/ 2.2380D-32, & 1.0012D-30, & 4.6240D-36, & 1.7800D-18, & 3.2217D-27, & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0) /) / ( m_p_cgs )**2 coolB = (/ 2.0, & 1.5, & 2.867, & -0.65, & 0.5, & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.) /) ncool = 5 ! ! 04-jan-10/fred: ! Reset above to original RBN parameters. Altered coolB(5) and coolT(5,6) in ! RBNr to ensure continuity and increase cooling at high temperatures later in ! diffuse remnant cores. ! RBNr: new lower terms for smooth cooling below 300K ! and extended range to 1E13 in SSrr to deal with temperature spiking ! else if (cooling_select == 'RBNr') then if (lroot) print*,'initialize_interstellar: RBN cooling fct (revised)' coolT_cgs = (/ 10.0, & 2000.0, & 8000.0, & 1E5, & 1E6, & 1E17, & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0) /) coolH_cgs = (/ 2.2380D-32, & 1.0012D-30, & 4.6240D-36, & 1.7783524D-18,& 2.238814D-25, & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0) /) / ( m_p_cgs )**2 coolB = (/ 2.0, & 1.5, & 2.867, & -0.65, & 0.5, & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.) /) ncool=5 else if (cooling_select == 'SS') then ! ! These are the SS et al (2002) coefficients multiplied by m_proton**2 ! coolT_cgs = (/ 10.0, & 141.0, & 313.0, & 6102.0, & 1.0E5, & 1.0E17, & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0) /) coolH_cgs = (/ 3.42D16, & 9.10D18, & 1.11D20, & 2.00D8, & 7.962D29, & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0) /) coolB = (/ 2.12, & 1.0, & 0.56, & 3.67, & -0.65, & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.) /) ncool=5 ! else if (cooling_select == 'SSr') then if (lroot) print*,'initialize_interstellar: revised SS cooling fct' coolT_cgs = (/ 10.0, & 141.0, & 313.0, & 6102.0, & 1E5, & 1E9, & 1E17, & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0E0) /) coolH_cgs = (/ 3.70D16, & 9.46D18, & 1.185D20, & 2.00D8, & 7.96D29, & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0) /) coolB = (/ 2.12, & 1.0, & 0.56, & 3.67, & -0.65, & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.) /) ncool=5 ! ! Revised to make continuous ! else if (cooling_select == 'SSrr') then if (lroot) print*,'initialize_interstellar: revised SS cooling fct' coolT_cgs = (/ 10.0, & 141.0, & 313.0, & 6102.0, & 1E5, & 4E7, & 1E17, & tiny(0E0), & tiny(0E0), & tiny(0E0), & tiny(0.0) /) coolH_cgs = (/ 3.703109927416290D16, & 9.455658188464892D18, & 1.185035244783337D20, & 1.9994576479D8, & 7.96D29, & 1.440602814622207D21, & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0), & tiny(0D0) /) coolB = (/ 2.12, & 1.0, & 0.56, & 3.67, & -0.65, & 0.5, & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.), & tiny(0.) /) ncool=6 ! ! 26-Jan-10/fred ! Combines Sanchez-Salcedo (2002) with Slyz et al (2005) above 1E5K ! as Gressel simulation (2008) with constants revised for continuity ! else if (cooling_select == 'WSW') then if (lroot) print*,'initialize_interstellar: WSW cooling fct' coolT_cgs = (/ 10.0, & 141.0, & 313.0, & 6102.0, & 1E5, & 2.88E5, & 4.73E5, & 2.11E6, & 3.98E6, & 2.0E7, & 1.0E17 /) coolH_cgs = (/ 3.703109927416290D16, & 9.455658188464892D18, & 1.185035244783337D20, & 1.102120336D10, & 1.236602671D27, & 2.390722374D42, & 4.003272698D26, & 1.527286104D44, & 1.608087849D22, & 9.228575532D20, & tiny(0D0) /) coolB = (/ 2.12, & 1.0, & 0.56, & 3.21, & -0.20, & -3.0, & -0.22, & -3.00, & 0.33, & 0.50, & tiny(0.) /) ncool=10 ! ! As above but with higher minimum temperature 90K instead of 10K ! else if (cooling_select == 'WSWr') then if (lroot) print*,'initialize_interstellar: WSWr cooling fct' coolT_cgs = (/ 90.0, & 141.0, & 313.0, & 6102.0, & 1E5, & 2.88E5, & 4.73E5, & 2.11E6, & 3.98E6, & 2.0E7, & 1.0E17 /) coolH_cgs = (/ 3.703109927416290D16, & 9.455658188464892D18, & 1.185035244783337D20, & 1.102120336D10, & 1.236602671D27, & 2.390722374D42, & 4.003272698D26, & 1.527286104D44, & 1.608087849D22, & 9.228575532D20, & tiny(0D0) /) coolB = (/ 2.12, & 1.0, & 0.56, & 3.21, & -0.20, & -3.0, & -0.22, & -3.00, & 0.33, & 0.5, & tiny(0.) /) ncool=10 else if (cooling_select == 'off') then if (lroot) print*,'initialize_interstellar: no cooling applied' coolT_cgs=tiny(0.0) coolH_cgs=tiny(0.0) coolB=tiny(0.) endif ! ! BEGIN TEMPORARY if (any(coolH_cgs(1:ncool) == 0) & .or. any(coolT_cgs(1:ncool+1) == 0)) then call fatal_error('initialize_interstellar', & 'Calculating lncoolH and lncoolT: One of the cooling coefficient is zero') endif ! END TEMPORARY lncoolH(1:ncool) = real(log(coolH_cgs(1:ncool)) - log(unit_Lambda) & + log(unit_temperature**coolB(1:ncool)) & + log(coolingfunction_scalefactor)) lncoolT(1:ncool+1) = real(log(coolT_cgs(1:ncool+1) / unit_temperature)) ! heating_rate_code=heating_rate*real(unit_length/unit_velocity**3) ! if (unit_system=='cgs') then if (TT_SN_max==impossible) TT_SN_max=TT_SN_max_cgs / unit_temperature if (rho_SN_min==impossible) rho_SN_min=rho_SN_min_cgs / unit_density if (rho_SN_max==impossible) rho_SN_max=rho_SN_max_cgs / unit_density TT_SN_min=TT_SN_min_cgs / unit_temperature TT_cutoff=TT_cutoff_cgs / unit_temperature TT_cutoff1=1.0/TT_cutoff if (SNI_area_rate==impossible) & SNI_area_rate=SNI_area_rate_cgs * unit_length**2 * unit_time if (SNII_area_rate==impossible) & SNII_area_rate=7.5*SNI_area_rate_cgs * unit_length**2 * unit_time if (h_SNI==impossible) h_SNI=h_SNI_cgs / unit_length SNII_mass_rate=SNII_mass_rate_cgs*unit_time SNI_mass_rate=SNI_mass_rate_cgs*unit_time h_SNII=h_SNII_cgs / unit_length solar_mass=solar_mass_cgs / unit_mass if (lroot) & print*,'initialize_interstellar: solar_mass (code) =', solar_mass if (cloud_rho==impossible) cloud_rho=cloud_rho_cgs / unit_density if (cloud_TT==impossible) cloud_TT=cloud_TT_cgs / unit_temperature r_SNI =r_SNI_yrkpc2 * (unit_time/yr_cgs) * (unit_length/kpc_cgs)**2 r_SNII=r_SNII_yrkpc2 * (unit_time/yr_cgs) * (unit_length/kpc_cgs)**2 T0UV=T0UV_cgs / unit_temperature cUV=cUV_cgs * unit_temperature if (GammaUV==impossible) & GammaUV=GammaUV_cgs * real(unit_length/unit_velocity**3) if (ampl_SN==impossible) ampl_SN=ampl_SN_cgs / unit_energy if (kampl_SN==impossible) then if (.not.lSN_velocity) then kampl_SN=0.0 else ampl_SN=0.5*ampl_SN kampl_SN=ampl_SN endif endif if (rho_min == impossible) rho_min=rho_min_cgs/unit_temperature if (T0hs == impossible) T0hs=T0hs_cgs/unit_temperature if (rho0ts == impossible) rho0ts=rho0ts_cgs/unit_density if (lroot) & print*,'initialize_interstellar: ampl_SN, kampl_SN = ', & ampl_SN, kampl_SN if (cloud_tau==impossible) cloud_tau=cloud_tau_cgs / unit_time if (mass_SN==impossible) mass_SN=mass_SN_cgs / unit_mass if (mass_SN_progenitor==impossible) & mass_SN_progenitor=mass_SN_progenitor_cgs / unit_mass if (width_SN==impossible) width_SN= & max(width_SN_cgs / real(unit_length),dxmax*2.5) if (SNR_damping_time==impossible) & SNR_damping_time=SNR_damping_time_cgs / unit_time if (SNR_damping_rate==impossible) & SNR_damping_rate=SNR_damping_rate_cgs / unit_time if (SN_clustering_radius==impossible) & SN_clustering_radius=SN_clustering_radius_cgs / unit_length else call stop_it('initialize_interstellar: SI unit conversions not implemented') endif ! ! Inverse volume share of mass lost through the boundary substitute for ! galactic fountain given no inflow boundary vertical condition ! if (ladd_massflux) addflux_dim1=1./(Lxyz(1)*Lxyz(2)*Lxyz(3)) ! if (heating_select == 'thermal-hs') then call thermal_hs(f,zrho) call heat_interstellar(f,heat_z,zrho) endif ! ! Cooling cutoff in shocks ! if (heatcool_shock_cutoff_rate/=0.) then lheatcool_shock_cutoff=.true. heatcool_shock_cutoff_rate1=1.0/heatcool_shock_cutoff_rate else lheatcool_shock_cutoff=.false. endif ! ! SNRdamping factor ! if (SNR_damping/=0.) lSNR_damping=.true. ! ! Slopeyness used for tanh rounding profiles etc. ! sigma_SN=dxmax*3 sigma_SN1=1./sigma_SN ! preSN(:,:)=0 ! t_interval_SNI = 1./(SNI_area_rate * Lxyz(1) * Lxyz(2)) t_interval_SNII = 1./(SNII_area_rate * Lxyz(1) * Lxyz(2)) average_SNI_heating = & r_SNI *ampl_SN/(sqrt(pi)*h_SNI )*heatingfunction_scalefactor average_SNII_heating= & r_SNII*ampl_SN/(sqrt(pi)*h_SNII)*heatingfunction_scalefactor if (lroot) print*,'initialize_interstellar: t_interval_SNI =', & t_interval_SNI,Lxyz(1),Lxyz(2),SNI_area_rate ! if (lroot .and. (ip<14)) then print*,'initialize_interstellar: nseed,seed',nseed,seed(1:nseed) print*,'initialize_interstellar: finished' endif ! if (lroot .and. lstart) then open(1,file=trim(datadir)//'/sn_series.dat',position='append') write(1,'("#",4A)') & '---it----------t--------itype-iproc----l-----m----n---', & '-----x------------y------------z-------', & '----rho-----------TT-----------EE---------t_sedov----', & '--radius------site_mass------maxTT----t_interval---' close(1) endif ! ! Write unit_Lambda to pc_constants file ! if (lroot) then print*,"initialize_interstellar: t_next_SNI, t_next_SNII=", & t_next_SNI, t_next_SNII open (1,file=trim(datadir)//'/pc_constants.pro',position="append") write (1,'(a,1pd26.16)') 'unit_Lambda=',unit_Lambda close (1) endif ! endsubroutine initialize_interstellar !***************************************************************************** subroutine input_persistent_interstellar(id,done) ! ! Read in the stored time of the next SNI ! ! 13-Dec-2011/Bourdin.KIS: reworked ! use IO, only: read_persist, lun_input, lcollective_IO ! integer :: id logical :: done ! integer :: i ! if (lcollective_IO) call fatal_error ('input_persistent_interstellar', & "The interstellar persistent variables can't be read collectively!") ! select case (id) ! for backwards-compatibility (deprecated): case (id_record_ISM_T_NEXT_OLD) read (lun_input) t_next_SNI, t_next_SNII done = .true. case (id_record_ISM_POS_NEXT_OLD) read (lun_input) x_next_SNII, y_next_SNII done = .true. case (id_record_ISM_TOGGLE_OLD) read (lun_input) lSNI, lSNII done = .true. case (id_record_ISM_SNRS_OLD) ! Forget any existing SNRs. SNRs(:)%state = SNstate_invalid read (lun_input) nSNR do i = 1, nSNR read (lun_input) SNRs(i) SNR_index(i) = i enddo done = .true. ! currently active tags: case (id_record_ISM_T_NEXT_SNI) if (read_persist ('ISM_T_NEXT_SNI', t_next_SNI)) return done = .true. case (id_record_ISM_T_NEXT_SNII) if (read_persist ('ISM_T_NEXT_SNII', t_next_SNII)) return done = .true. case (id_record_ISM_X_CLUSTER) if (read_persist ('ISM_X_CLUSTER', x_next_SNII)) return done = .true. case (id_record_ISM_Y_CLUSTER) if (read_persist ('ISM_Y_CLUSTER', y_next_SNII)) return done = .true. case (id_record_ISM_TOGGLE_SNI) if (read_persist ('ISM_TOGGLE_SNI', lSNI)) return done = .true. case (id_record_ISM_TOGGLE_SNII) if (read_persist ('ISM_TOGGLE_SNII', lSNII)) return done = .true. case (id_record_ISM_BOLD_MASS) if (read_persist ('ISM_BOLD_MASS', boldmass)) return done = .true. case (id_record_ISM_SNRS) ! Forget any existing SNRs. SNRs(:)%state = SNstate_invalid read (lun_input) nSNR do i = 1, nSNR read (lun_input) SNRs(i) SNR_index(i) = i enddo done = .true. endselect ! if (lroot) & print *, 'input_persistent_interstellar: ', t_next_SNI, t_next_SNII ! endsubroutine input_persistent_interstellar !***************************************************************************** logical function output_persistent_interstellar() ! ! Writes out the time of the next SNI ! ! 13-Dec-2011/Bourdin.KIS: reworked ! use IO, only: write_persist, write_persist_id, lun_output, lcollective_IO ! integer :: i ! if (lcollective_IO) call fatal_error ('output_persistent_interstellar', & "The interstellar persistent variables can't be written collectively!") ! output_persistent_interstellar = .true. ! if (write_persist ('ISM_T_NEXT_SNI', id_record_ISM_T_NEXT_SNI, t_next_SNI)) return if (write_persist ('ISM_T_NEXT_SNII', id_record_ISM_T_NEXT_SNII, t_next_SNII)) return if (write_persist ('ISM_X_CLUSTER', id_record_ISM_X_CLUSTER, x_next_SNII)) return if (write_persist ('ISM_Y_CLUSTER', id_record_ISM_Y_CLUSTER, y_next_SNII)) return if (write_persist ('ISM_TOGGLE_SNI', id_record_ISM_TOGGLE_SNI, lSNI)) return if (write_persist ('ISM_TOGGLE_SNII', id_record_ISM_TOGGLE_SNII, lSNII)) return if (write_persist ('ISM_BOLD_MASS', id_record_ISM_BOLD_MASS, boldmass)) return ! ! For self-defined data types, things are not so easy to implement. if (write_persist_id ('ISM_SNRS', id_record_ISM_SNRS)) return write (lun_output) nSNR do i = 1, nSNR write (lun_output) SNRs(SNR_index(i)) enddo ! output_persistent_interstellar = .false. ! endfunction output_persistent_interstellar !***************************************************************************** subroutine rprint_interstellar(lreset,lwrite) ! ! Reads and registers print parameters relevant to interstellar ! ! 01-jun-02/axel: adapted from magnetic fields ! use Diagnostics, only: parse_name ! integer :: iname logical :: lreset logical, intent(in), optional :: lwrite ! ! Reset everything in case of reset ! (This needs to be consistent with what is defined above!) ! if (lreset) then SNR_damping=0. idiag_taucmin=0 idiag_Hmax=0 idiag_Lamm=0 idiag_nrhom=0 idiag_rhoLm=0 idiag_Gamm=0 ! ! TT_SN_max=impossible ! rho_SN_min=impossible ! SNI_area_rate=impossible ! h_SNI=impossible ! GammaUV=impossible ! width_SN=impossible ! ampl_SN=impossible ! mass_SN=impossible ! mass_SN_progenitor=impossible ! cloud_tau=impossible endif ! lpenc_requested(i_ee)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_TT1)=.true. ! ! iname runs through all possible names that may be listed in print.in ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'taucmin',idiag_taucmin) call parse_name(iname,cname(iname),cform(iname),'Hmax',idiag_Hmax) call parse_name(iname,cname(iname),cform(iname),'Lamm',idiag_Lamm) call parse_name(iname,cname(iname),cform(iname),'nrhom',idiag_nrhom) call parse_name(iname,cname(iname),cform(iname),'rhoLm',idiag_rhoLm) call parse_name(iname,cname(iname),cform(iname),'Gamm',idiag_Gamm) enddo ! ! check for those quantities for which we want video slices ! if (lwrite_slices) then where(cnamev=='ism_cool'.or.cnamev=='ism_cool2') cformv='DEFINED' endif ! endsubroutine rprint_interstellar !***************************************************************************** subroutine get_slices_interstellar(f,slices) ! ! Write slices for animation of Interstellar variables. ! ! 26-jul-06/tony: coded ! use Slices_methods, only: assign_slices_scal real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! ! Loop over slices ! select case (trim(slices%name)) ! ! Shock profile ! case ('ism_cool') call assign_slices_scal(slices,f,icooling) case ('ism_cool2') call assign_slices_scal(slices,f,inetheat) ! endselect ! endsubroutine get_slices_interstellar !*********************************************************************** subroutine read_interstellar_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=interstellar_init_pars, IOSTAT=iostat) ! endsubroutine read_interstellar_init_pars !*********************************************************************** subroutine write_interstellar_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=interstellar_init_pars) ! endsubroutine write_interstellar_init_pars !*********************************************************************** subroutine read_interstellar_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=interstellar_run_pars, IOSTAT=iostat) ! endsubroutine read_interstellar_run_pars !*********************************************************************** subroutine write_interstellar_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=interstellar_run_pars) ! endsubroutine write_interstellar_run_pars !*********************************************************************** subroutine init_interstellar(f) ! ! Initialise some explosions etc. ! 24-nov-2002/tony: coded ! use General, only: itoa ! real, dimension (mx,my,mz,mfarray) :: f ! logical :: lnothing=.true. character (len=intlen) :: iinit_str integer :: i,j,iSNR ! intent(inout) :: f ! do j=1,ninit ! if (initinterstellar(j)/='nothing') then ! lnothing=.false. iinit_str=itoa(j) ! ! Select different initial conditions ! select case (initinterstellar(j)) ! case ('single') iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN call position_SN_testposition(f,SNRs(iSNR)) call explode_SN(f,SNRs(iSNR)) lSNI=.false. lSNII=.false. case ('sedov') iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN center_SN_x=0. center_SN_y=0. center_SN_z=0. call position_SN_testposition(f,SNRs(iSNR)) call explode_SN(f,SNRs(iSNR)) lSNI=.false. lSNII=.false. case ('courant-friedricks') iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN center_SN_x=0. center_SN_y=0. center_SN_z=-0.015 call position_SN_testposition(f,SNRs(iSNR)) call explode_SN(f,SNRs(iSNR)) iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=t SNRs(iSNR)%radius=width_SN center_SN_x=0. center_SN_y=0. center_SN_z=0.015 call position_SN_testposition(f,SNRs(iSNR)) call explode_SN(f,SNRs(iSNR)) lSNI=.false. lSNII=.false. case ('kompaneets') iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=0 SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN center_SN_x=0. center_SN_y=0. center_SN_z=0. call position_SN_testposition(f,SNRs(iSNR)) call explode_SN(f,SNRs(iSNR)) lSNI=.false. lSNII=.false. case ('multiple') do i = 1,initial_SNI iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0. SNRs(iSNR)%t=0. SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN if (uniform_zdist_SNI) then call position_SN_uniformz(f,SNRs(i)) else call position_SN_gaussianz(f,h_SNII,SNRs(i)) endif call explode_SN(f,SNRs(i)) enddo ! case default ! ! Catch unknown values ! write(unit=errormsg,fmt=*) 'No such value for initinterstellar(' & //trim(iinit_str)//'): ',trim(initinterstellar(j)) call fatal_error('init_interstellar',errormsg) ! endselect ! if (lroot) print*,'init_interstellar: initinterstellar(' & //trim(iinit_str)//') = ',trim(initinterstellar(j)) endif ! enddo ! if (lnothing.and.lroot) print*,'init_interstellar: nothing' ! call tidy_SNRs ! endsubroutine init_interstellar !***************************************************************************** subroutine pencil_criteria_interstellar() ! ! All pencils that the Interstellar module depends on are specified here. ! ! 26-mar-05/tony: coded ! 11-mar-06/axel: added idiag_nrhom ! lpenc_requested(i_ee)=.true. lpenc_requested(i_lnrho)=.true. lpenc_requested(i_lnTT)=.true. lpenc_requested(i_TT1)=.true. lpenc_requested(i_rho1)=.true. ! if (lheatcool_shock_cutoff) lpenc_requested(i_gshock)=.true. ! ! Diagnostic pencils ! ! AB: ! if (idiag_nrhom/=0) lpenc_diagnos(i_rho)=.true. ! endsubroutine pencil_criteria_interstellar !***************************************************************************** subroutine interstellar_before_boundary(f) ! ! This routine calculates and applies the optically thin cooling function ! together with UV heating. ! ! 01-aug-06/tony: coded ! use Diagnostics, only: max_mn_name, sum_mn_name use EquationOfState, only: gamma, gamma1, eoscalc ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! real, dimension (nx) :: heat,cool,lnTT, lnrho real, dimension (nx) :: damp_profile real :: minqty integer :: i,iSNR ! if (.not.(lcooltime_smooth.or.lcooltime_despike)) return ! ! Identifier ! if (headtt) print*,'interstellar_before_boundary: ENTER' ! ! Precalculate radiative cooling function ! do n=n1,n2 do m=m1,m2 if (ldensity_nolog) then lnrho = log(f(l1:l2,m,n,irho)) else lnrho = f(l1:l2,m,n,ilnrho) endif call eoscalc(f,nx,lnTT=lnTT) ! call calc_cool_func(cool,lnTT,lnrho) call calc_heat(heat,lnTT) ! if (nSNR==0) then ! ! 05-sep-10/fred ! NB The applied net heating/cooling: (heat-cool)/temp is stored in ! inetheat. The radiative cooling rho*Lambda is stored in icooling. ! Both are diagnostic. ! if (ltemperature) then f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool)*gamma elseif (pretend_lnTT) then f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool)*gamma else f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool) endif else damp_profile=1. do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state==SNstate_damping) then call proximity_SN(SNRs(iSNR)) minqty=0.5*(1.+tanh((t-SNRs(iSNR)%t_damping)/SNR_damping_rate)) damp_profile = damp_profile * ( (1.-minqty) * 0.5 * (1.+ & tanh((sqrt(dr2_SN)-(SNRs(iSNR)%radius*2.))*sigma_SN1-2.)) & + minqty ) endif enddo if (ltemperature) then f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool)*gamma*damp_profile elseif (pretend_lnTT) then f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool)*gamma*damp_profile else f(l1:l2,m,n,inetheat)=exp(-lnTT)*(heat-cool)*damp_profile endif endif enddo enddo ! ! endsubroutine interstellar_before_boundary !***************************************************************************** subroutine thermal_hs(f,zrho) ! ! This routine calculates a vertical profile for density for an appropriate ! isothermal entropy designed to balance the vertical 'Ferriere' gravity. ! T0hs and rho0ts are chosen to ensure uv-heating approx 0.0147 at z=0. ! Initial thermal & hydrostatice equilibrium is achieved by ensuring ! Lambda*rho(z)=Gamma(z). ! ! Requires gravz_profile='Ferriere' in gravity_simple.f90, ! init_lnrho & initss='thermal-hs' in density & entropy.f90. ! Constants g_A..D from gravz_profile. ! ! 22-mar-10/fred: coded ! 12-aug-10/fred: updated ! use SharedVariables, only: put_shared_variable use EquationOfState , only: getmu ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension(mz), intent(out) :: zrho ! real :: logrho real :: muhs real :: g_A, g_C real, parameter :: g_A_cgs=4.4E-9, g_C_cgs=1.7E-9 real :: g_B ,g_D real, parameter :: g_B_cgs=6.172E20 , g_D_cgs=3.086E21 integer :: ierr ! ! Identifier ! if (lroot.and.headtt.and.ip<14) print*,'thermal_hs: ENTER' ! ! Set up physical units. ! if (unit_system=='cgs') then g_A = g_A_cgs/unit_velocity*unit_time g_B = g_B_cgs/unit_length g_C = g_C_cgs/unit_velocity*unit_time g_D = g_D_cgs/unit_length else if (unit_system=='SI') then call fatal_error('initialize_entopy', & 'SI unit conversions not inplemented') endif ! ! Uses gravity profile from K. Ferriere, ApJ 497, 759, 1998, eq (34) ! at solar radius. ! call getmu(f,muhs) ! if (lroot) print*, 'thermal-hs: '// & 'hydrostatic thermal equilibrium density and entropy profiles' ! do n=1,mz if (lthermal_hse) then logrho = log(rho0ts)+(g_A*g_B*m_u*muhs/k_B/T0hs)*(log(T0hs)- & log(T0hs/(g_A*g_B)* & (g_A*sqrt(g_B**2+(z(n))**2)+0.5*g_C*(z(n))**2/g_D))) else logrho = log(rho0ts)-0.015*(- & g_A*g_B+ & g_A*sqrt(g_B**2+(z(n))**2)+0.5*g_C*(z(n))**2/g_D) endif if (logrho < -40.0) logrho=-40.0 zrho(n)=exp(logrho) enddo ! ! Share zrho and T0hs for use with entropy to initialize density and ! temperature in thermal_hs_equilibrium_ism in entropy ! call put_shared_variable('zrho', zrho, ierr) if (ierr/=0) call fatal_error('thermal_hs', & 'there was a problem when putting zrho') call put_shared_variable('T0hs', T0hs, ierr) if (ierr/=0) call fatal_error('thermal_hs', & 'there was a problem when putting T0hs') ! endsubroutine thermal_hs !***************************************************************************** subroutine heat_interstellar(f,zheat,zrho) ! ! This routine calculates a vertical profile for uv-heating designed to ! satisfy an initial condition with heating and cooling balanced for an ! isothermal hydrostatic equilibrium. ! Requires: gravz_profile='Ferriere' in gravity_simple.f90 ! initlnrho='thermal-hs' in density.f90 ! initss='thermal-hs' in entropy.f90 ! heating_select='thermal-hs' in interstellar.f90 ! Using here a similar method to O. Gressel 2008 (PhD) lthermal_hse=T ! or similar to Joung & Mac Low Apj 653 Dec 2006 without hse ! ! 22-mar-10/fred: ! adapted from galactic-hs,ferriere-hs ! 12-aug-10/fred: ! included zrho & T0hs from thermal_hs ! use EquationOfState , only: getmu ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f real, dimension(mz), intent(in) :: zrho real, dimension(mz), intent(out) :: zheat ! real :: g_A, g_C real, parameter :: g_A_cgs=4.4E-9, g_C_cgs=1.7E-9 real :: g_B ,g_D, H_z real, parameter :: g_B_cgs=6.172E20 , g_D_cgs=3.086E21, & H_z_cgs=9.258E20 real, dimension(mz) :: lambda=0.0, lnTT, TT integer :: j ! ! Identifier ! if (lroot.and.headtt.and.ip<14) print*,'heat_interstellar: ENTER' ! if (lroot) print*, & 'heat_interstellar: calculating z-dependent uv-heating'// & 'function for initial hydrostatic and thermal equilibrium' ! ! Set up physical units. ! if (unit_system=='cgs') then g_A = g_A_cgs/unit_velocity*unit_time g_C = g_C_cgs/unit_velocity*unit_time g_D = g_D_cgs/unit_length g_B = g_B_cgs/unit_length H_z = H_z_cgs/unit_length else if (unit_system=='SI') then call fatal_error('initialize_entopy', & 'SI unit conversions not inplemented') endif ! do n=1,mz TT(n)=T0hs/(g_A*g_B)* & (g_A*sqrt(g_B**2+(z(n))**2)+0.5*g_C*(z(n))**2/g_D) lnTT(n)=log(TT(n)) zheat(n)=GammaUV*exp(-abs(z(n))/H_z) enddo lam_loop: do j=1,ncool if (lncoolT(j) >= lncoolT(j+1)) exit lam_loop where (lncoolT(j)<=lnTT.and.lnTTTT_cutoff1) heatcool=TT_cutoff1*(heat-cool)*gamma endif elseif (pretend_lnTT) then heatcool=p%TT1*(heat-cool)*gamma if (lTT_cutoff) then where (p%TT1>TT_cutoff1) heatcool=TT_cutoff1*(heat-cool)*gamma endif else heatcool=p%TT1*(heat-cool) if (lTT_cutoff) then where (p%TT1>TT_cutoff1) heatcool=TT_cutoff1*(heat-cool) endif endif ! ! Prevent unresolved heating/cooling in early SNR core. ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state==SNstate_damping) then call proximity_SN(SNRs(iSNR)) minqty=0.5*(1.+tanh((t-SNRs(iSNR)%t_damping)/SNR_damping_rate)) damp_profile = & minqty + (1.-minqty) * 0.5 * & (1.+tanh((sqrt(dr2_SN)-(SNRs(iSNR)%radius*2.))*sigma_SN1-2.)) heatcool=heatcool*damp_profile cool=cool*damp_profile endif enddo ! ! Prevent unresolved heating/cooling in shocks. This is recommended as ! early cooling in the shock prematurely inhibits the strength of the ! shock wave and also drives down the timestep. Fred ! if (lheatcool_shock_cutoff) then call dot2(p%gshock,gsh2) ! damp_profile=exp(-(gsh2*heatcool_shock_cutoff_rate1)) ! cool=cool*damp_profile heat=heat*damp_profile heatcool=heatcool*damp_profile endif ! ! Save result in diagnostic aux variable ! cool=rho*Lambda, heatcool=(Gamma-rho*Lambda)/TT ! f(l1:l2,m,n,icooling)=cool f(l1:l2,m,n,inetheat)=heatcool ! ! Average SN heating (due to SNI and SNII) ! The amplitudes of both types is assumed the same (=ampl_SN) ! if (laverage_SN_heating) then heat=heat+average_SNI_heating *exp(-(z(n)/h_SNI )**2) heat=heat+average_SNII_heating*exp(-(z(n)/h_SNII)**2) endif ! ! Prepare diagnostic output ! Since these variables are divided by Temp when applied it is useful to ! monitor the actual applied values for diagnostics so TT1 included. ! if (ldiagnos) then if (idiag_Hmax/=0) then netheat=heatcool where (heatcool<0.0) netheat=0.0 call max_mn_name(netheat/p%ee,idiag_Hmax) endif if (idiag_taucmin/=0) then netcool=-heatcool where (heatcool>=0.0) netcool=1.0 call max_mn_name(netcool/p%ee,idiag_taucmin,lreciprocal=.true.) endif if (idiag_Lamm/=0) & call sum_mn_name(p%rho1*cool*p%TT1,idiag_Lamm) if (idiag_nrhom/=0) & call sum_mn_name(cool/p%ee,idiag_nrhom) ! -- call sum_mn_name(cool*p%rho/p%ee,idiag_nrhom) ! AB: the factor rho is already included in cool, so cool=rho*Lambda if (idiag_rhoLm/=0) & call sum_mn_name(p%TT1*cool,idiag_rhoLm) if (idiag_Gamm/=0) & call sum_mn_name(p%TT1*heat,idiag_Gamm) endif ! ! Limit timestep by the cooling time (having subtracted any heating) ! dt1_max=max(dt1_max,cdt_tauc*(cool)/ee,cdt_tauc*(heat)/ee) ! if (lfirst.and.ldt) then dt1_max=max(dt1_max,(-heatcool)/(p%ee*cdt_tauc)) where (heatcool>0.0) Hmax=Hmax+heatcool dt1_max=max(dt1_max,Hmax/(p%ee*cdt_tauc)) endif ! ! Apply heating/cooling to temperature/entropy variable ! if (ltemperature) then df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT)+heatcool else df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss)+heatcool endif ! endsubroutine calc_heat_cool_interstellar !***************************************************************************** subroutine calc_cool_func(cool,lnTT,lnrho) ! ! This routine calculates the temperature dependent radiative cooling. ! Applies Rosen et al., ApJ, 413, 137, 1993 ('RBN') OR ! Sanchez-Salcedo et al. ApJ, 577, 768, 2002 ('SS') OR ! Slyz et al MNRAS, 356 2005 ('WSW') fit of Wolfire with Sarazin & White ! ApJ, 443:152-168, 1985 and ApJ, 320:32-48, 1987 ! ! ! Cooling is Lambda*rho^2, with (eq 7) & Lambda=coolH(i)*TT**coolB(i), ! for coolT(i) <= TT < coolT(i+1). Nb: our coefficients coolH(i) differ from ! those in Rosen et al. by a factor (mu mp)^2, with mu=1.2, since Rosen works ! in number density, n. (their cooling = Lambda*n^2, rho=mu mp n) ! The factor Lambdaunits converts from cgs units to code units. ! ! [Currently, coolT(1) is not modified, but this may be necessary ! to avoid creating gas too cold to resolve.] ! real, dimension (nx), intent(out) :: cool real, dimension (nx), intent(in) :: lnTT, lnrho integer :: i ! cool=0.0 cool_loop: do i=1,ncool if (lncoolT(i) >= lncoolT(i+1)) exit cool_loop where (lncoolT(i) <= lnTT .and. lnTT < lncoolT(i+1)) cool=cool+exp(lncoolH(i)+lnrho+lnTT*coolB(i)) endwhere enddo cool_loop endsubroutine calc_cool_func !***************************************************************************** subroutine calc_heat(heat,lnTT) ! ! This routine adds UV heating, cf. Wolfire et al., ApJ, 443, 152, 1995 ! with the values above, this gives about 0.012 erg/g/s (T < ~1.E4 K) ! Nb: need rho0 from density_[init/run]_pars, to implement the arm/interarm ! scaling. ! ! Control with heating_select in interstellar_init_pars/run_pars. ! Default heating_rate GammaUV = 0.015. ! real, dimension (nx), intent(out) :: heat real, dimension (nx), intent(in) :: lnTT ! ! Constant heating with a rate heating_rate[erg/g/s]. ! if (heating_select == 'cst') then heat = heating_rate_code else if (heating_select == 'wolfire') then heat(1:nx)=GammaUV*0.5*(1.0+tanh(cUV*(T0UV-exp(lnTT)))) else if (heating_select == 'wolfire_min') then heat(1:nx)=GammaUV*0.5*(1.0+tanh(cUV*(T0UV-exp(lnTT)))) heat = max(heat,heating_rate_code) ! ! If using thermal-hs in initial entropy this must also be specified for ! thermal equilibrium and applies for vertically stratified density supported ! by vertical gravity profile 'Ferriere'. ! else if (heating_select == 'thermal-hs') then heat(1:nx) = heat_z(n)*0.5*(1.0+tanh(cUV*(T0UV-exp(lnTT)))) else if (heating_select == 'off') then heat = 0. endif ! endsubroutine calc_heat !***************************************************************************** subroutine check_SN(f) ! ! Checks for SNe, and implements appropriately: ! relevant subroutines in entropy.f90 ! real, dimension(mx,my,mz,mfarray) :: f ! ! Only allow SNII if no SNI this step (may not be worth keeping). ! logical :: l_SNI=.false. ! intent(inout) :: f ! ! Identifier ! if (headtt) print*,'check_SN: ENTER' ! ! Do separately for SNI (simple scheme) and SNII (Boris' scheme). ! if (t < t_settle) return call calc_snr_damping_factor(f) call tidy_SNRs if (lSNI) call check_SNI (f,l_SNI) if (lSNII) call check_SNII(f,l_SNI) call calc_snr_damping_add_heat(f) ! endsubroutine check_SN !***************************************************************************** subroutine check_SNI(f,l_SNI) ! ! If time for next SNI, then implement, and calculate time of subsequent SNI. ! real, dimension(mx,my,mz,mfarray) :: f logical :: l_SNI integer :: try_count, iSNR, ierr ! intent(inout) :: f,l_SNI ! ! Identifier ! if (headtt) print*,'check_SNI: ENTER' ! l_SNI=.false. if (t >= t_next_SNI) then iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0.0 SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=1 SNRs(iSNR)%radius=width_SN try_count=10 ! do while (try_count>0) ierr=iEXPLOSION_OK try_count=try_count-1 ! if (uniform_zdist_SNI) then call position_SN_uniformz(f,SNRs(iSNR)) else call position_SN_gaussianz(f,h_SNI,SNRs(iSNR)) endif ! if (lforce_locate_SNI.and.(SNRs(iSNR)%site%rho < rho_SN_min).or. & (SNRs(iSNR)%site%TT > TT_SN_max)) then call find_nearest_SNI(f,SNRs(iSNR)) endif ! if (.not.lSN_scale_rad) then if ((SNRs(iSNR)%site%rho < rho_SN_min) .or. & (SNRs(iSNR)%site%TT > TT_SN_max)) then cycle endif else ! ! Avoid sites with dense mass shown to excessively cool remnants, given the ! high thermal conductivity we are forced to adopt. ! if (SNRs(iSNR)%site%rho > rho_SN_max) then cycle endif endif ! call explode_SN(f,SNRs(iSNR),ierr,preSN) if (ierr==iEXPLOSION_OK) then l_SNI=.true. if (lSN_mass_rate) then call set_interval(f,t_interval_SNI,l_SNI) endif call set_next_SNI exit endif enddo ! if (try_count==0) then if (lroot) print*, & "check_SNI: 10 RETRIES OCCURED - skipping SNI insertion" endif ! ! Free up slots in case loop fails repeatedly over many time steps. ! call free_SNR(iSNR) ! ! Reset ierr or else explode_SN may terminate erroneously on subsequent ! timesteps never to be reset. ! ierr=iEXPLOSION_OK endif ! endsubroutine check_SNI !***************************************************************************** subroutine check_SNIIb(f,l_SNI) ! ! If time for next SNI, then implement, and calculate time of subsequent SNI. ! real, dimension(mx,my,mz,mfarray) :: f integer :: try_count, iSNR, ierr logical :: l_SNI ! intent(inout) :: f, l_SNI ! ! Identifier ! if (headtt) print*,'check_SNIIb: ENTER' ! if (t >= t_next_SNII) then iSNR=get_free_SNR() SNRs(iSNR)%site%TT=1E20 SNRs(iSNR)%site%rho=0.0 SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=2 SNRs(iSNR)%radius=width_SN try_count=10 ! do while (try_count>0) ierr=iEXPLOSION_OK try_count=try_count-1 ! if (uniform_zdist_SNI) then call position_SN_uniformz(f,SNRs(iSNR)) else call position_SN_gaussianz(f,h_SNII,SNRs(iSNR)) endif ! if (lforce_locate_SNI.and.(SNRs(iSNR)%site%rho < rho_SN_min).or. & (SNRs(iSNR)%site%TT > TT_SN_max)) then call find_nearest_SNI(f,SNRs(iSNR)) endif ! if (.not.lSN_scale_rad) then if ((SNRs(iSNR)%site%rho < rho_SN_min) .or. & (SNRs(iSNR)%site%TT > TT_SN_max)) then cycle endif else ! ! Avoid sites with dense mass shown to excessively cool remnants, given the ! high thermal conductivity we are forced to adopt. ! if (SNRs(iSNR)%site%rho > rho_SN_max) then cycle endif endif ! call explode_SN(f,SNRs(iSNR),ierr,preSN) if (ierr==iEXPLOSION_OK) then if (lSN_mass_rate) then call set_interval(f,t_interval_SNI,l_SNI) endif call set_next_SNII exit endif enddo ! if (try_count==0) then if (lroot) print*, & "check_SNIIb: 10 RETRIES OCCURED - skipping SNII insertion" endif ! ! Free up slots in case loop fails repeatedly over many time steps. ! call free_SNR(iSNR) ! ! Reset ierr or else explode_SN may terminate erroneously on subsequent ! timesteps never to be reset. ! ierr=iEXPLOSION_OK endif l_SNI=.true. ! endsubroutine check_SNIIb !*********************************************************************** subroutine set_next_SNI() ! use General, only: random_number_wrapper ! real, dimension(1) :: franSN ! ! Pre-determine time for next SNI. ! if (lroot.and.ip<14) print*, & "check_SNI: Old t_next_SNI=", t_next_SNI call random_number_wrapper(franSN) ! ! Vary the time interval with a uniform random distribution between ! 0.8 and 1.2 times the average rate required. ! t_next_SNI=t + (1.0 + 0.4*(franSN(1)-0.5)) * t_interval_SNI if (lroot.and.ip<20) print*, & 'check_SNI: Next SNI at time = ' ,t_next_SNI ! endsubroutine set_next_SNI !***************************************************************************** subroutine set_next_SNII() ! use General, only: random_number_wrapper ! real, dimension(1) :: franSN ! ! Pre-determine time for next SNII ! Check_SNII has a selfregulating random rate governed by the parameters of ! cloud mass and cloud temperature, but this is very hard to regulate to test ! different regimes, so this acts as a contraint on the rate. ! if (lroot.and.ip<14) print*, & "check_SNII: Old t_next_SNII=", t_next_SNII call random_number_wrapper(franSN) ! ! Vary the time interval with a uniform random distribution between ! 0.4 and 1.6 times the average rate required. ! t_next_SNII=t + (1.0 + 1.2*(franSN(1)-0.5)) * t_interval_SNII if (lroot.and.ip<20) print*, & 'check_SNII: Next SNII at time = ' ,t_next_SNII ! endsubroutine set_next_SNII !***************************************************************************** subroutine set_interval(f,t_interval,l_SNI) ! use Mpicomm, only: mpireduce_sum, mpibcast_real use EquationOfState, only: getmu ! real, dimension(mx,my,mz,mfarray) :: f real :: t_interval, surface_massII, mu integer :: iz real, dimension(nx,ny,nz) :: disk_massII real :: MmpiII, msumtmpII logical :: l_SNI ! intent(IN) :: f, l_SNI intent(OUT) :: t_interval ! ! Identifier ! if (headtt) print*,'set_interval: ENTER' ! ! Adapt expected time interval until next SN depending on the ISM mass ! within 2x h_SN of midplane ! ! SNII rate=5.E-12 mass(H1+HII)/solar_mass ! van den Bergh/Tammann Annu. Rev Astron. Astrophys. 1991 29:363-407 ! SNI rate=4.7E-14 mass(H1+HII)/solar_mass + 0.35 x SNII rate ! Mannucci et al A&A 433, 807-814 (2005) ! if (ldensity_nolog) then disk_massII=f(l1:l2,m1:m2,n1:n2,irho) else disk_massII=exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) endif ! do iz=1,nz if (abs(z(iz+nghost))>2.0*h_SNII) disk_massII(1:nx,1:ny,iz)=0.0 enddo ! surface_massII=sum(disk_massII) msumtmpII=surface_massII call mpireduce_sum(msumtmpII,MmpiII) call mpibcast_real(MmpiII) surface_massII=MmpiII*dv ! call getmu(f,mu) ! if (l_SNI) then ! t_interval=solar_mass/(SNI_mass_rate+0.35*SNII_mass_rate)/ & ! surface_massII/mu t_interval=7.5*solar_mass/SNII_mass_rate/surface_massII/mu if (lroot.and.ip<20) print*, & 'set_interval: expected interval for SNI =',t_interval else t_interval=solar_mass/surface_massII/SNII_mass_rate/mu if (lroot.and.ip<20) print*, & 'set_interval: expected interval for SNII =',t_interval endif ! endsubroutine set_interval !***************************************************************************** subroutine check_SNII(f,l_SNI) ! ! Check for SNII, via self-regulating scheme. ! ! 03-feb-10/fred: Tested and working correctly. ! use General, only: random_number_wrapper use Mpicomm, only: mpireduce_sum, mpibcast_real use EquationOfState, only: eoscalc, ilnrho_ss, irho_ss ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(nx) :: rho, rho_cloud, lnTT, TT, yH real :: cloud_mass, cloud_mass_dim, freq_SNII, prob_SNII real :: franSN, fsum1, fsum1_tmp, fmpi1 real, dimension(ncpus) :: cloud_mass_byproc integer :: icpu, m, n, iSNR, ierr logical :: l_SNI real :: dtsn ! intent(inout) :: f,l_SNI ! ! Identifier ! if (lroot.and.headtt.and.ip<14) print*,'check_SNII: ENTER' ! if (l_SNI) return ! Only do if no SNI this step. ! if (lSNII_gaussian) then ! Skip location by mass. call check_SNIIb(f,l_SNI) return endif ! iSNR=get_free_SNR() ! ! Determine and sum all cells comprising dense cooler clouds where type II ! SNe are prevalent. Only check if t_next_SNII exceeded (see set_next_SNII). ! if (t >= t_next_SNII) then cloud_mass=0.0 ! ! Calculate the total mass in locations where the temperature is below ! cloud_TT and the density is above cloud_rho, i.e. cold and dense. ! do n=n1,n2 do m=m1,m2 if (ldensity_nolog) then rho(1:nx)=f(l1:l2,m,n,irho) call eoscalc(irho_ss,f(l1:l2,m,n,irho),f(l1:l2,m,n,iss)& ,yH=yH,lnTT=lnTT) else rho(1:nx)=exp(f(l1:l2,m,n,ilnrho)) call eoscalc(ilnrho_ss,f(l1:l2,m,n,ilnrho),f(l1:l2,m,n,iss)& ,yH=yH,lnTT=lnTT) endif TT(1:nx)=exp(lnTT(1:nx)) rho_cloud(1:nx)=0.0 where (rho(1:nx) >= cloud_rho .and. TT(1:nx) <= cloud_TT) & rho_cloud(1:nx) = rho(1:nx) cloud_mass=cloud_mass+sum(rho_cloud(1:nx)) enddo enddo ! ! Sum the total over all processors and multiply by dv to find total mass. ! fsum1_tmp=cloud_mass call mpireduce_sum(fsum1_tmp,fsum1) call mpibcast_real(fsum1) cloud_mass_dim=fsum1*dv ! if (ip<14) print*, & 'check_SNII: cloud_mass,it,iproc=',cloud_mass,it,iproc ! if (lroot .and. ip < 14) & print*, 'check_SNII: cloud_mass_dim,fsum(1),dv:', & cloud_mass_dim,fsum1,dv ! ! Additional contraint on the interval between SNII events. The total time ! elapsed since last SNII is dtsn. Probability of next event increases with ! time and availabilty of cold dense cloud material. Calculate probability. ! (SNI distribution is independent of mass distribution.) ! This probability is close to 1 for solar neighbourhood rate so this check ! could be discarded, but worth keeping as may be interesting tool. ! dtsn=t-last_SN_t freq_SNII=frac_heavy*frac_converted*cloud_mass_dim/ & mass_SN_progenitor/cloud_tau prob_SNII=freq_SNII*dtsn*5. call random_number_wrapper(franSN) ! if (lroot.and.ip<20) then if (cloud_mass_dim>0.0.and.franSN<=2.0*prob_SNII) then print*,'check_SNII: freq,prob,rnd,dtsn:', & freq_SNII,prob_SNII,franSN,dtsn print*,'check_SNII: frac_heavy,frac_converted,cloud_mass_dim,', & 'mass_SN,cloud_tau',& frac_heavy,frac_converted,cloud_mass_dim,mass_SN,cloud_tau endif endif ! ! If likelihood of SNII greater than random number locate SNII. ! if (franSN <= prob_SNII) then ! ! The position_SN_bycloudmass needs the cloud_masses for each processor. ! Communicate and store them here, to avoid recalculation. ! cloud_mass_byproc(:)=0.0 ! ! Use non-root broadcasts for the communication... ! do icpu=1,ncpus fmpi1=cloud_mass call mpibcast_real(fmpi1,icpu-1) cloud_mass_byproc(icpu)=fmpi1 enddo ! ! Locate the next explosion. ! if (lroot.and.ip<14) print*, & 'check_SNII: cloud_mass_byproc:',cloud_mass_byproc call position_SN_bycloudmass& (f,cloud_mass_byproc,SNRs(iSNR),preSN,ierr) ! ! If location too hot reset ierr and return to program. ! if (ierr == iEXPLOSION_TOO_HOT) then call free_SNR(iSNR) ierr=iEXPLOSION_OK return endif ! ! Try to explode SNII and if successful reset time of most recent (last_SN_t) ! and next (t_next_SNII) explosion. ! SNRs(iSNR)%t=t SNRs(iSNR)%SN_type=2 call explode_SN(f,SNRs(iSNR),ierr,preSN) if (ierr==iEXPLOSION_OK) then if (lSN_mass_rate) then call set_interval(f,t_interval_SNII,l_SNI) endif call set_next_SNII last_SN_t=t endif ! endif endif ! ! If returned unexploded stop code running out of free slots & reset ierr. ! ierr=iEXPLOSION_OK call free_SNR(iSNR) l_SNI=.true. ! endsubroutine check_SNII !***************************************************************************** subroutine position_SN_testposition(f,SNR) ! ! Determine position for next SN (w/ fixed scale-height). ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant), intent(inout) :: SNR ! real :: z00, x00, y00 integer :: i ! if (headtt) print*,'position_SN_testposition: ENTER' ! ! Calculate the global (nzgrid) lower z-coordinate. ! if (lperi(1)) then; x00=xyz0(1)-.5*dx; else; x00=xyz0(1); endif if (lperi(2)) then; y00=xyz0(2)-.5*dy; else; y00=xyz0(2); endif if (lperi(3)) then; z00=xyz0(3)-.5*dz; else; z00=xyz0(3); endif ! ! Pick SN position (SNR%l,SNR%m,SNR%n). ! if (lroot) then if (center_SN_x==impossible) then i=max(int(nxgrid/2)+1,1) else i=int((center_SN_x-x00)/dx)+1 endif SNR%l=i+nghost ! if (center_SN_y==impossible) then i=max(int(nygrid/2)+1,1) else i=int((center_SN_y-y00)/dy)+1 endif SNR%ipy=(i-1)/ny ! uses integer division SNR%m=i-(SNR%ipy*ny)+nghost ! if (center_SN_z==impossible) then i=max(int(nzgrid/2)+1,1) else i=int((center_SN_z-z00)/dz)+1 endif SNR%ipz=(i-1)/nz ! uses integer division SNR%n=i-(SNR%ipz*nz)+nghost SNR%iproc=SNR%ipz*nprocy + SNR%ipy endif call share_SN_parameters(f,SNR) ! endsubroutine position_SN_testposition !***************************************************************************** subroutine position_SN_gaussianz(f,h_SN,SNR) ! ! Determine position for next SN (w/ fixed scale-height). ! 15-jul-14/luiz+fred: added SN horizontal clustering algorithm for SNII ! 70% probability SNII located within 300pc of previous ! use General, only: random_number_wrapper use Mpicomm, only: mpiallreduce_max, mpireduce_min, mpireduce_max,& mpibcast_real ! real, intent(in), dimension(mx,my,mz,mfarray) :: f real, intent(in) :: h_SN type (SNRemnant), intent(inout) :: SNR ! ! parameters required to determine the vertical centre of mass of the disk ! real, dimension(nprocz) :: tmp3 real, dimension(nz) :: rhotmp real :: zdisk, rhomax, maxrho, rhosum real :: mpirho, mpiz real, dimension(ncpus):: tmp2 integer :: yzproc, itmp, icpu, lm_range, previous_SNl, previous_SNm ! ! parameters for random location of SN - about zdisk ! real, dimension(nzgrid) :: cum_prob_SN real :: zn, z00, x00, y00 real, dimension(3) :: fran3 integer :: i, nzskip=10 !prevent SN from being too close to boundaries ! if (headtt) print*,'position_SN_gaussianz: ENTER' ! ! Calculate the global (nzgrid) lower z-coordinate. ! if (lperi(1)) then; x00=xyz0(1)+.5*dx; else; x00=xyz0(1); endif if (lperi(2)) then; y00=xyz0(2)+.5*dy; else; y00=xyz0(2); endif if (lperi(3)) then; z00=xyz0(3)+.5*dz; else; z00=xyz0(3); endif ! ! The disk oscillates. to keep the random dist centred at the disk find ! zmode where the peak mean density(z) resides and shift gaussian up/down ! rhomax=0.0 zdisk=0.0 ! if (h_SN==h_SNII) then if (ldensity_nolog) then rhosum=sum(f(l1:l2,m1:m2,n1:n2,irho)) else rhosum=sum(exp(f(l1:l2,m1:m2,n1:n2,ilnrho))) endif ! do icpu=1,ncpus mpirho=rhosum call mpibcast_real(mpirho,icpu-1) tmp2(icpu)=mpirho enddo ! do i=1,nprocz tmp3(i)=sum(tmp2((i-1)*nprocy+1:i*nprocy)) enddo ! rhomax=maxval(tmp3) do i=1,nprocz if (tmp3(i)==rhomax) itmp=(i-1)*nprocy enddo rhomax=maxval(tmp2(itmp+1:itmp+nprocy)) do icpu=1,nprocy if (tmp2(icpu+itmp) == rhomax) yzproc=icpu+itmp-1 enddo ! if (iproc==yzproc) then rhomax=0.0 do i=n1,n2 if (ldensity_nolog) then rhotmp(i-nghost)=sum(f(l1:l2,m1:m2,i,irho)) else rhotmp(i-nghost)=sum(exp(f(l1:l2,m1:m2,i,ilnrho))) endif maxrho=max(rhomax,rhotmp(i-nghost)) rhomax=maxrho if (rhotmp(i-nghost) == rhomax) zdisk = z(i) enddo mpiz=zdisk endif call mpibcast_real(mpiz,yzproc) zdisk = mpiz endif ! ! Pick SN position (SNR%l,SNR%m,SNR%n). ! call random_number_wrapper(fran3) ! ! Get 3 random numbers on all processors to keep rnd. generators in sync. ! if (lroot) then if (lOB_cluster .and. h_SN==h_SNII) then previous_SNl = int(( x_next_SNII - xyz0(1) )/Lxyz(1))*nxgrid +1 previous_SNm = int(( y_next_SNII - xyz0(2) )/Lxyz(2))*nygrid +1 lm_range = 2*SN_clustering_radius*nxgrid/Lxyz(1) if (fran3(1) < p_OB) then ! checks whether the SN is in a cluster i=int(fran3(1)*lm_range/p_OB)+previous_SNl+1 SNR%ipx=(i-1)/nx ! uses integer division SNR%l=i-(SNR%ipx*nx)+nghost ! i=int(fran3(1)*lm_range/p_OB)+previous_SNm+1 SNR%ipy=(i-1)/ny ! uses integer division SNR%m=i-(SNR%ipy*ny)+nghost else ! outside cluster i=int(fran3(1)*(nxgrid-lm_range)/(1.0-p_OB))+previous_SNl+1 if (i>nxgrid) i=i-nxgrid SNR%ipx=(i-1)/nx ! uses integer division SNR%l=i-(SNR%ipx*nx)+nghost ! i=int(fran3(1)*(nygrid-lm_range)/(1.0-p_OB))+previous_SNl+1 if (i>nygrid) i=i-nygrid SNR%ipy=(i-1)/ny ! uses integer division SNR%m=i-(SNR%ipy*ny)+nghost endif else ! clustering not used i=int(fran3(1)*nxgrid)+1 SNR%ipx=(i-1)/nx ! uses integer division SNR%l=i-(SNR%ipx*nx)+nghost ! i=int(fran3(2)*nygrid)+1 SNR%ipy=(i-1)/ny ! uses integer division SNR%m=i-(SNR%ipy*ny)+nghost endif ! ! Cumulative probability function in z calculated each time for moving zdisk. ! print*,'position_SN_gaussianz: zdisk =',zdisk cum_prob_SN=0.0 do i=nzskip+1,nzgrid-nzskip zn=z00+(i-1)*dz cum_prob_SN(i)=cum_prob_SN(i-1)+exp(-0.5*((zn-zdisk)/h_SN)**2) enddo cum_prob_SN = cum_prob_SN / max(cum_prob_SN(nzgrid-nzskip), tini) ! ! The following should never be needed, but just in case floating point ! errors ever lead to cum_prob_SNI(nzgrid-nzskip) < rnd < 1. ! cum_prob_SN(nzgrid-nzskip+1:nzgrid)=1.0 ! do i=nzskip+1,nzgrid-nzskip if (cum_prob_SN(i-1)<=fran3(3) .and. fran3(3)rho_SN_min).and.(SNR%site%TT < TT_SN_max)) then deltar2_test=((ii-SNR%l)**2+(m-SNR%m)**2+(n-SNR%n)**2) if (deltar2_test < deltar2) then nfound=1 deltar2=deltar2_test new_lmn=(/ nfound, ii, m, n /) elseif (deltar2==deltar2_test) then nfound=nfound+1 endif endif enddo ! enddo enddo ! if (nfound==0) then new_lmn=(/ nfound, SNR%l, SNR%m, SNR%n /) elseif (nfound>1) then chosen_site=int(nfound*fran_location(1)+0.5) nfound=0 search_two: do n=n1,n2 do m=m1,m2 if (ldensity_nolog) then rho_test=f(l1:l2,m,n,irho) else rho_test=exp(f(l1:l2,m,n,ilnrho)) endif call eoscalc(f,nx,lnTT=lnTT_test) TT_test=exp(lnTT_test) ! do ii=l1,l2 if ((SNR%site%rho>rho_SN_min).and.(SNR%site%TT0) then SNR%l=new_lmn(2) SNR%m=new_lmn(3) SNR%n=new_lmn(4) call share_SN_parameters(f,SNR) endif ! endsubroutine find_nearest_SNI !***************************************************************************** subroutine share_SN_parameters(f,SNR) ! ! Handle common SN positioning processor communications. ! ! 27-aug-2003/tony: coded ! use EquationOfState, only: eoscalc,ilnrho_lnTT use Mpicomm, only: mpibcast_int, mpibcast_real ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant), intent(inout) :: SNR ! real, dimension(nx) :: lnTT real, dimension(6) :: fmpi5 integer, dimension(4) :: impi4 ! ! Broadcast position to all processors from root; also broadcast SNR%iproc, ! needed for later broadcast of SNR%site%rho. ! impi4=(/ SNR%iproc, SNR%l, SNR%m, SNR%n /) call mpibcast_int(impi4,4) SNR%iproc=impi4(1) SNR%l=impi4(2) SNR%m=impi4(3) SNR%n=impi4(4) ! ! With current SN scheme, we need rho at the SN location. ! if (iproc==SNR%iproc) then if (ldensity_nolog) then SNR%site%lnrho=log(f(SNR%l,SNR%m,SNR%n,irho)) else SNR%site%lnrho=f(SNR%l,SNR%m,SNR%n,ilnrho) endif SNR%site%rho=exp(SNR%site%lnrho); ! ! 10-Jun-10/fred: ! Adjust radius according to density of explosion site to concentrate energy ! when locations are dense. ! SNR%radius=width_SN if (lSN_scale_rad) & SNR%radius=(0.75*solar_mass/SNR%site%rho*pi_1*N_mass)**(1.0/3.0) SNR%radius=max(SNR%radius,1.75*dxmax) ! minimum grid resolution ! m=SNR%m n=SNR%n call eoscalc(f,nx,lnTT=lnTT) SNR%site%lnTT=lnTT(SNR%l-l1+1) SNR%x=0.; SNR%y=0.; SNR%z=0. if (nxgrid/=1) SNR%x=x(SNR%l) +0.5*dx*(-1.)**SNR%l if (nygrid/=1) SNR%y=y(SNR%m) +0.5*dy*(-1.)**SNR%m if (nzgrid/=1) SNR%z=z(SNR%n) +0.5*dz*(-1.)**SNR%n ! ! Better initialise these to something on the other processors ! else SNR%site%lnrho=0. SNR%site%lnTT=0. SNR%x=0. SNR%y=0. SNR%z=0. endif ! ! Broadcast to all processors. ! fmpi5=(/ SNR%x, SNR%y, SNR%z, SNR%site%lnrho, SNR%site%lnTT, SNR%radius /) call mpibcast_real(fmpi5,6,SNR%iproc) ! SNR%x=fmpi5(1); SNR%y=fmpi5(2); SNR%z=fmpi5(3); SNR%site%lnrho=fmpi5(4); SNR%site%lnTT=fmpi5(5); SNR%radius=fmpi5(6) ! SNR%site%rho=exp(SNR%site%lnrho); ! call eoscalc(ilnrho_lnTT,SNR%site%lnrho,SNR%site%lnTT, & yH=SNR%site%yH,ss=SNR%site%ss,ee=SNR%site%ee) SNR%site%TT=exp(SNR%site%lnTT) ! if (lroot.and.ip<24) print*, & 'share_SN_parameters: SNR%iproc,x_SN,y_SN,z_SN,SNR%l,SNR%m,SNR%n,=', & SNR%iproc,SNR%x,SNR%y,SNR%z,SNR%l,SNR%m,SNR%n if (lroot.and.ip<24) print*, & 'share_SN_parameters: SNR%site%rho,SNR%site%ss,SNR%site%TT,SNR%radius=', & SNR%site%rho,SNR%site%ss,SNR%site%TT,SNR%radius ! endsubroutine share_SN_parameters !***************************************************************************** subroutine explode_SN(f,SNR,ierr,preSN) ! ! Implement SN (of either type), at pre-calculated position. ! ! ??-nov-02/grs : coded from GalaxyCode ! 20-may-03/tony: pencil formulation and broken into subroutines ! use EquationOfState, only: ilnrho_ee, eoscalc, getdensity, eosperturb ,& ilnrho_ss, irho_ss use Mpicomm, only: mpireduce_max, mpibcast_real, & mpireduce_sum, mpibcast_int use General, only: keep_compiler_quiet ! real, intent(inout), dimension(mx,my,mz,mfarray) :: f type (SNRemnant), intent(inout) :: SNR integer, intent(inout), optional, dimension(4,npreSN) :: preSN integer, optional :: ierr ! real :: c_SN,cmass_SN,cvelocity_SN real :: mass_shell real :: rho_SN_lowest real :: width_energy, width_mass, width_velocity real :: cavity_depth, r_cavity, rhom, ekintot real :: rhom_new, ekintot_new real :: rho_SN_new,lnrho_SN_new,yH_SN_new,lnTT_SN_new,ee_SN_new real :: TT_SN_new, uu_sedov ! real, dimension(nx) :: deltarho, deltaEE real, dimension(nx,3) :: deltauu real, dimension(2) :: dmpi2, dmpi2_tmp real, dimension(nx) :: lnrho, yH, lnTT, TT, rho_old, ee_old, site_rho real, dimension(nx,3) :: uu real :: maxlnTT, site_mass, maxTT=0.,mmpi, mpi_tmp real :: radiusA, radiusB, t_interval_SN integer :: i ! logical :: lmove_mass=.false. ! SNR%state=SNstate_exploding ! ! identifier ! if (lroot.and.ip<12) print*,'explode_SN: SN type =',SNR%SN_type ! ! Calculate explosion site mean density. ! call get_properties(f,SNR,rhom,ekintot) SNR%rhom=rhom ! ! Rescale injection radius by mass if required. Iterate a few times to ! improve match of mass to radius. ! if (lSN_scale_rad) then do i=1,20 SNR%radius=(0.75*solar_mass/SNR%rhom*pi_1*N_mass)**(1.0/3.0) SNR%radius=max(SNR%radius,1.75*dxmax) call get_properties(f,SNR,rhom,ekintot) SNR%rhom=rhom enddo SNR%radius=(0.75*solar_mass/SNR%rhom*pi_1*N_mass)**(1.0/3.0) SNR%radius=max(SNR%radius,1.75*dxmax) endif call get_properties(f,SNR,rhom,ekintot) SNR%rhom=rhom ! ! Calculate effective Sedov evolution time diagnostic and used in damping. ! SNR%t_sedov=sqrt((SNR%radius/xsi_sedov)**5*SNR%rhom/(kampl_SN+ampl_SN)) uu_sedov = 0.4*SNR%radius/SNR%t_sedov ! ! This may no longer be required. The appropriate radial velocity is now ! calculated using cvelocity_SN from kinetic energy injection kampl_SN. ! if ((uu_sedov_max > 0.).and.(uu_sedov > uu_sedov_max)) then do i=1,10 radiusA=SNR%radius radiusB=(0.16/uu_sedov_max**2*xsi_sedov**5* & (kampl_SN+ampl_SN)/SNR%rhom)**(1./3.) ! if (abs(radiusB-radiusA) < dxmax) then SNR%radius=max(radiusA,radiusB) call get_properties(f,SNR,rhom,ekintot) SNR%rhom=rhom SNR%t_sedov = sqrt((SNR%radius/xsi_sedov)**5* & SNR%rhom/(kampl_SN+ampl_SN)) uu_sedov = 0.4*SNR%radius/SNR%t_sedov exit endif ! SNR%radius=0.5*(radiusA+radiusB) call get_properties(f,SNR,rhom,ekintot) SNR%rhom=rhom SNR%t_sedov = sqrt((SNR%radius/xsi_sedov)**5* & SNR%rhom/(kampl_SN+ampl_SN)) uu_sedov = 0.4*SNR%radius/SNR%t_sedov enddo if (SNR%radius>2*width_SN) then if (present(ierr)) then ierr=iEXPLOSION_TOO_RARIFIED endif return endif if (lroot.and.ip<14) print*,"explode_SN: Tweaked width ",SNR%radius endif ! width_energy = SNR%radius*energy_width_ratio width_mass = SNR%radius*mass_width_ratio width_velocity = SNR%radius*velocity_width_ratio ! ! Energy insertion normalization. ! if (thermal_profile=="gaussian3") then c_SN=ampl_SN/(cnorm_SN(dimensionality)*width_energy**dimensionality) ! elseif (thermal_profile=="gaussian2") then c_SN=ampl_SN/(cnorm_gaussian2_SN(dimensionality)* & width_energy**dimensionality) ! elseif (thermal_profile=="gaussian") then c_SN=ampl_SN/(cnorm_gaussian_SN(dimensionality)* & width_energy**dimensionality) ! elseif (thermal_profile=="quadratic") then c_SN=ampl_SN/(cnorm_para_SN(dimensionality)* & width_energy**dimensionality) ! elseif (thermal_profile=="quadratictanh") then c_SN=ampl_SN/(cnorm_para_SN(dimensionality)* & width_energy**dimensionality) ! elseif (thermal_profile=="quartictanh") then c_SN=ampl_SN/(cnorm_quar_SN(dimensionality)* & width_energy**dimensionality) ! elseif (thermal_profile=="tanh") then if (dimensionality==1) then c_SN=ampl_SN/( 2.*width_energy ) elseif (dimensionality==2) then c_SN=ampl_SN/( pi*(width_energy)**2 ) elseif (dimensionality==3) then c_SN=ampl_SN/( 4./3.*pi*(width_energy)**3 ) endif endif ! if (lroot.and.ip<14) print*,'explode_SN: c_SN =',c_SN ! ! Mass insertion normalization. ! if (lSN_mass) then if (mass_profile=="gaussian3") then cmass_SN=mass_SN/(cnorm_SN(dimensionality)* & width_mass**dimensionality) ! elseif (mass_profile=="gaussian2") then cmass_SN=mass_SN/(cnorm_gaussian2_SN(dimensionality)* & width_mass**dimensionality) ! elseif (mass_profile=="gaussian") then cmass_SN=mass_SN/(cnorm_gaussian_SN(dimensionality)* & width_mass**dimensionality) ! elseif (mass_profile=="quadratic") then cmass_SN=mass_SN/(cnorm_para_SN(dimensionality)* & width_mass**dimensionality) ! elseif (mass_profile=="tanh") then if (dimensionality==1) then cmass_SN=mass_SN/( 2.*width_mass ) elseif (dimensionality==2) then cmass_SN=mass_SN/( pi*(width_mass)**2 ) elseif (dimensionality==3) then cmass_SN=mass_SN/( 4./3.*pi*(width_mass)**3 ) endif endif ! if (lroot.and.ip<14) print*,'explode_SN: cmass_SN =',cmass_SN else cmass_SN=0. endif ! ! Calculate cross over point between mass addition and removal if mass ! movement is used. ! r_cavity = width_energy* & (dimensionality*log(outer_shell_proportion/inner_shell_proportion)/& ((1./inner_shell_proportion**6)- & (1./outer_shell_proportion**6)))**(1./6.) if (lroot.and.ip<14) print*, & 'explode_SN: dimensionality,r_cavity',dimensionality,r_cavity if (lroot.and.ip<14) print*,'explode_SN: shell_(inner, outer)_prop.=', & inner_shell_proportion,outer_shell_proportion if (lroot.and.ip<14) print*, & 'explode_SN: width_energy,c_SN,SNR%site%rho=', & width_energy,c_SN,SNR%site%rho ! ! Now deal with (if nec.) mass relocation ! if (lroot.and.ip<14) print*, & 'explode_SN: rho_new,SNR%site%ee=',SNR%site%rho,SNR%site%ee ee_SN_new = (SNR%site%ee+frac_eth*c_SN/(SNR%site%rho+cmass_SN)) if (lroot.and.ip<14) print*, & 'explode_SN: rho_SN_new,ee_SN_new=',SNR%site%rho+cmass_SN,ee_SN_new call eoscalc(ilnrho_ee,real(log(SNR%site%rho+cmass_SN)),ee_SN_new, & lnTT=lnTT_SN_new,yH=yH_SN_new) TT_SN_new=exp(lnTT_SN_new) ! ! Velocity insertion normalization. ! 26-aug-10/fred: ! E_k=int(0.5*rho*vel^2)=approx 2pi*rhom*V0^2*int(r^2*v(r)dr). ! Total energy = kinetic (kampl_SN) + thermal (ampl_SN). ! if (lSN_velocity) then if (velocity_profile=="gaussian3") then cvelocity_SN= & sqrt(kampl_SN*pi_1/SNR%rhom/0.4177713791/width_velocity**3) ! elseif (velocity_profile=="gaussian2") then cvelocity_SN= & sqrt(kampl_SN*pi_1/SNR%rhom/0.3643185655/width_velocity**3) ! elseif (velocity_profile=="gaussian") then cvelocity_SN= & sqrt(kampl_SN*pi_1/SNR%rhom/0.313328534/width_velocity**3) ! elseif (velocity_profile=="r16thgaussian") then cvelocity_SN= sqrt & (kampl_SN*pi_1/SNR%rhom/0.3012691725/width_velocity**3.125) !25/8 ! elseif (velocity_profile=="r16thgaussian3") then cvelocity_SN=sqrt & (kampl_SN/pi/SNR%rhom/0.3956910052/width_velocity**3.125) !25/8 ! else cvelocity_SN=uu_sedov if (lroot) print*, & 'explode_SN: cvelocity_SN is uu_sedov, velocity profile = ', & velocity_profile endif if (lroot.and.ip<14) print*,'explode_SN: cvelocity_SN =',cvelocity_SN ! else cvelocity_SN=0. endif ! if (lroot.and.ip<14) print*, & 'explode_SN: SNR%site%TT, TT_SN_new, TT_SN_min, SNR%site%ee =', & SNR%site%TT,TT_SN_new,TT_SN_min, SNR%site%ee if (lroot.and.ip<14) print*,'explode_SN: yH_SN_new =',yH_SN_new ! if ((TT_SN_new < TT_SN_min).or.(mass_movement=='constant')) then if (lroot.and.ip<20) print*,'explode_SN: SN will be too cold!' ! lmove_mass=.not.(mass_movement == 'off') ! lmove_mass=.false. ! use to switch off for debug... ! ! The bit that BREAKS the pencil formulation... ! Must know the total moved mass BEFORE attempting mass relocation. ! ! ASSUME: SN will fully ionize the gas at its centre ! if (lmove_mass) then if (lroot.and.ip<16) print*,'explode_SN: moving mass to compensate.' call getdensity(real((SNR%site%ee*SNR%site%rho)+frac_eth*c_SN), & TT_SN_min,1.,rho_SN_new) if (mass_movement=='rho-cavity') then call get_lowest_rho(f,SNR,r_cavity,rho_SN_lowest) cavity_depth=SNR%site%rho-rho_SN_new if (cavity_depth > rho_SN_lowest-rho_min) then cavity_depth=rho_SN_lowest-rho_min if (cavity_depth <= 0.) then cavity_depth=0. lmove_mass=.false. endif if (lroot.and.ip<16) print*,"Reduced cavity from:,", & SNR%site%rho-rho_SN_new," to: ",cavity_depth rho_SN_new=SNR%site%rho-cavity_depth lnrho_SN_new=log(rho_SN_new) endif elseif (mass_movement=='Galaxycode') then lnrho_SN_new=log(rho_SN_new-cmass_SN) cavity_depth=max(SNR%site%lnrho-lnrho_SN_new,0.) cavity_profile="gaussian3log" elseif (mass_movement=='constant') then lnrho_SN_new=log(cmass_SN) cavity_depth=cmass_SN cavity_profile="tanh" endif endif ! if (lmove_mass) then ee_SN_new=(SNR%site%ee*SNR%site%rho+frac_eth*c_SN)/rho_SN_new ! call eoscalc(ilnrho_ee,lnrho_SN_new,ee_SN_new, & lnTT=lnTT_SN_new,yH=yH_SN_new) TT_SN_new=exp(lnTT_SN_new) ! if (lroot.and.ip<16) print*, & 'explode_SN: Relocate mass... TT_SN_new, rho_SN_new=', & TT_SN_new,rho_SN_new ! if (mass_movement=='rho_cavity') then ! ! Do nothing. ! elseif (mass_movement=='Galaxycode') then call calc_cavity_mass_lnrho & (f,SNR,width_energy,cavity_depth,mass_shell) if (lroot.and.ip<16) print*,'explode_SN: mass_shell=',mass_shell elseif (mass_movement=='constant') then call calc_cavity_mass_lnrho & (f,SNR,width_energy,cavity_depth,mass_shell) if (lroot.and.ip<16) print*,'explode_SN: mass_shell=',mass_shell endif endif endif ! ! Validate the explosion. ! site_mass=0.0 maxlnTT=-10.0 do n=n1,n2 do m=m1,m2 SNR%state=SNstate_waiting ! ! Calculate the distances to the SN origin for all points in the current ! pencil and store in the dr2_SN global array. ! call proximity_SN(SNR) ! ! Get the old energy ! if (ldensity_nolog) then lnrho=log(f(l1:l2,m,n,irho)) rho_old=exp(lnrho) else lnrho=f(l1:l2,m,n,ilnrho) rho_old=exp(lnrho) endif site_rho=rho_old ! ! Calculate the ambient mass for the remnant. ! where (dr2_SN > SNR%radius**2.0) site_rho = 0.0 site_mass=site_mass+sum(site_rho) deltarho=0. ! call eoscalc(irho_ss,rho_old,f(l1:l2,m,n,iss),& yH=yH,lnTT=lnTT,ee=ee_old) TT=exp(lnTT) ! ! Apply perturbations ! call injectenergy_SN(deltaEE,width_energy,c_SN,SNR%EE) ! if (lmove_mass) then if (mass_movement=='rho_cavity') then if (lSN_mass) then call make_cavity_rho(deltarho,width_energy,cavity_depth, & cnorm_SN(dimensionality),SNR%MM) else call make_cavity_rho(deltarho,width_energy,cavity_depth, & cnorm_SN(dimensionality),SNR%MM) endif lnrho=log(rho_old(1:nx)+deltarho(1:nx)) elseif (mass_movement=='Galaxycode') then if (lSN_mass) then call make_cavity_lnrho(lnrho,width_energy,cavity_depth, & (mass_shell+mass_SN),cnorm_SN(dimensionality),SNR%MM) else call make_cavity_lnrho(lnrho,width_energy,cavity_depth, & mass_shell,cnorm_SN(dimensionality),SNR%MM) endif elseif (mass_movement=='constant') then call make_cavity_lnrho(lnrho,width_mass,cmass_SN, & mass_shell,cnorm_SN(dimensionality),SNR%MM) endif else if (lSN_mass) then call injectmass_SN(deltarho,width_mass,cmass_SN,SNR%MM) lnrho=log(rho_old(1:nx)+deltarho(1:nx)) endif endif ! if (lSN_eth) then call eoscalc(ilnrho_ee,lnrho,real( & (ee_old*rho_old+deltaEE*frac_eth)/exp(lnrho)), lnTT=lnTT) where (dr2_SN > SNR%radius**2.0) lnTT=-10.0 maxTT=maxval(exp(lnTT)) maxlnTT=max(log(maxTT),maxlnTT) call mpibcast_real(maxlnTT,SNR%iproc) maxTT=exp(maxlnTT) ! ! Broadcast maxlnTT from remnant to all processors so all take the same path ! after these checks. ! if (maxTT>TT_SN_max) then if (present(ierr)) then ierr=iEXPLOSION_TOO_HOT endif return endif ! endif enddo enddo ! if (present(ierr)) then call mpibcast_int(ierr,SNR%iproc) if (ierr==iEXPLOSION_TOO_HOT) return endif ! if (lSN_velocity)then call get_props_check(f,SNR,rhom,ekintot_new,cvelocity_SN,cmass_SN) if (ekintot_new-ekintot > 2.5*kampl_SN) then if (present(ierr)) then ierr=iEXPLOSION_TOO_UNEVEN endif endif endif ! if (present(ierr)) then call mpibcast_int(ierr,SNR%iproc) if (ierr==iEXPLOSION_TOO_UNEVEN) return endif ! mpi_tmp=site_mass call mpireduce_sum(mpi_tmp,mmpi) call mpibcast_real(mmpi) site_mass=mmpi*dv ! SNR%EE=0. SNR%MM=0. !EE_SN2=0. do n=n1,n2 do m=m1,m2 ! ! Calculate the distances to the SN origin for all points in the current ! pencil and store in the dr2_SN global array. ! call proximity_SN(SNR) ! Get the old energy. if (ldensity_nolog) then lnrho=log(f(l1:l2,m,n,irho)) rho_old=exp(lnrho) else lnrho=f(l1:l2,m,n,ilnrho) rho_old=exp(lnrho) endif deltarho=0. ! call eoscalc(irho_ss,rho_old,f(l1:l2,m,n,iss),& yH=yH,lnTT=lnTT,ee=ee_old) TT=exp(lnTT) ! ! Apply perturbations. ! call injectenergy_SN(deltaEE,width_energy,c_SN,SNR%EE) if (lmove_mass) then if (mass_movement=='rho_cavity') then if (lSN_mass) then call make_cavity_rho(deltarho,width_energy,cavity_depth, & cnorm_SN(dimensionality),SNR%MM) else call make_cavity_rho(deltarho,width_energy,cavity_depth, & cnorm_SN(dimensionality),SNR%MM) endif lnrho=log(rho_old(1:nx)+deltarho(1:nx)) elseif (mass_movement=='Galaxycode') then if (lSN_mass) then call make_cavity_lnrho(lnrho,width_energy,cavity_depth, & (mass_shell+mass_SN),cnorm_SN(dimensionality),SNR%MM) else call make_cavity_lnrho(lnrho,width_energy,cavity_depth, & mass_shell,cnorm_SN(dimensionality),SNR%MM) endif elseif (mass_movement=='constant') then call make_cavity_lnrho(lnrho,width_mass,cmass_SN, & mass_shell,cnorm_SN(dimensionality),SNR%MM) endif else if (lSN_mass) then call injectmass_SN(deltarho,width_mass,cmass_SN,SNR%MM) lnrho=log(rho_old(1:nx)+deltarho(1:nx)) endif endif ! if (lSN_velocity) then uu=f(l1:l2,m,n,iux:iuz) call injectvelocity_SN(deltauu,width_velocity,cvelocity_SN) f(l1:l2,m,n,iux:iuz)=uu+deltauu endif ! TT=exp(lnTT) ! if (lcosmicray.and.lSN_ecr) then f(l1:l2,m,n,iecr) = f(l1:l2,m,n,iecr) + (deltaEE * frac_ecr) ! Optionally add cosmicray flux, consistent with addition to ecr, via ! delta fcr = -K grad (delta ecr) ! Currently only set up for the 'gaussian3' profile in ecr. ! Still experimental/in testing if (lcosmicrayflux .and. lSN_fcr) then if (thermal_profile=='gaussian3') then do i=1,3 ! Currently hardwire factor of 0.05 as isotropic ecr diffusivity ! (Should instead use Kpara and Kperp properly.) f(l1:l2,m,n,ifcr+i-1) = f(l1:l2,m,n,ifcr+i-1) & + deltaEE*frac_ecr & * 0.05*(6.*(sqrt(dr2_SN)**5)/width_energy**6) & * outward_normal_SN(1:nx,i) enddo else call fatal_error("interstellar.explode_SN", & "fcr insertion only set up for thermal_profile='gaussian3'") endif endif endif ! ! Save changes to f-array. ! if (ldensity_nolog) then f(l1:l2,m,n,irho)=exp(lnrho) else f(l1:l2,m,n,ilnrho)=lnrho endif if (lSN_eth) then call eosperturb & (f,nx,ee=real((ee_old*rho_old+deltaEE*frac_eth)/exp(lnrho))) endif if (ldensity_nolog) then call eoscalc(irho_ss,f(l1:l2,m,n,irho),f(l1:l2,m,n,iss),& yH=yH,lnTT=lnTT) else call eoscalc(ilnrho_ss,f(l1:l2,m,n,ilnrho),f(l1:l2,m,n,iss),& yH=yH,lnTT=lnTT) endif lnTT=log(TT) if (lentropy.and.ilnTT/=0) f(l1:l2,m,n,ilnTT)=lnTT if (iyH/=0) f(l1:l2,m,n,iyH)=yH ! enddo enddo ! call get_properties(f,SNR,rhom_new,ekintot_new) if (lroot) print*,"TOTAL KINETIC ENERGY CHANGE:",ekintot_new-ekintot ! ! Sum and share diagnostics etc. amongst processors. ! dmpi2_tmp=(/ SNR%MM, SNR%EE /) call mpireduce_sum(dmpi2_tmp,dmpi2,2) call mpibcast_real(dmpi2,2) SNR%MM=dmpi2(1)*dv SNR%EE=dmpi2(2)*dv+ekintot_new-ekintot !include added kinetic energy ! if (lroot.and.ip<20) print*, & 'explode_SN: SNR%MM=',SNR%MM if (SNR%SN_type==1) then t_interval_SN=t_interval_SNI else t_interval_SN=t_interval_SNII endif ! if (lroot) then open(1,file=trim(datadir)//'/sn_series.dat',position='append') print*, 'explode_SN: step, time = ', it,t print*, 'explode_SN: dv = ', dv print*, 'explode_SN: SN type = ', SNR%SN_type print*, 'explode_SN: proc, l, m, n = ', SNR%iproc, SNR%l,SNR%m,SNR%n print*, 'explode_SN: x, y, z = ', SNR%x,SNR%y,SNR%z print*, 'explode_SN:remnant radius = ', SNR%radius print*, 'explode_SN: rho, TT = ', SNR%site%rho,SNR%site%TT print*, 'explode_SN: maximum TT = ', maxTT print*, 'explode_SN: Mean density = ', SNR%rhom print*, 'explode_SN: Total energy = ', SNR%EE print*, 'explode_SN: Added mass = ', SNR%MM print*, 'explode_SN: Ambient mass = ', site_mass print*, 'explode_SN: Sedov time = ', SNR%t_sedov print*, 'explode_SN: Shell velocity = ', uu_sedov write(1,'(i10,E13.5,5i6,11E13.5)') & it, t, SNR%SN_type, SNR%iproc, SNR%l, SNR%m, SNR%n, & SNR%x, SNR%y, SNR%z, SNR%site%rho, SNR%site%TT, SNR%EE, & SNR%t_sedov, SNR%radius, site_mass, maxTT, t_interval_SN close(1) endif ! if (present(preSN)) then do i=2,npreSN preSN(:,i-1)= preSN(:,i) enddo preSN(1,npreSN)= SNR%l preSN(2,npreSN)= SNR%m preSN(3,npreSN)= SNR%n preSN(4,npreSN)= SNR%iproc endif ! if (lSNR_damping) then SNR%state=SNstate_damping SNR%t_damping=SNR_damping_time-SNR%t_sedov+t else SNR%state=SNstate_finished endif ! if (present(ierr)) then ierr=iEXPLOSION_OK endif ! endsubroutine explode_SN !***************************************************************************** subroutine calc_snr_damping_factor(f) ! use Sub, only: dot,dot2 use Mpicomm ! real, intent(inout), dimension(mx,my,mz,mfarray) :: f real, dimension(nx,3) :: r_vec, uu real, dimension(nx) :: uur real, dimension(nx) :: r2, lnrho real :: radius2, fac, fac2 integer :: i, iSNR ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state/=SNstate_damping) cycle ! if ((t-SNRs(iSNR)%t_damping)/SNR_damping_rate >= 2.5) then if (lroot) print*,"No more damping!!" SNRs(iSNR)%state=SNstate_finished return endif ! fac=0. radius2=SNRs(iSNR)%radius**2 do n=n1,n2 do m=m1,m2 call proximity_SN(SNRs(iSNR)) uu = f(l1:l2,m,n,iux:iuz) if (ldensity_nolog) then lnrho = log(f(l1:l2,m,n,irho)) else lnrho = f(l1:l2,m,n,ilnrho) endif r_vec=0. r_vec(:,1) = x(l1:l2) - SNRs(iSNR)%x r_vec(:,2) = y(m) - SNRs(iSNR)%y r_vec(:,3) = z(n) - SNRs(iSNR)%z call dot2(r_vec,r2) call dot(r_vec,uu,uur) where (uur > 0.) uur=-uur/r2 endwhere where (r2 > radius2) uur=0. endwhere fac=max(fac,maxval(uur)) enddo enddo ! call mpiallreduce_max(fac,fac2) ! SNRs(iSNR)%damping_factor = fac2 * 2E-4 & * (1. - tanh((t-SNRs(iSNR)%t_damping)/SNR_damping_rate)) enddo ! endsubroutine calc_snr_damping_factor !***************************************************************************** subroutine calc_snr_unshock(penc) ! real, dimension(mx), intent(inout) :: penc real, dimension(mx) :: dr2_SN_mx integer :: i, iSNR ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state/=SNstate_damping) cycle call proximity_SN_mx(SNRs(iSNR),dr2_SN_mx) penc = penc*(1.-exp(-(dr2_SN_mx/SNRs(iSNR)%radius**2))) enddo ! endsubroutine calc_snr_unshock !***************************************************************************** subroutine calc_snr_damping(p) ! use Sub, only: multsv, multsv_add, dot ! type (pencil_case) :: p real, dimension(nx) :: profile real, dimension(nx,3) :: tmp, tmp2, gprofile integer :: i, iSNR ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state/=SNstate_damping) cycle if (lfirst) SNRs(iSNR)%energy_loss=0. ! call proximity_SN(SNRs(iSNR)) profile=SNRs(iSNR)%damping_factor*exp(-(dr2_SN/SNRs(iSNR)%radius**2)) gprofile=-SNRs(iSNR)%damping_factor* & spread(dr2_SN**2/(SNRs(iSNR)%radius**(2.5))* & exp(-(dr2_SN/SNRs(iSNR)%radius**2)),2,3) gprofile(:,1)= gprofile(:,1) * x(l1:l2) gprofile(:,2)= gprofile(:,2) * y(m) gprofile(:,3)= gprofile(:,3) * z(n) if (ldensity) then call multsv(p%divu,p%glnrho,tmp2) tmp=tmp2 + p%graddivu call multsv(profile,tmp,tmp2) call multsv_add(tmp2,p%divu,gprofile,tmp) if (lfirst.and.ldt) p%diffus_total=p%diffus_total+sum(profile) SNRs(iSNR)%energy_loss=SNRs(iSNR)%energy_loss-sum(profile*p%divu**2) p%fvisc=p%fvisc+tmp endif ! ! Have to divide by dxmin**2 to compensate for the * dxmin**2 in the shock ! code. This was commented out, but needs checking as I do not currently use ! damping so have generally left all damping references unchanged except for ! format. Fred ! ! penc=max(penc,SNR_damping* & ! exp(-(dr2_SN_mx/SNRs(iSNR)%radius**2))/dxmin**2) ! enddo ! endsubroutine calc_snr_damping !***************************************************************************** subroutine calc_snr_damp_int(int_dt) ! real :: int_dt integer :: i,iSNR ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state/=SNstate_damping) cycle SNRs(iSNR)%heat_energy = SNRs(iSNR)%heat_energy & + int_dt*SNRs(iSNR)%energy_loss*dv*int_dt enddo ! endsubroutine calc_snr_damp_int !*********************************************************************** subroutine calc_snr_damping_add_heat(f) ! use Mpicomm use EquationOfState, only: eoscalc, eosperturb ! real, intent(inout), dimension(mx,my,mz,mfarray) :: f real, dimension(nx) :: profile, ee_old, rho, lnrho real :: factor, fmpi, fmpi_tmp integer :: i, iSNR ! do i=1,nSNR iSNR=SNR_index(i) if (SNRs(iSNR)%state/=SNstate_damping) cycle ! fmpi_tmp=SNRs(iSNR)%heat_energy call mpireduce_sum(fmpi_tmp,fmpi) call mpibcast_real(fmpi) SNRs(iSNR)%heat_energy=fmpi ! factor = SNRs(iSNR)%heat_energy & / (cnorm_gaussian_SN(dimensionality)*SNRs(iSNR)%radius**dimensionality) do n=n1,n2 do m=m1,m2 call proximity_SN(SNRs(iSNR)) call eoscalc(f,nx,ee=ee_old) if (ldensity_nolog) then lnrho=log(f(l1:l2,m,n,irho)) rho=exp(lnrho) else lnrho=f(l1:l2,m,n,ilnrho) rho=exp(lnrho) endif profile = factor*exp(-(dr2_SN/SNRs(iSNR)%radius**2)) call eosperturb(f,nx,ee=real((ee_old*rho+profile)/exp(lnrho))) enddo enddo ! SNRs(iSNR)%heat_energy=0. ! if (t >= SNRs(iSNR)%t_damping) then SNRs(iSNR)%state=SNstate_finished endif enddo ! endsubroutine calc_snr_damping_add_heat !***************************************************************************** subroutine get_properties(f,remnant,rhom,ekintot) ! ! Calculate integral of mass cavity profile and total kinetic energy. ! ! 22-may-03/tony: coded ! 10-oct-09/axel: return zero density if the volume is zero ! use Sub use Mpicomm ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant) :: remnant real :: radius2 real :: rhom, ekintot real, dimension(nx) :: rho, u2 real, dimension(nx,3) :: uu integer, dimension(nx) :: mask real, dimension(3) :: tmp,tmp2 logical :: precise_sqrt=.true. ! ! Obtain distance to SN and sum all points inside SNR radius and ! divide by number of points. ! radius2 = (remnant%radius)**2 tmp=0.0 do n=n1,n2 do m=m1,m2 call proximity_SN(remnant) if (ldensity_nolog) then rho=f(l1:l2,m,n,irho) else rho=exp(f(l1:l2,m,n,ilnrho)) endif ! ! Necessary to compute ekintot in double precision here as very low rho ! can multiply very low u2 to make NaN. fred. ! uu=f(l1:l2,m,n,iuu:iuu+2) where (abs(f(l1:l2,m,n,iuu:iuu+2)) radius2) rho(1:nx)=0. mask(1:nx)=0 endwhere tmp(1)=tmp(1)+sum(rho) tmp(2)=tmp(2)+sum(mask) enddo enddo ! ! Calculate mean density inside the remnant and return error if the volume is ! zero. ! call mpireduce_sum(tmp,tmp2,3) call mpibcast_real(tmp2,3) ekintot=0.5*tmp2(3)*dv if (abs(tmp2(2)) < tini) then write(0,*) 'tmp2 = ', tmp2 call fatal_error("interstellar.get_properties","Dividing by zero?") else rhom=tmp2(1)/tmp2(2) endif ! endsubroutine get_properties !***************************************************************************** subroutine get_props_check(f,remnant,rhom,ekintot,cvelocity_SN,cmass_SN) ! ! Calculate integral of mass cavity profile and total kinetic energy. ! ! 22-may-03/tony: coded ! 10-oct-09/axel: return zero density if the volume is zero ! use Sub use Mpicomm ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant) :: remnant real, intent(in) :: cvelocity_SN, cmass_SN real :: radius2 real :: rhom, ekintot real :: width_mass, width_velocity real, dimension(nx) :: rho, u2, deltarho real, dimension(nx,3) :: uu real, dimension(nx,3) :: deltauu integer, dimension(nx) :: mask real, dimension(3) :: tmp,tmp2 logical :: precise_sqrt=.true. ! ! Obtain distance to SN and sum all points inside SNR radius and ! divide by number of points. ! width_mass = remnant%radius*mass_width_ratio width_velocity = remnant%radius*velocity_width_ratio radius2 = (remnant%radius)**2 tmp=0.0 do n=n1,n2 do m=m1,m2 call proximity_SN(remnant) if (ldensity_nolog) then rho=f(l1:l2,m,n,irho) else rho=exp(f(l1:l2,m,n,ilnrho)) endif if (lSN_mass) then call injectmass_SN(deltarho,width_mass,cmass_SN,remnant%MM) rho=rho+deltarho endif ! ! Necessary to compute ekintot in double precision here as very low rho ! can multiply very low u2 to make NaN. fred. ! uu=f(l1:l2,m,n,iuu:iuu+2) if (lSN_velocity) then call injectvelocity_SN(deltauu,width_velocity,cvelocity_SN) uu=uu+deltauu endif where (abs(f(l1:l2,m,n,iuu:iuu+2)) radius2) rho(1:nx)=0. mask(1:nx)=0 endwhere tmp(1)=tmp(1)+sum(rho) tmp(2)=tmp(2)+sum(mask) enddo enddo ! ! Calculate mean density inside the remnant and return zero if the volume is ! zero. ! tmp2=tmp call mpireduce_sum(tmp,tmp2,3) call mpibcast_real(tmp2,3) ekintot=0.5*tmp2(3)*dv if (abs(tmp2(2)) < tini) then write(0,*) 'tmp2 = ', tmp2 call fatal_error("interstellar.get_props_check","Dividing by zero?") else rhom=tmp2(1)/tmp2(2) endif ! endsubroutine get_props_check !***************************************************************************** subroutine get_lowest_rho(f,SNR,radius,rho_lowest) ! ! Calculate integral of mass cavity profile. ! ! 22-may-03/tony: coded ! use Mpicomm ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant), intent(in) :: SNR real, intent(in) :: radius real, intent(out) :: rho_lowest real :: tmp real :: radius2 real, dimension(nx) :: rho ! ! Find lowest rho value in the surronding cavity. ! rho_lowest=1E10 radius2 = radius**2 do n=n1,n2 do m=m1,m2 call proximity_SN(SNR) if (ldensity_nolog) then rho=f(l1:l2,m,n,irho) else rho=exp(f(l1:l2,m,n,ilnrho)) endif where (dr2_SN(1:nx) > radius2) rho=1E10 rho_lowest=min(rho_lowest,minval(rho(1:nx))) enddo enddo ! tmp=-exp(rho_lowest) call mpireduce_max(tmp,rho_lowest) call mpibcast_real(rho_lowest) ! endsubroutine get_lowest_rho !***************************************************************************** subroutine proximity_SN(SNR) ! ! Calculate pencil of distance to SN explosion site. ! ! 20-may-03/tony: extracted from explode_SN code written by grs ! 22-may-03/tony: pencil formulation ! type (SNRemnant), intent(in) :: SNR ! real,dimension(nx) :: dx_SN, dr_SN real :: dy_SN real :: dz_SN ! ! Obtain distance to SN ! dx_SN=x(l1:l2)-SNR%x if (lperi(1)) then where (dx_SN > Lx/2.) dx_SN=dx_SN-Lx where (dx_SN < -Lx/2.) dx_SN=dx_SN+Lx endif ! dy_SN=y(m)-SNR%y if (lperi(2)) then if (dy_SN > Ly/2.) dy_SN=dy_SN-Ly if (dy_SN < -Ly/2.) dy_SN=dy_SN+Ly endif ! dz_SN=z(n)-SNR%z if (lperi(3)) then if (dz_SN > Lz/2.) dz_SN=dz_SN-Lz if (dz_SN < -Lz/2.) dz_SN=dz_SN+Lz endif ! dr2_SN=dx_SN**2 + dy_SN**2 + dz_SN**2 ! if (lSN_velocity) then dr_SN=sqrt(dr2_SN) dr_SN=max(dr_SN(1:nx),tiny(0.0)) ! ! Avoid dr_SN = 0 above to avoid div by zero below. ! outward_normal_SN(:,1)=dx_SN/dr_SN where (dr2_SN(1:nx) == 0.) outward_normal_SN(:,1)=0.0 outward_normal_SN(:,2)=dy_SN/dr_SN where (dr2_SN(1:nx) == 0.) outward_normal_SN(:,2)=0.0 outward_normal_SN(:,3)=dz_SN/dr_SN where (dr2_SN(1:nx) == 0.) outward_normal_SN(:,3)=0.0 endif ! endsubroutine proximity_SN !***************************************************************************** subroutine proximity_SN_mx(SNR,dr2_SN_mx) ! ! Calculate pencil of distance to SN explosion site. ! ! 20-may-03/tony: extracted from explode_SN code written by grs ! 22-may-03/tony: pencil formulation ! type (SNRemnant), intent(in) :: SNR real,dimension(mx), intent(out) :: dr2_SN_mx real,dimension(mx) :: dx_SN real :: dy_SN real :: dz_SN ! ! Obtain distance to SN. ! dx_SN=x-SNR%x if (lperi(1)) then where (dx_SN > Lx/2.) dx_SN=dx_SN-Lx where (dx_SN < -Lx/2.) dx_SN=dx_SN+Lx endif ! dy_SN=y(m)-SNR%y if (lperi(2)) then if (dy_SN > Ly/2.) dy_SN=dy_SN-Ly if (dy_SN < -Ly/2.) dy_SN=dy_SN+Ly endif ! dz_SN=z(n)-SNR%z if (lperi(3)) then if (dz_SN > Lz/2.) dz_SN=dz_SN-Lz if (dz_SN < -Lz/2.) dz_SN=dz_SN+Lz endif ! dr2_SN_mx=dx_SN**2 + dy_SN**2 + dz_SN**2 ! endsubroutine proximity_SN_mx !***************************************************************************** subroutine calc_cavity_mass_lnrho(f,SNR,width,depth,mass_removed) ! ! Calculate integral of mass cavity profile. ! ! 22-may-03/tony: coded ! use Mpicomm, only: mpibcast_real, mpireduce_sum ! real, intent(in), dimension(mx,my,mz,mfarray) :: f type (SNRemnant), intent(in) :: SNR real, intent(in) :: width, depth real, intent(out) :: mass_removed real, dimension(nx) :: lnrho, lnrho_old real, dimension(nx) :: rho real :: dmpi1, dmpi1_tmp real, dimension(nx) :: profile_cavity ! ! Obtain distance to SN ! mass_removed=0. do n=n1,n2 do m=m1,m2 call proximity_SN(SNR) ! if (ldensity_nolog) then lnrho_old=log(f(l1:l2,m,n,irho)) else lnrho_old=f(l1:l2,m,n,ilnrho) endif if (cavity_profile=="gaussian3log") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**3)) lnrho=lnrho_old - profile_cavity mass_removed=mass_removed+sum(exp(lnrho_old)-exp(lnrho)) elseif (cavity_profile=="gaussian3") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**3)) lnrho=lnrho_old - profile_cavity mass_removed=mass_removed+sum(exp(lnrho_old)-exp(lnrho)) elseif (cavity_profile=="gaussian2") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**2)) lnrho=lnrho_old - profile_cavity mass_removed=mass_removed+sum(exp(lnrho_old)-exp(lnrho)) elseif (cavity_profile=="gaussian") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2))) lnrho=lnrho_old - profile_cavity mass_removed=mass_removed+sum(exp(lnrho_old)-exp(lnrho)) elseif (cavity_profile=="tanh") then profile_cavity=(1.-tanh( (width-sqrt(dr2_SN(1:nx)) ) *sigma_SN1 ))*0.5 rho=exp(lnrho_old)*profile_cavity mass_removed=mass_removed+sum(exp(lnrho_old)-rho) endif ! enddo enddo dmpi1_tmp=mass_removed call mpireduce_sum(dmpi1_tmp,dmpi1) call mpibcast_real(dmpi1) mass_removed=dmpi1*dv ! endsubroutine calc_cavity_mass_lnrho !*********************************************************************** subroutine make_cavity_rho(deltarho,width,depth, & cnorm_dim,MMtot_SN) ! real, intent(in) :: width, depth, cnorm_dim real, intent(inout) :: MMtot_SN real, intent(out), dimension(nx) :: deltarho ! real, dimension(nx) :: profile_shell_outer,profile_shell_inner real :: width_shell_outer, width_shell_inner, c_shell ! width_shell_outer=outer_shell_proportion*width width_shell_inner=inner_shell_proportion*width ! ! deltarho(1:nx) = -depth*exp(-(dr2_SN(1:nx)/width**2)**3) ! c_shell=-depth*cnorm_dim/((1./width_shell_outer**dimensionality)- & (1./width_shell_inner**dimensionality)) ! ! Add missing mass back into shell. ! profile_shell_outer(1:nx)= & exp(-(dr2_SN(1:nx)/width_shell_outer**2)**3)/ & cnorm_dim/width_shell_outer**dimensionality profile_shell_inner(1:nx)= & exp(-(dr2_SN(1:nx)/width_shell_inner**2)**3)/ & cnorm_dim/width_shell_inner**dimensionality deltarho(1:nx)=c_shell* & (profile_shell_outer(1:nx) - profile_shell_inner(1:nx)) MMtot_SN=MMtot_SN + sum(deltarho(1:nx)) ! endsubroutine make_cavity_rho !***************************************************************************** subroutine make_cavity_lnrho(lnrho,width,depth,mass_shell, & cnorm_dim,MMtot_SN) ! real, intent(in) :: width, depth, mass_shell, cnorm_dim real, intent(inout) :: MMtot_SN real, intent(inout), dimension(nx) :: lnrho ! real, dimension(nx) :: profile_shell_outer,profile_cavity real, dimension(nx) :: profile_shell_inner real :: width_shell_outer, width_shell_inner, c_shell real :: mass_before, mass_after ! width_shell_outer=outer_shell_proportion*width width_shell_inner=inner_shell_proportion*width ! c_shell = mass_shell/(cnorm_dim*(width_shell_outer**dimensionality - & width_shell_inner**dimensionality)) ! profile_shell_outer(1:nx)=exp(-(dr2_SN(1:nx)/width_shell_outer**2)**3) profile_shell_inner(1:nx)=exp(-(dr2_SN(1:nx)/width_shell_inner**2)**3) ! mass_before=sum(exp(lnrho(1:nx))) if (cavity_profile=="gaussian3log") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**3)) lnrho = lnrho(1:nx) - profile_cavity lnrho = log(exp(lnrho(1:nx))+c_shell* & (profile_shell_outer(1:nx)-profile_shell_inner(1:nx))) elseif (cavity_profile=="gaussian3") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**3)) lnrho = lnrho(1:nx)-profile_cavity lnrho = log(exp(lnrho(1:nx))+c_shell* & (profile_shell_outer(1:nx)-profile_shell_inner(1:nx))) elseif (cavity_profile=="gaussian2") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2)**2)) lnrho = lnrho(1:nx)-profile_cavity lnrho = log(exp(lnrho(1:nx))+c_shell* & (profile_shell_outer(1:nx)-profile_shell_inner(1:nx))) elseif (cavity_profile=="gaussian") then profile_cavity=(depth*exp(-(dr2_SN(1:nx)/width**2))) lnrho = lnrho(1:nx)-profile_cavity lnrho = log(exp(lnrho(1:nx))+c_shell* & (profile_shell_outer(1:nx)-profile_shell_inner(1:nx))) elseif (cavity_profile=="tanh") then profile_cavity=(1.-tanh((width-sqrt(dr2_SN(1:nx)))*sigma_SN1))*0.5 lnrho = log(exp(lnrho(1:nx))*profile_cavity+c_shell* & (profile_shell_outer(1:nx)-profile_shell_inner(1:nx))+ & depth*(1.-tanh((sqrt(dr2_SN(1:nx))-width)*sigma_SN1))*0.5) endif mass_after=sum(exp(lnrho(1:nx))) MMtot_SN=MMtot_SN + (mass_after-mass_before) ! endsubroutine make_cavity_lnrho !***************************************************************************** subroutine injectenergy_SN(deltaEE,width,c_SN,EEtot_SN) ! real, intent(in) :: width,c_SN real, intent(inout) :: EEtot_SN real, intent(out), dimension(nx) :: deltaEE ! real, dimension(nx) :: profile_SN ! ! Whether mass is moved or not, inject energy. ! if (thermal_profile=="gaussian3") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**3) elseif (thermal_profile=="gaussian2") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**2) elseif (thermal_profile=="gaussian") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)) elseif (thermal_profile=="quadratic") then profile_SN=max(1.0-(dr2_SN(1:nx)/width**2),0.0) elseif (thermal_profile=="quadratictanh") then profile_SN=max(1.0-(dr2_SN(1:nx)/width**2),0.0)* & 0.5*(1.-tanh((sqrt(dr2_SN)-width)*sigma_SN1)) elseif (thermal_profile=="quartictanh") then profile_SN=max(1.0-(dr2_SN(1:nx)/width**2)**2,0.0)* & 0.5*(1.-tanh((sqrt(dr2_SN)-width)*sigma_SN1)) elseif (thermal_profile=="tanh") then profile_SN=(1.-tanh((sqrt(dr2_SN(1:nx))-width)*sigma_SN1))*0.5 endif ! deltaEE(1:nx)=c_SN*profile_SN(1:nx) ! spatial energy density EEtot_SN=EEtot_SN+sum(deltaEE(1:nx)) ! endsubroutine injectenergy_SN !***************************************************************************** subroutine injectmass_SN(deltarho,width,cmass_SN,MMtot_SN) ! real, intent(in) :: width,cmass_SN real, intent(inout) :: MMtot_SN real, intent(out), dimension(nx) :: deltarho ! real, dimension(nx) :: profile_SN ! ! Inject mass. ! if (mass_profile=="gaussian3") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**3) elseif (mass_profile=="gaussian2") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**2) elseif (mass_profile=="gaussian") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)) elseif (mass_profile=="quadratic") then profile_SN=max(1.0-(dr2_SN(1:nx)/width**2),0.0) elseif (mass_profile=="tanh") then ! ! This is normally handled in the mass movement section ! profile_SN=(1.-tanh((sqrt(dr2_SN(1:nx))-width)*sigma_SN1))*0.5 endif ! deltarho(1:nx)=cmass_SN*profile_SN(1:nx) ! spatial mass density MMtot_SN=MMtot_SN+sum(deltarho(1:nx)) ! endsubroutine injectmass_SN !*********************************************************************** subroutine injectvelocity_SN(deltauu,width,cvelocity_SN) ! real, intent(in) :: width,cvelocity_SN real, intent(out), dimension(nx,3) :: deltauu ! real, dimension(nx) :: profile_SN ! integer :: j ! ! Calculate deltauu. ! if (velocity_profile=="quintictanh") then profile_SN=((sqrt(dr2_SN)/width)**5)*0.5* & (1.-tanh((sqrt(dr2_SN)-(1.1*width))/(0.08*width))) ! elseif (velocity_profile=="lineartanh") then profile_SN=max(sqrt(dr2_SN)/width,1.0)*0.5* & (1.-tanh((sqrt(dr2_SN)-width)*sigma_SN1-2.)) ! elseif (velocity_profile=="quadratictanh") then profile_SN=min((dr2_SN/(width**2)),0.5* & (1.-tanh((sqrt(dr2_SN)-width)*sigma_SN1-2.))) ! elseif (velocity_profile=="gaussian") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)) ! elseif (velocity_profile=="gaussian2") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**2) ! elseif (velocity_profile=="gaussian3") then profile_SN=exp(-(dr2_SN(1:nx)/width**2)**3) ! elseif (velocity_profile=="r8thgaussian3") then profile_SN=(dr2_SN(1:nx))**0.0625*exp(-(dr2_SN(1:nx)/width**2)**3) ! elseif (velocity_profile=="r8thgaussian") then profile_SN=(dr2_SN(1:nx))**0.0625*exp(-(dr2_SN(1:nx)/width**2)) ! elseif (velocity_profile=="r16thgaussian3") then profile_SN=(dr2_SN(1:nx))**0.03125*exp(-(dr2_SN(1:nx)/width**2)**3) ! elseif (velocity_profile=="r16thgaussian") then profile_SN=(dr2_SN(1:nx))**0.03125*exp(-(dr2_SN(1:nx)/width**2)) ! elseif (velocity_profile=="cubictanh") then profile_SN=(sqrt(dr2_SN/width)**3)* & (1.-tanh((sqrt(dr2_SN)-(1.1*width))*sigma_SN1)) ! elseif (velocity_profile=="gaussian3der") then profile_SN=(((sqrt(dr2_SN)**5)/width**6*(1./35.))* & exp(-(dr2_SN(1:nx)/width**2)**3)) ! elseif (velocity_profile=="quadratic") then profile_SN=dr2_SN(1:nx)/width**2 where (dr2_SN>(width**2)) profile_SN=0. endif ! do j=1,3 deltauu(1:nx,j)=cvelocity_SN*profile_SN(1:nx)* & outward_normal_SN(1:nx,j) ! spatial mass density enddo ! endsubroutine injectvelocity_SN !***************************************************************************** function get_free_SNR() ! integer :: get_free_SNR integer :: i,iSNR ! if (nSNR>=mSNR) then call fatal_error("get_free_SNR", & "Run out of SNR slots... Increase mSNR.") endif ! iSNR=-1 do i=1,mSNR if (SNRs(i)%state==SNstate_invalid) then iSNR=i exit endif enddo ! if (iSNR<0) then call fatal_error("get_free_SNR", & "Could not find an empty SNR slot. Slots were not properly freed.") endif ! nSNR=nSNR+1 SNRs(iSNR)%state=SNstate_waiting SNR_index(nSNR)=iSNR get_free_SNR=iSNR ! endfunction get_free_SNR !***************************************************************************** subroutine free_SNR(iSNR) ! integer :: i,iSNR ! if (SNRs(iSNR)%state==SNstate_invalid) then if (lroot) print*,"Tried to free an already invalid SNR" return endif ! nSNR=nSNR-1 SNRs(iSNR)%state=SNstate_invalid ! do i=iSNR,nSNR SNR_index(i)=SNR_index(i+1) enddo ! endsubroutine free_SNR !***************************************************************************** subroutine tidy_SNRs ! integer :: i ! do i=1,mSNR if (SNRs(i)%state==SNstate_finished) call free_SNR(i) enddo ! endsubroutine tidy_SNRs !***************************************************************************** subroutine addmassflux(f) ! ! This routine calculates the mass flux through the vertical boundary. ! As no inflow boundary condition precludes galactic fountain this adds ! the mass flux proportionately throughout the volume to substitute mass ! which would otherwise be replaced over time by the galactic fountain. ! ! 23-Nov-11/fred: coded ! real, intent(inout), dimension(mx,my,mz,mfarray) :: f ! real :: prec_factor=1.0E-7 real :: add_ratio integer :: l,m,n ! ! Skip this subroutine if not selected eg before turbulent pressure settles ! if (.not. ladd_massflux) return ! ! Only add boundary mass at intervals to reduce MPI operations and to ! ensure flux replacement are large enough to reduce losses at the limit ! of machine accuracy. At same frequency as SNII. ! if (t >= t_next_mass) then ! ! Determine multiplier required to restore mass to level before boundary ! losses. addrate (default=1.0) can be increased to raise mass levels if ! required. ! add_ratio=1.0+prec_factor*addrate ! ! Add mass proportionally to the existing density throughout the ! volume to replace that lost through boundary. ! add_ratio needs to be large enough for single precision to record small ! changes, so accumulate small mass losses in boldmass until large enough. ! if (ldensity_nolog) then f(l1:l2,m1:m2,n1:n2,irho)= & dble(f(l1:l2,m1:m2,n1:n2,irho))*add_ratio else f(l1:l2,m1:m2,n1:n2,ilnrho)= & dble(f(l1:l2,m1:m2,n1:n2,ilnrho))+log(add_ratio) endif t_next_mass=t_next_SNII endif ! endsubroutine addmassflux !***************************************************************************** ! subroutine addmassflux(f) !! !! This routine calculates the mass flux through the vertical boundary. !! As no inflow boundary condition precludes galactic fountain this adds !! the mass flux proportionately throughout the volume to substitute mass !! which would otherwise be replaced over time by the galactic fountain. !! !! 12-Jul-10/fred: coded !! This older routine is numerically expensive and above was adopted using a !! simple factor to keep the mean density steady following hydrodynamic !! steady turbulence ! ! use Mpicomm, only: mpireduce_sum, mpibcast_real, & ! mpibcast_real, mpireduce_max !! ! real, intent(inout), dimension(mx,my,mz,mfarray) :: f !! ! real :: prec_factor=1.0E-6 ! real :: sum_tmp, rmpi ! real :: bflux, bmass, add_ratio, rhosum ! real :: bfmpi, sum_1tmp, nmpi, sum_3tmp ! real :: newmass, oldmass ! integer :: l,m,n !! !! Skip this subroutine if not selected eg before turbulent pressure settles !! ! if (.not. ladd_massflux) return !! !! Calculate the total flux through the vertical boundaries to determine !! mass loss to the system. Sum the total mass in the domain. !! ! if (ldensity_nolog) then ! rhosum = sum(dble(f(l1:l2,m1:m2,n1:n2,irho))) ! else ! rhosum = sum(exp(dble(f(l1:l2,m1:m2,n1:n2,ilnrho)))) ! endif ! sum_tmp=rhosum ! call mpireduce_sum(sum_tmp,rmpi) ! call mpibcast_real(rmpi) ! rhosum=rmpi ! oldmass=rhosum*dv !! !! Calculate mass loss through the vertical boundaries rho*u_z !! ! bflux=0.0 ! do n=n1,n2 ! if (z(n) == xyz0(3)) then ! do l=l1,l2 ! do m=m1,m2 ! if (ldensity_nolog) then ! bflux=bflux-dble(f(l,m,n,irho))*dble(f(l,m,n,iuz)) ! else ! bflux=bflux-exp(dble(f(l,m,n,ilnrho)))*dble(f(l,m,n,iuz)) ! endif ! enddo ! enddo ! endif ! if (z(n) == xyz1(3)) then ! do l=l1,l2 ! do m=m1,m2 ! if (ldensity_nolog) then ! bflux=bflux+dble(f(l,m,n,irho))*dble(f(l,m,n,iuz)) ! else ! bflux=bflux+exp(dble(f(l,m,n,ilnrho)))*dble(f(l,m,n,iuz)) ! endif ! enddo ! enddo ! endif ! enddo ! sum_1tmp=bflux !! ! if (ip<45.and.bflux /=0.0) print*,'addmassflux: bflux on iproc =', & ! bflux, iproc !! !! Sum over all processors and communicate total to all. !! ! call mpireduce_sum(sum_1tmp,bfmpi) ! call mpibcast_real(bfmpi) ! bflux=bfmpi ! if (lroot.and.ip<45) print*,'addmassflux: bflux after mpi sum =', bflux ! if (bflux>0.0) then !! !! Multiply mass flux by area element and timestep to determine lost mass. !! Add unused flux mass from previous timesteps. !! ! bmass=bflux*dt*dx*dy+boldmass !! !! Determine multiplier required to restore mass to level before boundary !! losses. addrate (default=1.0) can be increased to raise mass levels if !! required. !! ! add_ratio=(bmass*addrate+oldmass)/oldmass !! ! if (lroot.and.ip<45) print*, & ! 'addmassflux: bmass, add_ratio, timestep =', & ! bmass, add_ratio, dt !! !! Add mass proportionally to the existing density throughout the !! volume to replace that lost through boundary. !! add_ratio needs to be large enough for single precision to record small !! changes, so accumulate small mass losses in boldmass until large enough. !! ! if (add_ratio>=prec_factor+1.0) then ! if (ldensity_nolog) then ! f(l1:l2,m1:m2,n1:n2,irho)= & ! dble(f(l1:l2,m1:m2,n1:n2,irho))*add_ratio ! else ! f(l1:l2,m1:m2,n1:n2,ilnrho)= & ! dble(f(l1:l2,m1:m2,n1:n2,ilnrho))+log(add_ratio) ! endif !! !! For debugging purposes newmass can be calculated and compared to !! bmass+oldmass, which should be equal. !! ! if (ldensity_nolog) then ! rhosum=sum(dble(f(l1:l2,m1:m2,n1:n2,irho))) ! else ! rhosum=sum(exp(dble(f(l1:l2,m1:m2,n1:n2,ilnrho)))) ! endif ! sum_3tmp=rhosum ! call mpireduce_sum(sum_3tmp,nmpi) ! call mpibcast_real(nmpi) ! rhosum=nmpi ! newmass=nmpi*dv ! if (lroot.and.ip<45) print*,'addmassflux: oldmass, newmass =', & ! oldmass, newmass ! if (lroot.and.ip<70) print*, & ! 'addmassflux: added mass vs mass flux=', & ! newmass-oldmass, bmass ! boldmass=0.0 ! else ! boldmass=bmass ! endif ! endif !! ! endsubroutine addmassflux !!***************************************************************************** endmodule Interstellar