c c......................................................................... c PROGRAM rates c version: RATES 1.0 c release date of this version: Sept 1, 2014 c code author: Chris Churchill (cwc@nmsu.edu) c reference: Churchill, Klimek, Medina, & Vander Vliet c (2014, ApJ, submitted) c main driver routine for computing the equilibrium ionzation c conditions of a gas cloud assuming constant density and constant c temperature for a given metallicity and ionizing SED c includes photoionization, Auger ionization, direct collisional c ionization, excitation auto-ionization, charge exchange c ionization, photo-recombination, dielectronc recombination, and c charge exchange recombination c in this version, no radiative transfer through the cloud model is c incorporated so the optically thin gas is assumed (no ionization c structure in the cloud model) c Three modes are available (defined in 'rates.inp' file) c MODE=1 generate a model for single nH and T c MODE=2 generate a model for single nH as function of T c MODE=3 generate a model as function of both nH and T c all modes assume a fixed metallicity and ionizing SED c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit none include 'rates.h' integer i,ii,j,k,kk,NT,NnH double precision toler double precision z,T,Zmet,ne double precision spop_age,spop_Zmet,spop_Msol,rstar double precision logTlo,logThi,dlogT,Tlo,pow double precision lognHlo,lognHhi,dlognH,nHlo include 'rates.com' include 'com-modes.com' c stamp header; open the .runlog file OPEN(unit=4,file='rates.runlog',status='unknown') WRITE(6,*) ' ' WRITE(6,*) '........................................' WRITE(6,*) '. .' WRITE(6,*) '. RATES V1.0 .' WRITE(6,*) '. .' WRITE(6,*) '. An Ionization Balance Code .' WRITE(6,*) '. .' WRITE(6,*) '. Version 1.0 September 2014 .' WRITE(6,*) '. .' WRITE(6,*) '........................................' WRITE(6,*) ' ' WRITE(4,*) ' ' WRITE(4,*) '........................................' WRITE(4,*) '. .' WRITE(4,*) '. RATES V1.0 .' WRITE(4,*) '. .' WRITE(4,*) '. An Ionization Balance Code .' WRITE(4,*) '. .' WRITE(6,*) '. Version 1.0 September 2014 .' WRITE(4,*) '. .' WRITE(4,*) '........................................' WRITE(4,*) ' ' c configure, communicate, and open files for writing c see comments in individual subroutines for details CALL configall(T,logTlo,logThi,dlogT, & lognHlo,lognHhi,dlognH,Zmet,z, & spop_Msol,spop_age,spop_Zmet,rstar,toler) CALL commconfig(Zmet,z,T,logTlo,logThi,dlogT,NT, & lognHlo,lognHhi,dlognH,toler, & spop_Msol,spop_age,spop_Zmet,rstar) CALL oprintfiles c compute photoionization and Auger rate coefficients; these are c populated only if photoionization is included IF (doPH) CALL alphas_nonT c run MODE=1 c single nH value, single T value, communicate results as directed IF (runmode.eq.1) then NT = 1 ! we are doing only 1 temperature CALL initabund(Zmet) ! compute abundance fractions CALL commcloud(z,T,Zmet,0,0) ! communicate cloud properties CALL alphas_T(T) ! compute T dependent rate coefficients CALL commalphas(T) ! output them to runlog file CALL ionbalance(toler) ! compute equilibrium solution CALL commsol ! output solution to runlog file CALL storefakjT(T,NT) ! populate arrays for printing to files CALL storenentot(NT) ! populate arrays for printing to files IF (printionfracs) CALL printfkjT(NT,Zmet) ! ionization fractions IF (printdensities) CALL printnkjT(NT,Zmet) ! ion number densities IF (printrates) CALL printakjT(NT,Zmet) ! rates coefficients ENDIF c run MODE=2 c single nH value, loop over T, communicate results as directed IF (runmode.eq.2) then CALL initabund(Zmet) NT = (logThi-logTlo)/dlogT + 1 ! set temperature grid Tlo = 10.0d0**logTlo ! set minimum temperature DO i=1,NT ! loop over temperature pow = float(i-1)*dlogT T = Tlo*(10.0d0**pow) CALL commcloud(z,T,Zmet,0,i) CALL alphas_T(T) CALL commalphas(T) CALL ionbalance(toler) CALL commsol CALL storefakjT(T,i) CALL storenentot(i) ENDDO ! next temperature IF (printionfracs) CALL printfkjT(NT,Zmet) IF (printdensities) CALL printnkjT(NT,Zmet) IF (printrates) CALL printakjT(NT,Zmet) ENDIF c run MODE=3 c loop over nH, loop over T, communicate as directed IF (runmode.eq.3) then NnH = (lognHhi-lognHlo)/dlognH + 1 ! set up density grid nHlo = 10.0d0**lognHlo ! set minimum density NT = (logThi-logTlo)/dlogT + 1 ! set up temperatre grid Tlo = 10.0d0**logTlo ! set minimum temperature DO ii=1,NnH ! loop over density pow = float(ii-1)*dlognH nH = nHlo*(10.0d0**pow) CALL initabund(Zmet) DO i=1,NT ! loop over temperature pow = float(i-1)*dlogT T = Tlo*(10.0d0**pow) CALL commcloud(z,T,Zmet,ii,i) CALL alphas_T(T) CALL commalphas(T) CALL ionbalance(toler) CALL commsol CALL storefakjT(T,i) CALL storenentot(i) ENDDO ! next temperature IF (printionfracs) CALL printfkjT(NT,Zmet) IF (printdensities) CALL printnkjT(NT,Zmet) IF (printrates) CALL printakjT(NT,Zmet) ENDDO ! next density END IF c close all the files and write out the SED if directed to do so CALL cprintfiles CLOSE(unit=4) IF (printseds) CALL printSED c fini STOP END c c......................................................................... c c INCLUDE MODULES (subroutines and functions) c c rate coefficients, i/o, and communication include 'rates-alphas.f' include 'rates-files.f' include 'rates-comm.f' c photoionization include 'rates-sed.f' include 'funcs-photo.f' include 'files-photo.f' include 'phfit-mod.f' include 'phfit-data.blk' c recombination include 'funcs-recomb.f' include 'files-recomb.f' c collisional ionization include 'funcs-coll.f' include 'files-coll.f' c charge transfer include 'kingdon-mod.f' include 'kingdon-data.blk' c auger ionization include 'files-auger.f' include 'funcs-auger.f' c general subs, rate matrix, and solver include 'rates-subs.f' include 'rates-solve4.f' c eof