! 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