! To extract a HCN/HNC linelist from the Einstein A, coeffients ! of Harris et al. 2002, ApJ, 578,657 and energy levels of ! Harris et al. MNRAS (Submitted). ! ! INPUT FILES: ! ! fort.99 An Einstein A coefficent file from Harris et al. 2002. ! enlvs-combined.txt The energy level file from Harris et al. 2005. ! ! OUTPUT FILE: ! ! linelist.txt The extracted linelist of the format outlined in ! Harris et al. 2005. ! ! The user my wish to chnage the following variables ! ! T: The temperature for which absolute intensity and the partition function ! are computed ! simin: The minimum intensity (cm mol^{-1}) to be output. ! numin: Minimum frequncy output (cm^{-1}) ! numax: Maximum frequncy output (cm^{-1}) ! zai : if true output an ab initio linelist, if false output the linelist ! which uses lab determined energy levels where possible. ! iout : Selects the format of the output linelist. ! 1: output the minimum data necisarry for an opacity calculation ! with Einstein A coeffients (nu, E", J", J', A). ! 2: outputs linelist in format given in Harris et al. (2005), plus the ! integrated absorption intesity at T (cm molecule^{-1}). ! 3: outputs data in compressed Einstien A coeffient file ! (index1, index2, A) ! 4: outputs data in compressed Einstien A coeffient file, with frequency ! (index1, index2, A, nu) ! 5: Outputs linelist in format given in Harris et al. (2005), plus the ! dimensionless weighted oscillator strength (gf). ! 6: output the minimum data necisarry for an opacity calculation ! with weighted oscillatior strengths (nu, E", gf). ! ANY OTHER VALUE: Outputs linelist in format given in Harris et al. (2005). ! zmin : if true output the minimum data necisarry for an opacity calculation ! (nu, E", J", J', A). ! zSi : if .true. and zmin is false, this outputs the integrated absorption ! intesity at T (cm molecule^{-1}) as well as the Einstein A coefficent. ! ! Author G. J. Harris August 2005. implicit none integer in, nen,out parameter (in=10,out=11,nen=168110) ! storeage space for energy level data. real(8), allocatable, dimension(:) :: eai,ex integer, allocatable, dimension(:) :: JJ,Pp,nn,v1,v2,ll,v3,iso ! data real(8) const,k,gfk ! user changable variables real(8) T, simin, numin,numax integer iout logical zai ! local variables real(8) A, Si,nu, El,eu,e1,e2,beta,expfact,Q,intS,gf integer ind1,ind2,iu,il,count character(2) at data const /1.3270914d-12/ data gfk /1.499197d0/ data k /0.6950386d0/ !---------------------------- ! user changeable variables T=3000.0d0 simin=1.0d-26 numin=2000.0d0 numax=3800.0d0 zai=.false. iout=5 ! --------------------------- beta=1.0d0/(k*T) allocate(eai(nen),ex(nen),JJ(nen),Pp(nen),iso(nen)) allocate(nn(nen),v1(nen),v2(nen),ll(nen),v3(nen)) ! read energy levels call read_energy(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3) ! compute partition function call pfcalc(eai,ex,jj,T,Q) write(6,*) "T = ", T, "Q = ", Q ! read transition strength file, line by line. open(unit=in,form='formatted',file='fort.99') open(unit=out,form='formatted',file='linelist.txt') ints=0.0d0 count=0 10 read(in,1001,end=92) ind2, ind1, A if(ex(ind1) .ge. 0.0d0 .and. ex(ind2) .ge. 0.0d0 .and. & .not. zai) then nu=ex(ind2)-ex(ind1) e1=ex(ind1) e2=ex(ind2) at='lb' else nu=eai(ind2)-eai(ind1) e1=eai(ind1) e2=eai(ind2) at='ai' endif if (nu .eq. 0.0d0) goto 10 if(nu .lt. 0.0d0) then il=ind2 iu=ind1 el=e2 eu=e1 nu=-nu else il=ind1 iu=ind2 el=e1 eu=e2 endif if(nu .lt. numin .or. nu .gt. numax) goto 10 expfact=exp(-El*beta)*(1.0d0-exp(-nu*beta)) SI=(const/(Q*(nu**2)))*expfact*A*dble((2*JJ(iu))+1) if(Si .ge. simin) then count=count+1 ints=ints+si if(iout .eq. 1) then write(out,1004) nu,el,JJ(il),JJ(iu),A elseif(iout .eq. 2) then write(out,1002) nu,el,JJ(il),pp(il),nn(il),JJ(iu),pp(iu),& nn(iu),A,Si,il,iu,& iso(il),v1(il),v2(IL),ll(il),v3(il),& iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at elseif(iout .eq. 3) then write(out,1001) iu, il, A elseif(iout .eq. 4) then write(out,1005) iu, il, A, nu elseif(iout .eq. 5) then gf=(gfk*dble((2*JJ(iu))+1)*A) / (nu*nu) write(out,1002) nu,el,JJ(il),pp(il),nn(il),JJ(iu),pp(iu),& nn(iu),A,gf,il,iu,& iso(il),v1(il),v2(IL),ll(il),v3(il),& iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at elseif(iout .eq. 6) then gf=(gfk*dble((2*JJ(iu))+1)*A) / (nu*nu) write(out,1006) nu,el,gf else write(out,1003) nu,el,JJ(il),pp(il),nn(il),JJ(iu),& pp(iu),nn(iu),A,il,iu,& iso(il),v1(il),v2(IL),ll(il),v3(il),& iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at endif endif goto 10 92 close(in) close(out) write(6,*) "Integrated intensity ", ints write(6,*) "number of lines output ", count deallocate(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3) 1001 format(2i7,1p,d14.6) 1002 format(f12.6,f13.6,i3,i2,i5,i3,i2,i5,1p,2e10.3,0p,2i7,1x,5i3,1x,5i3,1x,a2) 1003 format(f12.6,f13.6,i3,i2,i5,i3,i2,i5,1p,e10.3,0p,2i7,1x,5i3,1x,5i3,1x,a2) 1004 format(f12.6,f13.6,2i3,1p,e10.3) 1005 format(2i7,1p,d14.6,0p,f13.6) 1006 format(f12.6,f13.6,e10.3) end ! ----------------------------------------------------------------- subroutine read_energy(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3) ! the primary function of this subroutine is to read and store the ! energy level data. implicit none integer in, nen parameter (in=10,nen=168110) ! arguments real(8) Eai(*), Ex(*) integer JJ(*),Pp(*),nn(*),iso(*),v1(*),v2(*),ll(*),v3(*) ! local variables integer it, jt,pt,nt,isot,v1t,v2t,v3t,lt real(8) eait,ext integer i open(unit=in,form='formatted',file='enlvs-combined.txt') do 10 i=1, nen read(in,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ext ! write(6,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ext JJ(it)=Jt nn(it)=nt pp(it)=pt iso(it)=isot v1(it)=v1t v2(it)=v2t v3(it)=v3t ll(it)=lt eai(it)=eait ex(it)=ext if(it .gt. 1 .and. ex(it) .lt. 1.0d-10) ex(it)=-1.0d0 ! write(6,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ex(it) 10 continue close(in) return 1001 format(i6,i3,i2,i5,f13.6,5i3,f13.6,1p,e11.2,0p,2x,a1) end ! --------------------------------------------- subroutine pfcalc(eai,ex,jj,T,Q) implicit none integer nen parameter (nen=168110) ! arguments real(8), dimension(*) :: eai,ex integer, dimension(*) :: jj real(8) T,Q ! data real(8) k ! local variables real(8) beta, entmp integer i data k /0.6950333d0/ beta=1.0d0/(k*T) Q=0.0d0 do 10 i=1, nen if(ex(i) .lt. 0.0d0) then entmp=eai(i) else entmp=ex(i) endif Q=Q+(dble((2*jj(i))+1)*exp(-entmp*beta)) 10 continue return end ! ---------------------------------------------------