! This program converts those decompressed fort.7001.1xxxx files included in the list.of.tgz.2 back to co2list or co2hit line list
! convert.f90 does not have line profile, while this convertLS.f90 does have 4 line profile parameters used by HITRAN format
! real*16 and integer*8 are explicitly used, so do not enforce -r8 or -i4 in compilation.
	program convert_co2_list_line_shape
	implicit integer (i-n)
	implicit double precision (a-h,o-z)
	real*16 sint,reduced,a21,term1,term2
        integer*8 mval
	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,cpqr(-1:1)*1,cef(0:1)
	character lineE*105, line*120, ca21*11,ca22*8, ftgz*19, cidex1*5, cidex2*5,cfill*38
	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

! 2022-11
! define the threshold for line list ouput (sthr0), and for HITRAN format output with line profile (sthr_ls)
       sthr0=1.d-30; sthr_ls=1.d-27

!abundance from HITRAN site
        ab(1)=0.984204d0; ab(2)=0.0110574d0; ab(3)=0.00394707d0; 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
! initialize P/Q/R. and e/f, and cfill
        cpqr(-1)='P'; cpqr(0)='Q'; cpqr(1)='R'
        cef(0)='e'; cef(1)='f'
! cfill is a character string one can add to the end of each line, for HITRAN format only
        cfill(1:38)='     3655532429 9 9 711   123.0  121.0'


!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)=1995.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.2")
! read in first .tgz file to decompress into fort.7000 
	open(220,file='list.of.tgz.2',status='old')
221	read(220,'(A19)',err=222,end=222)ftgz
	nab=0
	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.40000.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-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. 


! the selection below can be customized to accommodate different needs at different temperature: 
! si(1:5)*ab(is) are the intensity at 296K, 1000K, 2000K, 3000K, and 4000K,respectively. 
        smax=maxval(si(1:5))*ab(is)  
        if(smax.gt.sthr0)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)
 
! 2022-11 Line Profile Model shared from HITRAN group (Dr. Hargreaves) in recent line list evaluation
! intensity cutoff (filter) can be customized for different temperature and different spectra ranges. 
! write to 8001 with line shape parameters from Rob, and use s296K instead of intensity at specified temperature, following HITEMP/HITRAN format

        if(smax.gt.sthr_ls)then

          mval=ji; if(jf-ji.eq.1)mval=jf

    ! Gamma_air from Hashemi et al. (2020) - Table 2
      ga0 = 1.906E-1; ga1 = 0.1268E+1; ga2 = -1.799E-2; ga3 = 4.027E-4
      ga4 = 1.415E+1; ga5 = 2.731E-2;  ga6 = -9.286E-4; ga7 = 6.254E-5
      gamma_air = ((ga0+ga1*mval+ga2*mval**2+ga3*mval**3)/(1.0+ga4*mval+ga5*mval**2+ga6*mval**3+ga7*mval**4))

    ! HITRAN2020 n_air values (Hashemi 2020 - Table 2 is mislabelled) 
      ena0 = -2.157E+3; ena1 = 6.104E+3; ena2 = 0.9686E+3; ena3 = 2.276E+1
      ena4 =  6.186E+3; ena5 = 1.987E+3; ena6 = 6.790E-13; ena7 = 3.627E-1
      if (mval.le.1)then
        en_air = 0.69d0
      else
        en_air = ((ena0+ena1*mval+ena2*mval**2+ena3*mval**3)/(1.d0+ena4*mval+ena5*mval**2+ena6*mval**3+ena7*mval**4))
      end if

    ! Gamma_self from Tan et al. (2022) - Table 2
      gs0 = 1.312E-1; gs1 =  1.320E-2; gs2 = -3.851E-4; gs3 = 4.312E-6
      gs4 = 1.396E-1; gs5 = -3.000E-3; gs6 =  2.635E-5; gs7 = 1.954E-7
      gamma_self = ((gs0+gs1*mval+gs2*mval**2+gs3*mval**3)/(1.0+gs4*mval+gs5*mval**2+gs6*mval**3+gs7*mval**4))
      !write(921,'(i3,3f16.5)')mval,gamma_air, gamma_self, en_air

! now we have gamma_air, gamma_self, and n_air. next is air pressure shift
        delta_air=0.d0

! def calc_co2_shifts( v1_f,v1_i, v2_f,v2_i, v3_f,v3_i,  Jf,Ji ):
!    # # # =============================================================================
!    # # #   FROM HARTMAN 2009 TABLE 1
!    # # #   Inputs are the upper and lower vibration quanta. Then J upper and lower.
!    # # # =============================================================================
!    Alp_i_V = [ -4.16e-3,    0.440e-3,   2.730e-3 ]
!    Bet_i_V = [  0.000,      0.194,      0.010    ]
!    Alp_i_R = [  0.362e-3,  -1.460e-3,   1.810e-3 ]
!    Bet_i_R = [  0.000,      0.038,      1.163    ]
!    A_i     = [  0.626,      0.323,      1.400    ]
!    xm      = max(Jf,Ji)
!    Delta_R = 0.00
!    Delta_V = 0.00
!    for p in range(0,3):
!        Delta_R = Delta_R + Alp_i_R[p]* exp(-xm*Bet_i_R[p])
!        Delta_V = Delta_V + Alp_i_V[p]* exp(-xm*Bet_i_V[p])
!    Fvib = A_i[0]*float(v1_f-v1_i)  +  A_i[1]*float(v2_f-v2_i)  +  A_i[2]*float(v3_f-v3_i)
!    Shift_CO2_air = float(Ji-Jf)*Delta_R+Fvib*Delta_V
!    return Shift_CO2_air
        alp_v1=-4.16e-3; alp_v2=0.44e-3; alp_v3=2.73e-3
        bet_v1=0.d0; bet_v2=0.194d0; bet_v3=0.010d0
        alp_r1=0.362e-3; alp_r2=-1.460e-3; alp_r3=1.810e-3
        bet_r1=0.d0; bet_r2=0.038d0; bet_r3=1.163d0
        a_i1=0.626d0; a_i2=0.323d0; a_i3=1.4d0

! note: CDSD polyad quantum numbers exist for J<=150 and E'<24,000 cm-1. need to consider other levels
!       kv1=ie(is,jf,iparf,irf,9)-ie(is,ji,ipari,iri,9)
!       kv2=ie(is,jf,iparf,irf,10)-ie(is,ji,ipari,iri,10)
!       kv3=ie(is,jf,iparf,irf,12)-ie(is,ji,ipari,iri,12)
	iv1a=ie(is,jf,iparf,irf,9); iv2a=ie(is,jf,iparf,irf,10); iv3a=ie(is,jf,iparf,irf,12)
        iv1b=ie(is,ji,ipari,iri,9); iv2b=ie(is,ji,ipari,iri,10); iv3b=ie(is,ji,ipari,iri,12)

! if CDSD polyad vib quanta not avaialble, use ie(...,1:3) for v1v2v3
! and we use ie(...,1:3) to repace the 1st,2nd, and 4th number in CDSD polyad vib quanta labels. 
	if(sum(ie(is,jf,iparf,irf,9:13)).eq.0)then
        if(ef.gt.2.3d4.or.jf.gt.150)then
	  iv1a=ie(is,jf,iparf,irf,1); iv2a=ie(is,jf,iparf,irf,2); iv3a=ie(is,jf,iparf,irf,3)
	  ie(is,jf,iparf,irf,9)=iv1a; ie(is,jf,iparf,irf,10)=iv2a; ie(is,jf,iparf,irf,12)=iv3a
	end if
	end if
        if(sum(ie(is,ji,ipari,iri,9:13)).eq.0)then
        if(ei.gt.2.3d4.or.ji.gt.150)then
          iv1b=ie(is,ji,ipari,iri,1); iv2b=ie(is,ji,ipari,iri,2); iv3b=ie(is,ji,ipari,iri,3)
	  ie(is,ji,ipari,iri,9)=iv1b; ie(is,ji,ipari,iri,10)=iv2b; ie(is,ji,ipari,iri,12)=iv3b
        end if
        end if
! now compute vib quanta changes for air shifts
	kv1=iv1a-iv1b
	kv2=iv2a-iv2b
	kv3=iv3a-iv3b
        
        maxjfji=max(jf,ji)

        delta_r=alp_r1*exp(-maxjfji*bet_r1)+alp_r2*exp(-maxjfji*bet_r2)+alp_r3*exp(-maxjfji*bet_r3)
        delta_v=alp_v1*exp(-maxjfji*bet_v1)+alp_v2*exp(-maxjfji*bet_v2)+alp_v3*exp(-maxjfji*bet_v3)
        fvib=a_i1*kv1 + a_i2*kv2 + a_i3*kv3

        delta_air=(ji-jf)*delta_r + fvib*delta_v

! write out HITRAN format with line profiles to fort.8001 --> co2hit.(xxxx)cm-1
! the cfill string defined at line#33 was made up to match HITRAN format
! it can be replaced by the JPSr information,see below:
        write(cfill,'(2x,i3,i2,i6,2x,i3,i2,i6,2x,5i2)')jf,iparf,irf,ji,ipari,iri,jw(1:5)
        write(8001,2222)iso,om2,si(1)*ab(is),a21,gamma_air,gamma_self,ei,en_air,delta_air, &
            ie(is,jf,iparf,irf,9:13),ie(is,ji,ipari,iri,9:13),cpqr(jf-ji),ji,cef(mod(ji+ipari,2)), cfill
2222  format(i3,f12.6,2es10.3,f5.4,f5.3,f10.4,f4.2,f8.6,5x,5i2,5x,5i2,20x,A1,i3,A1,A38)

        end if


        !if (iso.eq.23)write(7301,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)

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

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

! 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')
	call system('mv ./fort.8001 ./co2hit.'//cidex2//'cm-1')

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

	goto 221
222	continue
	close(220)
	
	end 
