program convert_co2_list
	implicit integer (i-n)
	implicit double precision (a-h,o-z)
	real*16 sint,reduced,a21,term1,term2
	parameter (Nline=100000000)   ! not really used here

	dimension zpe(4), part(4,0:5),jw(5),etp1(0:5),t(0:5),si(0:5)    !t(0) is the user-specified temperature, si(0) is the computed intensity at given temperature. 
	dimension e(4,0:250,1:2,0:40000,3),ie(4,0:250,1:2,0:40000,1:13)   ! 1-5 ames, 6-13 cdsd 
        character fn1*1,fn2*2,fn3*3,fn4*4,fn5*5,ciw*5, ciso(4)*3
	character lineE*105, line*120, ca21*11,ca22*8, ftgz*19, cidex1*5, cidex2*5
	dimension min1(1),min2(1),itp1(6)
        dimension ab(4)   ! this is the abundance
	integer*8 nab(4)  ! count the number of each iso lines

!abundance from HITRAN site
        ab(1)=0.984204d0; ab(2)=0.011057d0; ab(3)=0.003947d0; ab(4)=7.33989d-4
!label of isos
	ciso(1)='626'; ciso(2)='636'; ciso(3)='628'; ciso(4)='627'
! pi and temperature
	pi=dacos(-1.d0); tocm=219474.631482453d0
	t(1)=296.d0; 
	do it=2,5
	  t(it)=1000.d0*(it-1)
	end do

!read in zero point energies and partition function values (40Kcm-1 for 626, 36K cm-1 for 636/628/627)
	zpe=0.d0; part=0.d0
	open(20,file='./zpe_part',status='old')
	read(20,*)
	do i=1,4
	  read(20,*)zpe(i)
	end do
	read(20,*); read(20,*)
	do i=1,4
	  read(20,*)part(i,1:5)
	end do
	close(20)

! read in 4 iso rovibrational levels, create symbolic links first
	call system("ln -s './co2.626.J250.levels.Ames-CDSD.dat.v5' level.1")
        call system("ln -s './co2.636.J200.levels.Ames-CDSD.dat.v5' level.2")
        call system("ln -s './co2.628.J200.levels.Ames-CDSD.dat.v5' level.3")
        call system("ln -s './co2.627.J200.levels.Ames-CDSD.dat.v5' level.4")

	t(0)=1773.d0   ! temperature to run simulation
	e=0.d0; ie=0; beta=1.d0/(t(0)*3.16681d-6)
	do i=1,4
	write(fn1,'(i1)')i
	open(20,file='./level.'//fn1,status='old')
101     read(20,'(A105)',err=102,end=102)lineE
	read(lineE,104)j,ip1,isym1,ir,e(i,j,(3+ip1)/2,ir,1),ie(i,j,(3+ip1)/2,ir,1:5),e(i,j,(3+ip1)/2,ir,2:3) 
104	format(3i3,i6,f13.6,5i3,f8.6,f11.6)
	if(lineE(67:67).ne.' ')read(lineE(64:102),*)ie(i,j,(3+ip1)/2,ir,6:13)
	part(i,0)=part(i,0)+dble(2.d0*j+1)*dexp(-beta*e(i,j,(3+ip1)/2,ir,1)/tocm)
	goto 101
102	continue
	close(20)
	write(*,'(A,A,f10.1,A,f20.4)')ciso(i),' levels have been read in, partition value at ',t(0),' = ',part(i,0)
	end do	

! find the list of .tgz files
	call system("ls fort.7000.*tgz >list.of.tgz")
! read in first .tgz file to decompress into fort.7000 
	open(220,file='list.of.tgz',status='old')
221	read(220,'(A19)',err=222,end=222)ftgz
	nab=0			! to calc the total number of lines (instead of each 10cm-1), comment out this line.
	call system("tar xzf ./"//ftgz)
! now ready to read in the compressed list file (.tgz)
	open(7000,file='./fort.7000.'//ftgz(11:15),status='old')
301	read(7000,203,err=302,end=302)isp,jf,jj2,ica21,irf,iri
203     format(i1,i3,i2,i8,2i5)

! convert isp back to iso# (is) and parity of E' and E" 
	!iparf=(isp-1)/4+1
	iparf=1; if(isp.gt.4)iparf=2; ipari=3-iparf
	is=isp-4*(iparf-1)

! convert jj2 back to delta_J (jdif) and the line importance indicators, jw(1:5) for 5 temperatures
! it = 1 - 296K; 2 - 1000K; 3 - 2000K; 4 - 3000K; 5 - 4000K 
! jw(it) = 0 - can be ignored; jw(it)= 1 - better included
 
	jdif=jj2/32-1; ji=jf-jdif
	jj=jj2-32*(jdif+1)
	jw=0; 
	do i=1,5
	  ii=2**(i-1); ii2=ii*2
	  jw(6-i)=mod(jj,ii2)/ii 
	end do

! here you can insert a filter for specific temperature, if necessary
! e.g. ignore a line if it is not important at 2000K 
!       if(jw(3).eq.0)goto 301

! converted A21 back to original value and exponent 
	a21=dble(int(ica21/100))*1.d-5
	diexp=-1.d0*mod(ica21,100)+10.d0
	a21=a21*10.d0**(diexp)

! You can choose Ames-X01d PES based energies OR CDSD energies (when available, J=0-150, 0-24000 cm-1)
! the .tgz data were sorted by Ames-X01d based wavenumbers. You may re-sort them after using CDSD energies. 
! to use Ames original:
!	ef=e(is,jf,iparf,irf,1); ei=e(is,ji,ipari,iri,1)
! by default to use CDSD levels when available + Ames original for remains
	ef=e(is,jf,iparf,irf,1)+e(is,jf,iparf,irf,3)
	ei=e(is,ji,ipari,iri,1)+e(is,ji,ipari,iri,3)

! 2021-04-13 update: We tried to match part of X01d-PES levels with CDSD2019 (J<150, E'<15,000 cm-1)
!      if the 8th parameter in the ie() array was set to negative, which indicates the match was probably reliable (no guarantee though)
!      It is interesting to find out how much impact the inevitable inconsistency caused by such replacements will have at certain temperature 
!   for high-T test, suggest to start with CDSD energy replacement for ALL J<=150 and E<15000 cm-1 levels. 
!   if runs into some problems, we can enforce the filter below, i.e. to use the ie(,8)<0 & E<15000 cm-1 filter.
!   for potential fixes, we just need to update the energy level list files in future. 
!   but for now, not sure how much difference they will make
        !if(ie(is,jf,iparf,irf,8).gt.0)ef=e(is,jf,iparf,irf,1)
	!if(ie(is,ji,ipari,iri,8).gt.0)ei=e(is,ji,ipari,iri,1)
        !if(ef.gt.15000.d0)ef=e(is,jf,iparf,irf,1)
        !if(ei.gt.15000.d0)ei=e(is,ji,ipari,iri,1)

! Your E' cutoff for co2 626 (and other iso) can be put here, I suggest to try E'=37500 cm-1 and E=40000 cm-1
!       if(is.eq.1.and.ef.gt.36000.d0)goto 301
! by default it will go up to 40,000 cm-1 for 626.  (other 3 isos always up to 36000 cm-1)

!2021-04-13 update: although the 627 simulation appears normal at my end, you may disable 627 lines if it contributes too much noise
!       if(is.eq.4)goto 301

	om2=ef-ei
!2021-06 update: E(CDSD) replacement will generate some negative om2 at the lower end, i.e. original om2 < 30 cm-1, just ignore it for now
	if(om2.le.0.d0)goto 301 

!2021-04-13 update: >99% of the 626 (1E-36) and 636/628/627 (1E-32) intensities & A21 (J<=150,E'<15000 cm-1)
! have been updated with the latest values computed on the Ames-2 PES (published 2016) + 3pt-davtq5zO-30K DMS (to publish)
! no orders of magnitude intensity difference were found for bands, only dozens of outlier lines associated with a few resonance-affected levels. 

!now to compute intensity from zpe/part + om2 and elow, which is corresponding to 100% abundance
	si=0.d0
	do it=0,5    ! 0 for user-specified temperature, line #45 
          beta2=0.69503476d0; partf0=part(is,it); temp=t(it)
          term1=dexp(-ei/(beta2*temp))*(1.d0-dexp(-om2/(beta2*temp)))
	  term2=partf0*om2*om2*2.99792458d10*8.d0*pi
	  xi=A21*term1*dble(2*jf+1)/term2
	  si(it)=xi
	end do
	!stop
! then we write out abundance-adjusted intensity 

	iso=is+20
! old fort.7001 is the format I used in line list generation & mixing & rejecting unnecessary lines
!	write(7001,1213)iso,om2,a21,si(1:5)*ab(is),jf,iparf,irf,ji,ipari,iri,jw(1:5) 
!1213    format(i2,f12.5,6es12.5,i4,i2,i6,i4,i2,i6,5i2)
 
! Please choose your own format to write out whatever quantities that fit your needs 
! so here is just an example for full length, i.e. including everything I can think of
! iso-isotopologue #, 21-24 for 626/636/628/627
! om2 - transition wavenumber
! a21 - Einstein A21
! si(0:5)*ab(is) - 100% line intensities at user-specified and 5 pre-set temperatures are scaled by their terrestial abundances
! jf, iparf, irf - J, parity, root number for E'
! ji, ipari, iri - J, parity, root number for E"
! jw(1:5) - importance indicator for 5 pre-set temperatures, 0 - to ignore, 1 - to include
! ef, ei - E' and E"
! ie(:,:,:,:,1:3) are the v1v2v3 from vtet, ie(:,:,:,:,6:13) are the p (polyad #),c (wang's sym),n(root# in that polyad) and v1v2l2v3r from CDSD
! if you need wang's symmetry:  mod(j+ipar,2)=0 --> e;  mod(j+ipar,2)=1 --> f
! now the first intensity value after A21 is for the temperature user specified at line #45. 

	write(7001,1213)iso,om2,a21,si(0:5)*ab(is),jf,iparf,irf,ji,ipari,iri,jw(1:5), &
             ef,ei,ie(is,jf,iparf,irf,1:3),ie(is,jf,iparf,irf,6:13),ie(is,ji,ipari,iri,1:3),ie(is,ji,ipari,iri,6:13)
1213    format(i2,f12.5,7es12.5,i4,i2,i6,i4,i2,i6,5i2,2f12.5,3i3,1x,2i3,i5,5i3,2x,3i3,1x,2i3,i5,5i3)

! count the number of lines of 4 isotopologues
        nab(is)=nab(is)+1

	goto 301
302	continue
	close(7000)
	close(7001)

! convert .tgz index number to wavenumber (upper limit) for this datafile
	read(ftgz(12:15),*)idex; idex=idex*10; cidex2='00000'
	write(cidex2(1:5),'(i5)')idex
	if(idex.lt.10000)cidex2(1:1)='0'	
        if(idex.lt.1000)cidex2(1:2)='00'
        if(idex.lt.100)cidex2(1:3)='000'
        if(idex.lt.10)cidex2(1:4)='0000'

! save the decompressed list
	call system('mv ./fort.7001 ./co2list.'//cidex2//'cm-1')

	write(*,'(A,A,4i14)')ftgz, ' done', nab(1:4)

	goto 221
222	continue
	close(220)
	
	end