c Noise interpolation table: character*72 noisefile,csoutfile,sigoutfile real*8 zqso,minlambda,maxlambda,zabsmin,zabsmax c parameter (zqso=3.0,minlambda=3646,maxlambda=1216.*(1.+zqso)) c parameter (zabsmin = 1.5, zabsmax = 1e8)! Min/max redshift considered. c$$$ parameter (zqso=4.5,minlambda=3646,maxlambda=1216.*(1.+zqso)) c$$$ parameter (zabsmin = 1.5, zabsmax = 1e8)! Min/max redshift considered. parameter (zqso=4.5,zabsmin=1.5) parameter (minlambda=1216.*(1.+zabsmin),maxlambda=1216.*(1.+zqso)) c parameter (zabsmax = 1e8)! Min/max redshift considered. C c Output files: C parameter (csoutfile = 'contspec.dat') ! Spectrum C parameter (sigoutfile = 'contspec.sig.dat') ! Noise array C c noise file C parameter (noisefile = 'noise.dat') parameter (noisefile = 'noise.dat', csoutfile='contspec.dat', sigoutfile='contspec.sig.dat') c ----------------------- c Parameter declerations: c ----------------------- logical contspec, addnoise, readfilesplus real*8 sigtonoise, minnoise real*8 pixsize, fwhm real*8 minbother_blue, minbother_red, zresol real*8 ibfactor integer simtype logical doH1, doHe2 logical doC3, doC4, doN5, doO6, doSi3, doSi4, doFe2 real*8 Z_rel, Z_index, maxz_rel, zfilfactor real*8 ZC_rel, ZN_rel, ZO_rel, ZSi_rel, ZFe_rel logical read_ionbal, use_internal_h1, use_internal_he2 character*72 ionbal_rootfile character*72 ssoutfile character*72 simdatafile, datadir, ibdir character*27 tfile logical dotauweighted, dotot, docleanspec logical limsigma logical NoPecVel logical verbose integer nly, ncspectra, nsspectra, nsight, n1d, nveloc integer nproc, itdisp, itlow, ithigh, itnone, itstar, itdark integer itgas, itdone, nmax, nfilesmax, nmaxout integer iseed real*8 tpi, pi c ---------------------- c Parameter definitions: c ---------------------- parameter (simtype = 1) c simtype = 1 : Tom Theuns (TT) code (mod. version of HYDRA) c simtype = 2 : Katz, Hernquist, Weinberg (KHW) code c simtype = 3 : Renyue Cen code (Eulerian sightlines) c simtype = 4 : Springel & Hernquist code (Gadget) parameter (contspec = .false.) ! Continuous spectrum? logical squaregrid parameter (squaregrid=.false.) c Which elements should be included? parameter (doH1 = .true.) parameter (doHe2 = .false.) parameter (doC3 = .false.) parameter (doC4 = .false.) parameter (doN5 = .false.) parameter (doO6 = .false.) parameter (doSi3 = .false.) parameter (doSi4 = .false.) parameter (doFe2 = .false.) c Number of higher order Lyman lines to include: parameter (nly = 1) ! 0 =< nly =< nlyman ( = see below) c Metallicities (relative to solar) c Z = min(Z_rel * (1+delta)^Z_index, maxz_rel), where delta is the c density contrast in baryons (the min. is only take if Z_index is nonzero). c If Z_rel < 0 then the simulation metallicity is multiplied by a c factor abs(Z_rel) c parameter (Z_rel = -3.0) parameter (Z_rel=1.d-3) parameter (Z_index = 0.0) parameter (maxz_rel = 1.) c Filling factor of metal enriched gas (only used if Z_rel > 0 and TT) c A random fraction 0 =< zfilfactor =< 1 will be given a metallicity c (using the prescription given above) parameter (zfilfactor = 1.0) c Relative metallicities (relative to solar) parameter (ZC_rel = 1.0, ZN_rel = 1.0, ZO_rel = 1.0, & ZSi_rel = 1.0, ZFe_rel = 1.0) c Factor to rescale ionizing background with: I_UV = I_UV * ibfactor c parameter (ibfactor = 1.0) c Compute tau-weighted dens/temp? (saved in separate files) parameter (dotauweighted = .true.) c Cycle individual sight lines periodically such that maximum flux is c at the ends? logical docycling parameter (docycling = .true.) c Compute clean spectra? (no noise, no smoothing, no (self-)contamination) c Done for C4 and O6, saved in separate files. parameter (docleanspec = .false.) parameter (dotot = .true.) ! Compute total dens,temp,vel,zmet? parameter (verbose = .true.) ! Print lots of info? parameter (iseed = -138) ! random number seed c Used only if contspec: c ====================== parameter (ncspectra = 50) ! # continuous spectra c parameter (ncspectra = 1) ! # continuous spectra c parameter (minlambda = 3645.24, maxlambda = 7200.0, zqso = 3.62) c qso 1422 c parameter (minlambda = 3645.2,zqso = 3.62) c parameter (zqso=3.62) parameter (pixsize = 0.04, fwhm = 6.6) ! Pix size in A, FWHM in km/s c parameter (zresol = 0.05) ! z-binning of sim-files c 128^3 simulation in 10Mpc box parameter (zresol = 0.1) ! z-binning of sim-files c 64^3 simulation in 10Mpc box c parameter (zresol = 0.18) ! z-binning of sim-files c Noise properties: parameter (addnoise = .false.) c If either sigtonoise =< 0 or minnoise < 0, then sigma(lambda,flux) c is read in from the file "noisefile", defined below. c (sigma is the standard deviation of the Gaussian random variable) c If sigtonoise > 0 and minnoise >= 0, then c sigma = minnoise + (1/sigtonoise - minnoise) * flux c Note that we must have minnoise < 1 / sigtonoise parameter (sigtonoise = -50.) parameter (minnoise = -0.01) c Read files "filesplus"? (otherwise all files in "files" are opened c to read headers, and a file "filesplus" is created). parameter (readfilesplus = .true.) c ====================== c Used only if not contspec: c ========================== c parameter (nsspectra = 1024) ! # spectra parameter (nsspectra = 400) ! # spectra c parameter (simdatafile = 'd1000.1194') parameter (simdatafile = 'indata') parameter (ssoutfile = 'singlespec.dat') ! Output filename. c parameter (simdatafile = 'cosmoL128_z3.00_sun.tps.met') c parameter (simdatafile = 'cosmoL32_z3.00.tps') c parameter (simdatafile = 'cosmoL32_z3.00_sun.tps.met') c ========================== c NOT used if KHW code: c ===================== c parameter (datadir = '/work/schaye/designer/sightlines/sight/') c parameter (datadir = './SightlineSpec/') parameter (datadir = './') c parameter (datadir = '/work/schaye/designer/sightlines/linux') c parameter (datadir = '/work/schaye/feedbacksim/sightlines/linux/') c parameter (datadir = '/work/schaye/feedbacksim/new/linux/') c parameter (datadir = '/work/schaye/cen/') c Read ionization interpolation table? parameter (read_ionbal = .true.) c Use HI ionization computed by code? parameter (use_internal_h1 = .false.) c Use HeII ionization computed by code? (not for SH code) parameter (use_internal_he2 = .false.) c ===================== c Used only if sightlines (and thus also TT code): c ================================================ parameter (nsight = 6) ! Number of sight lines per redshift parameter (n1d = 256) ! Number of gas particles per dimension. c parameter (n1d = 128) ! Number of gas particles per dimension. c parameter (n1d = 64) ! Number of gas particles per dimension. c ================================================ c Used only if KHW code: c ====================== real*8 h100_khw, omega0_khw, xlambda0_khw, omegab0_khw, box100_khw parameter(h100_khw = 0.65) ! h parameter(omega0_khw = 0.4) ! Omega_m parameter(xlambda0_khw = 0.6) ! Omega_Lambda parameter(omegab0_khw = 0.02 / h100_khw**2) ! Omega_b parameter(box100_khw = 11.1111) ! box size (h^{-1} Mpc comoving) c ====================== c Directory containing ionization balance interpol. tables. c parameter (ibdir = '/work/schaye/cloudy94/ionizbal/hm01_Q+G/sun/') c parameter (ibdir = '/work/schaye/cloudy94/ionizbal/hm01_Q+G.d2.75/') parameter (ibdir = 'Backgrounds/') c parameter (ibdir = './') c Root name of files containing ionization balance interpolation tables. c (will be extended automatically, e.g. h1.ionbal_rootfile ) c parameter (ionbal_rootfile = 'hm96qso') c parameter (ionbal_rootfile = 'hm01') parameter (ionbal_rootfile = '') c parameter (ionbal_rootfile = 'plS100.d1e6') parameter (NoPecVel = .false.) ! Set pec. vels to zero? c --------------------- c Numerical parameters: c --------------------- c Min. max. optical depth for inclusion: parameter (minbother_blue = 1e-3) ! lambda_0 =< Ly-alpha c parameter (minbother_blue = 1e-8) ! lambda_0 =< Ly-alpha. When only doing Lyman-alpha, can be very small parameter (minbother_red = 1e-4) ! lambda_0 > Ly-alpha parameter (limsigma = .true.) ! Only follow profiles until negligible? c # Pixels when spectra bits are computed: c (this is also the final number if not contspec) c fine for 10Mpc box c parameter (nveloc = 2000) c for high Z and 40 Mpc box c parameter (nveloc=6*1024) parameter (nveloc = 2048) c parameter (nveloc=1600) c parameter (nveloc = 768) c # processors (I left this in for future parallelization, just keep c it set to 1). parameter (nproc = 1) c Particle types: parameter(itdisp=-3,itlow=-2,ithigh=2, & itnone=-2,itstar=-1,itdark=0,itgas=1,itdone=2) c -------------------- c Maximum array sizes: c -------------------- c Number of particles in simulation: c parameter(Nmax=2*n1d**3) c parameter(Nmax=2*64**3) c parameter(Nmax=2*128**3) c for full 3D simulation parameter(Nmax=2*n1d**3) c parameter(Nmax=50000) ! plenty for sightline simulation c For RC code: c parameter(Nmax=768*1000*6) c Note that TT designer simulation arrays contain first all dm c particles, then the rest --> could save memory. c Number of sightlines files: parameter(nfilesmax = 6000) c Ionization balance interpolation: integer ib_nz_max, ib_nt_max, ib_nd_max, ib_ntot_max parameter (ib_nz_max = 50, ib_nt_max = 80, ib_nd_max = 30) parameter (ib_ntot_max = ib_nz_max * ib_nt_max * ib_nd_max) c Noise interpolation: integer n_nl_max, n_nf_max ! Max # bins in wavelength and flux parameter (n_nl_max = 30, n_nf_max = 10) c ------------------------------- c Declerations and common blocks: c ------------------------------- c Code units: real*8 lunit,munit,tunit,vunit,dunit common /units/ lunit,munit,tunit,vunit,dunit c Ionization balance interpolation parameters: integer*4 ib_nz, ib_nt, ib_nd, ib_ntot real*4 ib_z(ib_nz_max), ib_t(ib_nt_max), ib_d(ib_nd_max) real*4 ib_H1(ib_ntot_max), ib_He2(ib_ntot_max) real*4 ib_C3(ib_ntot_max), ib_C4(ib_ntot_max) real*4 ib_N5(ib_ntot_max), ib_O6(ib_ntot_max) real*4 ib_Si3(ib_ntot_max), ib_Si4(ib_ntot_max), ib_Fe2(ib_ntot_max) integer nznt, iz1, iz2 ! Freq. used. real*4 dz1, dz2 ! Freq. used. common /ionbal/ ib_nz, ib_nt, ib_nd, ib_ntot, & ib_z, ib_t, ib_d, & nznt, iz1, iz2, dz1, dz2, ! Frequently used. & ib_H1, ib_He2, ib_C3, ib_C4, ib_N5, ib_O6, ib_Si3, & ib_Si4, ib_Fe2 common /ibfac/ ibfactor real*8 tau_H1(nveloc,nproc), tau_He2(nveloc,nproc) real*8 tau_C3(nveloc,nproc), tau_C4(nveloc,nproc) real*8 tau_N5(nveloc,nproc), tau_O6(nveloc,nproc) real*8 tau_Si3(nveloc,nproc), tau_Si4(nveloc,nproc) real*8 tau_Fe2(nveloc,nproc) real*8 RhoZ_H1(nveloc,nproc), TempZ_H1(nveloc,nproc) real*8 RhoZ_C4(nveloc,nproc), TempZ_C4(nveloc,nproc) real*8 RhoZ_C3(nveloc,nproc), TempZ_C3(nveloc,nproc) real*8 RhoZ_O6(nveloc,nproc), TempZ_O6(nveloc,nproc) real*8 rho_tot(nveloc,nproc), temp_tot(nveloc,nproc) real*8 veloc_tot(nveloc,nproc), met_tot(nveloc,nproc) common /spectra/ tau_H1, tau_He2, tau_C3, tau_C4, & tau_N5, tau_O6, tau_Si3, tau_Si4, tau_Fe2, & RhoZ_H1, TempZ_H1, RhoZ_C4, TempZ_C4, & RhoZ_O6, TempZ_O6, & RhoZ_C3, TempZ_C3, & rho_tot, temp_tot, veloc_tot, met_tot c HI, RhoZ, TempZ, rho_tot, temp_tot, veloc_tot, met_tot, c 1, 2, 3, 4, 5, 6, 7, c HeII, CIII, CIV, NV, OVI, SiIII, SiIV, FeII c 8, 9, 10, 11, 12, 13, 14, 15 c Max. size of output array to write to file for non-continuous c spectra: parameter (nmaxout = (nveloc + 1) * 8) ! Last digit is # arrays c Simulation parameters real*8 ztime, boxkms, rhocb, atime_sim, densscale logical sightlines common /boxpars/ ztime, boxkms, rhocb, atime_sim, densscale, & sightlines c Keep track of warning logical n1dwarned common /warnings/ n1dwarned ! Printed warning already? c ------------ c Atomic data: c ------------ c Multiplet transitions should be in order of decreasing osc strength! real*8 atom_munit parameter (atom_munit = 1.66054e-24) ! g integer nlyman parameter (nlyman=31) real*8 Lambda_H1(nlyman), f_H1(nlyman) data Lambda_H1 / 1215.6701, 1025.7223, 972.5368, 949.7431, & 937.8035, 930.7483, 926.2257, 923.1504, & 920.9631, 919.3514, 918.1294, 917.1806, & 916.429, 915.824, 915.329, 914.919, & 914.576, 914.286, 914.039, 913.826, & 913.641, 913.480, 913.339, 913.215, & 913.104, 913.006, 912.918, 912.839, & 912.768, 912.703, 912.645 / data f_H1 / 0.416400, 0.079120, 0.029000, 0.013940, & 0.007799, 0.004814, 0.003183, 0.002216, & 0.001605, 0.00120, 0.000921, 7.226e-4, & 0.000577, 0.000469, 0.000386, 0.000321, & 0.000270, 0.000230, 0.000197, 0.000170, & 0.000148, 0.000129, 0.000114, 0.000101, & 0.000089, 0.000080, 0.000071, 0.000064, & 0.000058, 0.000053, 0.000048 / real*8 massH parameter (massH = 1.00794 * atom_munit) real*8 Lambda_He2 data Lambda_He2 / 303.7822 / real*8 f_He2 data f_He2 / 0.416400 / real*8 massHe parameter (massHe = 4.002602 * atom_munit) c Solar metallicities all number relative to hydrogen. real*8 Lambda_C3, f_C3 data Lambda_C3 / 977.020 / data f_C3 / 0.7620 / real*8 Lambda_C4(2), f_C4(2) data Lambda_C4 / 1548.2041, 1550.7812 / c$$$ data f_C4 / 0.190800, 0.095220 / data f_C4 / 0.190800, 0.0 / real*8 massC, ZC_solar parameter (massC = 12.0107 * atom_munit, ZC_solar = 3.55e-4) real*8 Lambda_N5(2), f_N5(2) data Lambda_N5 / 1238.821, 1242.804 / c$$$ data f_N5 / 0.157000, 0.078230 / data f_N5 / 0.157000, 0.0 / real*8 massN, ZN_solar parameter (massN = 14.00674 * atom_munit, ZN_solar = 9.33e-5) real*8 Lambda_O6(2), f_O6(2) data Lambda_O6 / 1031.927, 1037.616 / c$$$ data f_O6 / 0.132900, 0.066090 / data f_O6 / 0.132900, 0.0 / real*8 massO, ZO_solar parameter (massO = 15.9994 * atom_munit, ZO_solar = 7.41e-4) real*8 Lambda_Si3, f_Si3 data Lambda_Si3 / 1206.500 / data f_Si3 / 1.6690 / real*8 Lambda_Si4(2), f_Si4(2) data Lambda_Si4 / 1393.76018, 1402.77291 / c$$$ data f_Si4 / 0.5140, 0.2553 / data f_Si4 / 0.5140, 0. / real*8 massSi, ZSi_solar parameter (massSi = 28.0855 * atom_munit, ZSi_solar = 3.55e-5) real*8 Lambda_Fe2(9), f_Fe2(9) data Lambda_Fe2 / 1144.9379, 1608.45085, 1063.176, 1096.8769, & 1260.533, 1121.975, 1081.8748, 1143.226, & 1125.448 / data f_Fe2 / 0.10500, 0.0619, 0.060000, 0.032000, & 0.025000, 0.020000, 0.014000, 0.013310, & 0.010990 / real*8 massFe, ZFe_solar parameter (massFe = 55.845 * atom_munit, ZFe_solar = 3.24e-5) real*8 ZC_abs, ZN_abs, ZO_abs, ZSi_abs, ZFe_abs parameter (ZC_abs = ZC_rel * ZC_solar) parameter (ZN_abs = ZN_rel * ZN_solar) parameter (ZO_abs = ZO_rel * ZO_solar) parameter (ZSi_abs = ZSi_rel * ZSi_solar) parameter (ZFe_abs = ZFe_rel * ZFe_solar) c ------------------------------ c Physical constants/parameters: c ------------------------------ real*8 Ymass,yNumber parameter (Ymass = 0.24) ! Helium mass fraction c n_He / n_H for primordial gas: parameter (yNumber = Ymass / (1.-Ymass) * massH/massHe) real*8 Boltz, Planck, LightSpeed, ThomsonCross, G parameter (Boltz = 1.3807e-16) ! erg/K parameter (Planck = 6.6261e-27) ! erg s parameter (LightSpeed = 2.9979e10) ! cm/s parameter (ThomsonCross = 6.6525e-25) ! cm^2 parameter (G = 6.673e-8) ! cm^3/g/s^2 real*8 Mpc, yr, Msun, m_H, H0 parameter (Mpc = 3.0856e24) ! cm parameter (yr = 3.1558e7) ! s parameter (Msun = 1.989e33) ! g parameter (m_H = 1.6726e-24) ! g (proton mass) parameter (H0 = 1e7/Mpc) ! s^{-1} (100 km/s/Mpc) c -------------------- c Numerical constants: c -------------------- parameter(tpi = 6.283185307179d0, pi = tpi/2.) c ------------------ c Derived constants: c ------------------ integer npix, maxnvpix, maxnvpixtot c # pixels in lambda space parameter (npix = (maxlambda - minlambda) / pixsize + 1) c The constants below are the maximum allowed array sizes. These c arrays are larger than required. Actual sizes are computed in c specwizard. c # pixels in velocity space parameter (maxnvpix = 3e5) c # pixels in velocity space for Fourier transform (power of 2) parameter (maxnvpixtot = 2**19)