! $Id$ ! ! This modules contains the routines for SNe-driven ISM simulations. ! Still in development. ! !***************** 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 RemnantFeature real :: x, y, z, t ! Time and location real :: MM, EE, CR ! Mass, energy & CR injected real :: rhom ! Local mean density at explosion time real :: radius ! Injection radius real :: t_sedov endtype ! type RemnantIndex integer :: l,m,n ! Grid position integer :: iproc,ipx,ipy,ipz integer :: SN_type integer :: state endtype ! type SNRemnant type (RemnantFeature) :: feat type (RemnantIndex) :: indx type (ExplosionSite) :: site endtype ! ! required for *put_persistant_interstellar, update integer value to match ! any changes to number of above types ! integer :: nSITE = 7, nFEAT = 10, nINDX = 9 ! ! 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_finished = 3 ! ! 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, mu ! ! 'Current' SN Explosion site parameters ! integer, parameter :: mSNR = 10 integer :: nSNR = 0 type (SNRemnant), dimension(mSNR) :: SNRs real, dimension(12) :: SNRsFEAT integer, dimension(9) :: SNRsINDX real, dimension(7) :: SNRsSITE integer, dimension(mSNR) :: SNR_index integer, parameter :: npreSN = 5 integer, dimension(4,npreSN) :: preSN ! ! 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_cluster=0.0, y_cluster=0.0, z_cluster=0.0, t_cluster=0.0 real :: t_interval_SNI=impossible, t_interval_SNII=impossible real :: zdisk !varying location of centre of mass of the disk ! ! 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, unit_Gamma ! ! SNe placement limitations (for code stability) ! Minimum resulting central temperature of a SN explosion. ! If this is not reached then consider moving mass to achieve this. ! 10-aug-10/fred: ! As per joung et al apj653 2005 min temp 1E6 to avoid excess radiative ! energy losses in early stages. ! real, parameter :: rho_SN_min_cgs=1E-28,rho_SN_max_cgs=5E-24 real, parameter :: TT_SN_min_cgs=1.E6, TT_SN_max_cgs=5E9 real :: rho_SN_min=impossible, rho_SN_max=impossible real :: TT_SN_min=impossible, TT_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 ! ! 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, kperp=0.05, kpara=0.025 ! ! 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 ! ! 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. ! ! Cooling & heating flags ! logical :: laverage_SN_heating = .false. logical :: lheating_UV = .true. ! ! Remnant location flags ! logical :: luniform_zdist_SNI = .false., lSNII_gaussian=.true. 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_time=impossible real :: SN_clustering_radius_cgs=2.5e20 ! cm (~80 pc) real :: SN_clustering_time_cgs=1.6e14 ! cm (~5 Myr) ! ! Adjust SNR%feat%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= 'gaussian3' character (len=labellen) :: mass_profile = 'gaussian3' ! ! Variables required for returning mass to disk given no inflow ! boundary condition used in addmassflux ! real :: addrate=1.0 real :: boldmass=0.0 logical :: ladd_massflux = .false. ! ! Gravity constansts - rquired for pre-2015 vertical heating profile ! real :: a_S, a_D real, parameter :: a_S_cgs=4.4e-9, a_D_cgs=1.7e-9 real :: z_S, z_D, H_z real, parameter :: z_S_cgs=6.172e20, z_D_cgs=3.086e21, H_z_cgs=9.258E20 ! ! 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, & frac_ecr, frac_eth, thermal_profile, velocity_profile, mass_profile, & luniform_zdist_SNI, inner_shell_proportion, outer_shell_proportion, & 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, kperp, kpara ! ! 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, & luniform_zdist_SNI, SNI_area_rate, SNII_area_rate, & inner_shell_proportion, outer_shell_proportion, & frac_ecr, frac_eth, thermal_profile,velocity_profile, mass_profile, & h_SNI, h_SNII, TT_SN_min, lSN_scale_rad, & mass_SN_progenitor, cloud_tau, cdt_tauc, cloud_rho, cloud_TT, & laverage_SN_heating, coolingfunction_scalefactor,& 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, & heatcool_shock_cutoff, heatcool_shock_cutoff_rate, ladd_massflux, & N_mass, addrate, T0hs, rho0ts, & lSNII_gaussian, rho_SN_max, lSN_mass_rate, lthermal_hse, lheatz_min, & p_OB, SN_clustering_time, SN_clustering_radius, lOB_cluster, kperp, kpara ! 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(:)%indx%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 ! f(:,:,:,icooling)=0.0 f(:,:,:,inetheat)=0.0 ! if (lroot) print*,'initialize_interstellar: t_next_SNI',t_next_SNI ! if (lroot.and.luniform_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 ! 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 (n Gamma(T)) unit_Lambda = unit_velocity**2 / unit_density / unit_time unit_Gamma = unit_velocity**3 / unit_length ! T0UV=T0UV_cgs / unit_temperature cUV=cUV_cgs * unit_temperature if (GammaUV==impossible) & GammaUV=GammaUV_cgs / unit_Gamma ! 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 if (TT_SN_max==impossible) TT_SN_max=TT_SN_max_cgs / unit_temperature if (TT_SN_min==impossible) TT_SN_min=TT_SN_min_cgs / unit_temperature 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 if (h_SNII==impossible) h_SNII=h_SNII_cgs / unit_length SNII_mass_rate=SNII_mass_rate_cgs*unit_time SNI_mass_rate=SNI_mass_rate_cgs*unit_time solar_mass=solar_mass_cgs / unit_mass if (lroot) & print*,'initialize_interstellar: solar_mass (code) =', solar_mass 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 if (ampl_SN==impossible) then ampl_SN=ampl_SN_cgs / unit_energy if (lcosmicray .and. lSN_ecr) ampl_SN=frac_eth*ampl_SN endif 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_rho==impossible) cloud_rho=cloud_rho_cgs / unit_density if (cloud_TT==impossible) cloud_TT=cloud_TT_cgs / unit_temperature 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 (SN_clustering_radius==impossible) & SN_clustering_radius=SN_clustering_radius_cgs / unit_length if (SN_clustering_time==impossible) & SN_clustering_time=SN_clustering_time_cgs / unit_time else call stop_it('initialize_interstellar: SI unit conversions not implemented') endif ! call getmu(f,mu) call select_cooling(cooling_select,lncoolT,lncoolH,coolB) ! if (lroot) print*,'initialize_interstellar: unit_Lambda',unit_Lambda if (lroot) print*,'initialize_interstellar: unit_Gamma',unit_Gamma ! heating_rate_code=heating_rate*real(unit_length/unit_velocity**3) ! if (heating_select == 'thermal-hs') then ! call thermal_hs(f,zrho) call heat_interstellar(f,heat_z) 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 ! ! 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 write (1,'(a,1pd26.16)') 'unit_Gamma=',unit_Gamma close (1) endif ! endsubroutine initialize_interstellar !***************************************************************************** subroutine select_cooling(cooling_select,lncoolT,lncoolH,coolB) ! ! Routine for selecting parameters for temperature dependent cooling ! Lambda. ! character (len=labellen), intent(IN) :: cooling_select real, dimension (:), intent(OUT) :: lncoolT, coolB double precision, dimension (:), intent(OUT) :: lncoolH ! ! 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)) ! endsubroutine select_cooling !***************************************************************************** subroutine input_persist_interstellar_id(id,done) ! ! Read in the stored time of the next SNI ! ! 13-Dec-2011/Bourdin.KIS: reworked ! 14-jul-2015/fred: removed obsolete Remnant persistant variable from current ! read and added new cluster variables. All now consistent with any io ! use IO, only: read_persist, lun_input, lcollective_IO ! integer :: id logical :: done ! integer :: i ! ! if (lcollective_IO) call fatal_error ('input_persist_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_cluster, y_cluster done = .true. case (id_record_ISM_TOGGLE_OLD) read (lun_input) lSNI, lSNII done = .true. case (id_record_ISM_BOLD_MASS) if (read_persist ('ISM_BOLD_MASS', boldmass)) return done = .true. case (id_record_ISM_SNRS_OLD) ! Forget any existing SNRs. SNRs(:)%indx%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_cluster)) return done = .true. case (id_record_ISM_Y_CLUSTER) if (read_persist ('ISM_Y_CLUSTER', y_cluster)) return done = .true. case (id_record_ISM_Z_CLUSTER) if (read_persist ('ISM_Z_CLUSTER', z_cluster)) return done = .true. case (id_record_ISM_T_CLUSTER) if (read_persist ('ISM_T_CLUSTER', t_cluster)) 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. endselect ! if (lroot) & print *, 'input_persist_interstellar: ', t_next_SNI, t_next_SNII ! endsubroutine input_persist_interstellar_id !***************************************************************************** subroutine input_persist_interstellar() ! ! Read in the stored time of the next SNI. ! ! 12-Oct-2019/PABourdin: coded ! use IO, only: read_persist ! logical :: error ! if (.not. l_persist_overwrite_lSNI) then call warning('input_persist_interstellar','lSNI from run.in overwritten. ' // & 'Set l_persist_overwrite_lSNI=T to update', 0) error = read_persist ('ISM_TOGGLE_SNI', lSNI) if (lroot .and. .not. error) print *, 'input_persist_interstellar: lSNI = ', lSNI endif ! if (.not. l_persist_overwrite_lSNII) then call warning('input_persist_interstellar','lSNII from run.in overwritten. ' // & 'Set l_persist_overwrite_lSNII=T to update', 0) error = read_persist ('ISM_TOGGLE_SNII', lSNII) if (lroot .and. .not. error) print *, 'input_persist_interstellar: lSNII = ', lSNII endif ! if (.not. l_persist_overwrite_tSNI) then call warning('input_persist_interstellar','t_next_SNI from run.in overwritten. ' // & 'Set l_persist_overwrite_tSNI=T to update', 0) error = read_persist ('ISM_T_NEXT_SNI', t_next_SNI) if (lroot .and. .not. error) print *, 'input_persist_interstellar: t_next_SNI = ', t_next_SNI endif ! if (.not. l_persist_overwrite_tSNII) then call warning('input_persist_interstellar','t_next_SNII from run.in overwritten. ' // & 'Set l_persist_overwrite_tSNII=T to update', 0) error = read_persist ('ISM_T_NEXT_SNII', t_next_SNII) if (lroot .and. .not. error) print *, 'input_persist_interstellar: t_next_SNII = ', t_next_SNII endif ! if (.not. l_persist_overwrite_xcluster) then call warning('input_persist_interstellar','x_cluster from run.in overwritten. ' // & 'Set l_persist_overwrite_xcluster=T to update', 0) error = read_persist ('ISM_X_CLUSTER', x_cluster) if (lOB_cluster .and. lroot .and. .not. error) print *, 'input_persist_interstellar: x_cluster = ', x_cluster endif ! if (.not. l_persist_overwrite_ycluster) then call warning('input_persist_interstellar','y_cluster from run.in overwritten. ' // & 'Set l_persist_overwrite_ycluster=T to update', 0) error = read_persist ('ISM_Y_CLUSTER', y_cluster) if (lOB_cluster .and. lroot .and. .not. error) print *, 'input_persist_interstellar: y_cluster = ', y_cluster endif ! if (.not. l_persist_overwrite_zcluster) then call warning('input_persist_interstellar','z_cluster from run.in overwritten. ' // & 'Set l_persist_overwrite_zcluster=T to update', 0) error = read_persist ('ISM_Z_CLUSTER', z_cluster) if (lOB_cluster .and. lroot .and. .not. error) print *, 'input_persist_interstellar: z_cluster = ', z_cluster endif ! if (.not. l_persist_overwrite_tcluster) then call warning('input_persist_interstellar','t_cluster from run.in overwritten. ' // & 'Set l_persist_overwrite_tcluster=T to update', 0) error = read_persist ('ISM_T_CLUSTER', t_cluster) if (lOB_cluster .and. lroot .and. .not. error) print *, 'input_persist_interstellar: t_cluster = ', t_cluster endif ! endsubroutine input_persist_interstellar !***************************************************************************** logical function output_persistent_interstellar() ! ! Writes out the time of the next SNI ! ! 13-Dec-2011/Bourdin.KIS: reworked ! 14-jul-2015/fred: removed obsolete Remnant persistant variable from current ! write and added new cluster variables. All now consistent with any io ! 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_cluster)) return if (write_persist ('ISM_Y_CLUSTER', id_record_ISM_Y_CLUSTER, y_cluster)) return if (write_persist ('ISM_Z_CLUSTER', id_record_ISM_Z_CLUSTER, z_cluster)) return if (write_persist ('ISM_T_CLUSTER', id_record_ISM_T_CLUSTER, t_cluster)) 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 ! 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 idiag_taucmin=0 idiag_Hmax=0 idiag_Lamm=0 idiag_nrhom=0 idiag_rhoLm=0 idiag_Gamm=0 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 ! where(cnamev=='ism_cool'.or.cnamev=='ism_netcool') cformv='DEFINED' ! 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_netcool') 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)%feat%t=t SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%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)%feat%t=t SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%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)%feat%t=t SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%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)%feat%t=t SNRs(iSNR)%feat%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)%feat%t=0 SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%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)%feat%t=0. SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%radius=width_SN if (luniform_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 ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine interstellar_before_boundary !***************************************************************************** subroutine heat_interstellar(f,zheat) ! ! This routine calculates a vertical profile for uv-heating designed to ! satisfy an initial condition with heating and cooling approximately balanced ! for an isothermal hydrostatic equilibrium. ! Requires: gravz_profile='Ferriere' in gravity_simple.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 ! 13-jul-15/fred: requires initial_condition/hs_equilibrium_ism.f90 ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f real, dimension(mz), intent(out) :: zheat ! real, dimension(mz) :: lambda=0.0, lnTT, zrho real :: logrho, 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 a_S = a_S_cgs/unit_velocity*unit_time a_D = a_D_cgs/unit_velocity*unit_time z_D = z_D_cgs/unit_length z_S = z_S_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 if (lthermal_hse) then logrho = log(rho0ts)+(a_S*z_S*m_u*mu/k_B/T0hs)*(log(T0hs)- & log(T0hs/(a_S*z_S)* & (a_S*sqrt(z_S**2+(z(n))**2)+0.5*a_D*(z(n))**2/z_D))) else logrho = log(rho0ts)-0.015*(- & a_S*z_S+ & a_S*sqrt(z_S**2+(z(n))**2)+0.5*a_D*(z(n))**2/z_D) endif if (logrho < -40.0) logrho=-40.0 zrho(n)=exp(logrho) TT=T0hs/(a_S*z_S)* & (a_S*sqrt(z_S**2+(z(n))**2)+0.5*a_D*(z(n))**2/z_D) lnTT(n)=log(TT) 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.lnTT=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 real, dimension (mz) :: TT ! ! 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 tidy_SNRs if (lSNI) call check_SNI (f,l_SNI) if (lSNII) call check_SNII(f,l_SNI) ! 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)%feat%t=t SNRs(iSNR)%indx%SN_type=1 SNRs(iSNR)%feat%radius=width_SN try_count=10 ! do while (try_count>0) ierr=iEXPLOSION_OK try_count=try_count-1 ! if (luniform_zdist_SNI) then call position_SN_uniformz(f,SNRs(iSNR)) else call position_SN_gaussianz(f,h_SNI,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)%feat%t=t SNRs(iSNR)%indx%SN_type=2 SNRs(iSNR)%feat%radius=width_SN try_count=10 ! do while (try_count>0) ierr=iEXPLOSION_OK try_count=try_count-1 ! if (luniform_zdist_SNI) then call position_SN_uniformz(f,SNRs(iSNR)) else call position_SN_gaussianz(f,h_SNII,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 ! real, dimension(mx,my,mz,mfarray) :: f real :: t_interval, surface_massII 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 ! 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. ! ! 13-jul-15/fred: ! Location by mass was found to lose too much energy due to the numerically ! necessarily high thermal conductivity coefficients. SNe in OBs mainly explode ! into diffuse bubbles left by their neighbours anyway. This routine left for ! reference and possible later applications for clustering and feedback through ! star formation, which to date has been neglected ! l_SNII_guassian skips this algorithm and switches to check_SNIIb. ! 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)%feat%t=t SNRs(iSNR)%indx%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%indx%l,SNR%indx%m,SNR%indx%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%indx%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%indx%ipy=(i-1)/ny ! uses integer division SNR%indx%m=i-(SNR%indx%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%indx%ipz=(i-1)/nz ! uses integer division SNR%indx%n=i-(SNR%indx%ipz*nz)+nghost SNR%indx%iproc=SNR%indx%ipz*nprocy + SNR%indx%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 :: rhomax, maxrho, rhosum real :: mpirho, mpiz real, dimension(ncpus):: tmp2 integer :: yzproc, itmp, icpu, lm_range integer :: previous_SNl, previous_SNm, previous_SNn ! ! 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%indx%l,SNR%indx%m,SNR%indx%n). ! call random_number_wrapper(fran3) ! ! Get 3 random numbers on all processors to keep rnd. generators in sync. ! ! 13-jul-15/fred: NB need to revisit OB clustering x,y not updated or time ! constrained. May need to include z also ! if (lroot) then if (lOB_cluster .and. h_SN==h_SNII) then if (t < t_cluster) then ! still using current cluster coords previous_SNl = int(( x_cluster - xyz0(1) )/Lxyz(1))*nxgrid +1 previous_SNm = int(( y_cluster - xyz0(2) )/Lxyz(2))*nygrid +1 previous_SNn = int(( z_cluster - xyz0(3) )/Lxyz(3))*nzgrid +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%indx%ipx=(i-1)/nx ! uses integer division SNR%indx%l=i-(SNR%indx%ipx*nx)+nghost ! i=int(fran3(1)*lm_range/p_OB)+previous_SNm+1 SNR%indx%ipy=(i-1)/ny ! uses integer division SNR%indx%m=i-(SNR%indx%ipy*ny)+nghost ! i=int(fran3(1)*lm_range/p_OB)+previous_SNn+1 SNR%indx%ipz=(i-1)/nz ! uses integer division SNR%indx%n=i-(SNR%indx%ipz*nz)+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%indx%ipx=(i-1)/nx ! uses integer division SNR%indx%l=i-(SNR%indx%ipx*nx)+nghost ! i=int(fran3(1)*(nygrid-lm_range)/(1.0-p_OB))+previous_SNl+1 if (i>nygrid) i=i-nygrid SNR%indx%ipy=(i-1)/ny ! uses integer division SNR%indx%m=i-(SNR%indx%ipy*ny)+nghost ! 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) 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) 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%feat%radius*mass_width_ratio width_velocity = remnant%feat%radius*velocity_width_ratio radius2 = (remnant%feat%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%feat%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%feat%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%feat%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%feat%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 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=="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) ! 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 !*********************************************************************** subroutine injectfcr_SN(deltafcr,width,cfcr_SN) ! real, intent(in) :: width,cfcr_SN real, intent(out), dimension(nx,3) :: deltafcr ! real, dimension(nx) :: profile_SN ! integer :: j ! ! Calculate deltauu. ! if (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) ! endif ! do j=1,3 deltafcr(1:nx,j)=cfcr_SN*profile_SN(1:nx)* & kperp * outward_normal_SN(1:nx,j) ! spatial mass density enddo ! endsubroutine injectfcr_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)%indx%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)%indx%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)%indx%state==SNstate_invalid) then if (lroot) print*,"Tried to free an already invalid SNR" return endif ! nSNR=nSNR-1 SNRs(iSNR)%indx%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)%indx%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/reduced 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 ! ! Skip this subroutine if not selected eg before turbulent pressure settles ! if (.not. ladd_massflux) return ! ! Only add boundary mass at intervals to ensure flux replacements 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=1.0 can be varied to regulate mass levels as required. ! Replaces previous MPI heavy algorithm to calculate and store actual boundary ! flux. Verified that mass loss matched boundary loss, rather than numerical, ! so sufficient to monitor rhom and adjust addrate to maintain mass. ! add_ratio=1.0+prec_factor*addrate ! ! Add mass proportionally to the existing density throughout the ! volume to replace that lost through boundary. ! 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 !!***************************************************************************** endmodule Interstellar