c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine read_line_parameters(dr,Nmol,mol_index,mol_niso, & iso_index,nfiles,filename,fformat,pct,rl,pr, & Hiso_abundance,iso_abundance, & Nlines, & dbr,lastf,lastl,raz) implicit none c c Purpose: to read spectral line parameters from the HITRAN database, c for a specified number of lines c c Inputs: c + dr: option for LBL data reorganization c + Nmol: number of molecules present at level i c + mol_index: indexes of molecules present at level i c + mol_niso: number of isotopes for each molecule c + iso_index: indexes (using HITRAN convention) of isotopes c + nfiles: number of spectrocopic database files to read c + filename: names of the "nfiles" files to read c + fformat: format of each one of the "nfiles" files c + pct: option "print computation times" c + rl: option for line selection c + pr: percentage of lines rejected when selection is enabled c + Hiso_abundance(mol,iso) is the abundance of isotope 'iso' for molecule 'mol' in terrestrial air c + iso_abundance(mol,iso) is the abundance of isotope 'iso' for molecule 'mol' for the considered gas mixture c c Outputs: c + Nlines: number of lines that have been read c + index: molecule number [unitless] c + isotope: isotopologue number [unitless] c + nu0: vacuum wavenumber [cm¯¹] c + Sref: line intensity [cm¯¹/molecule.cm¯²] at standard Tref c + gamma_air: air-broadened half-width [cm¯¹.atm¯¹] c + gamma_self: self-broadened half-width [cm¯¹.atm¯¹] c + Elow: lower-state energy [cm¯¹] c + n_air: temperature-dependance exponent for gamma_air [unitless] c + delta_air: air pressure-induced line shift [cm¯¹.atm¯¹] c + n_self: temperature-dependance exponent for gamma_self [unitless] c c Input/output c + dbr: 1 if all LBL databases have been read; 0 otherwise c + lastf: index of the last file opened in the case dbr=0 c + lastl: index of the last line read in the last file opened, in the case dbr=0 c include 'max.inc' include 'formats.inc' include 'arrays1.inc' integer line,ios,j,nl c Inputs integer dr integer Nmol integer mol_index(1:Nmol_max) integer mol_niso(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) integer nfiles character*(Nchar_mx) filename(1:nfiles_mx) integer fformat(1:nfiles_mx) logical pct,rl double precision pr double precision Hiso_abundance(1:Nmol_max,1:Niso_max) double precision iso_abundance(1:Nmol_max,1:Niso_max) c Outputs integer Nlines c In/out integer dbr,lastf,lastl,raz c tmp integer nlrif integer ifile,index_tmp,isotope_tmp,lline integer Ierr_tmp(1:6),Iref_tmp(1:6) double precision nu0_tmp,Sref_tmp,A_tmp,gamma_air_tmp double precision gamma_self_tmp,Elow_tmp,n_air_tmp,n_self_tmp double precision delta_air_tmp,gup_tmp,glow_tmp character*15 Vup_tmp,Vlow_tmp,Qup_tmp,Qlow_tmp integer Vup_tmp2,Vlow_tmp2 character*9 Qup_tmp2,Qlow_tmp2 integer molec,isotop double precision tmp c Test character*(Nchar_mx) fichier,nlines_file integer*8 date1,date2 c Test character*1 star_tmp,zero_ch character*(Nchar_mx) file_database integer iso,isof,l1,l1f,l integer strlen character*(Nchar_mx) label label='subroutine read_line_parameters' c Debug c write(*,*) '--------------------------' c write(*,*) label(1:strlen(label)),' :in' c write(*,*) 'lastf=',lastf c write(*,*) 'lastl=',lastl c write(*,*) '--------------------------' c Debug write(zero_ch,11) 0 c Debug c write(*,*) 'This is routine read_line_parameters' c write(*,*) 'nfiles=',nfiles c write(*,*) 'lastf=',lastf,' lastl=',lastl c Debug Nlines=0 do ifile=lastf,nfiles file_database=filename(ifile) nlines_file='./nlines' call get_nlines(file_database,nlines_file,nl) c Debug c write(*,*) 'file:',file_database(1:strlen(file_database)) c write(*,*) 'nl=',nl c Debug open(10,file=file_database(1:strlen(file_database)) & ,status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'File could not be opened:', & file_database(1:strlen(file_database)) stop else if (dr.eq.1) then write(*,*) 'reading and reordering data from:' else write(*,*) 'reading LBL database:' endif write(*,*) file_database(1:strlen(file_database)) c & ,fformat(ifile) endif c Debug c write(*,*) 'Nlines=',Nlines,' nl=',nl,' Nlines+nl=',Nlines+nl c Debug if ((ifile.eq.lastf).and.(lastl.gt.0)) then do line=1,lastl read(10,60) enddo lline=nl-lastl nlrif=lastl ! number of lines already read in current file else lline=nl nlrif=0 ! number of lines already read in current file endif do line=1,lline if (fformat(ifile).eq.1) then ! HITRAN file read(10,60) index_tmp ! molecule number [unitless] & ,isotope_tmp ! isotopologue number [unitless] & ,nu0_tmp ! vacuum wavenumber [cm¯¹] & ,Sref_tmp ! line intensity [cm¯¹/molecule.cm¯²] at standard Tref & ,A_tmp ! Einstein A coefficient [s¯¹] & ,gamma_air_tmp ! air-broadened half-width [cm¯¹.atm¯¹] & ,gamma_self_tmp ! self-broadened half-width [cm¯¹.atm¯¹] & ,Elow_tmp ! lower-state energy [cm¯¹] & ,n_air_tmp ! temperature-dependance exponent for gamma_air [unitless] & ,delta_air_tmp ! air pressure-induced line shift [cm¯¹.atm¯¹] c & ,Vup_tmp ! upper-state "global" quanta [unknown] c & ,Vlow_tmp ! lower-state "global" quanta [unknown] c & ,Qup_tmp ! upper-state "local" quanta [unknown] c & ,Qlow_tmp ! lower-state "local" quanta [unknown] c & ,(Ierr_tmp(j),j=1,6) ! uncertainty index [unknown] c & ,(Iref_tmp(j),j=1,6) ! reference index [unknown] c & ,star_tmp ! flag [character] c & ,gup_tmp ! statistical weight of the upper state [unknown] c & ,glow_tmp ! statistical weight of the lower state [unknown] n_self_tmp=n_air_tmp else if (fformat(ifile).eq.2) then ! HITEMP file read(10,70) index_tmp ! molecule number [unitless] & ,isotope_tmp ! isotopologue number [unitless] & ,nu0_tmp ! vacuum wavenumber [cm¯¹] & ,Sref_tmp ! line intensity [cm¯¹/molecule.cm¯²] at standard Tref & ,A_tmp ! Einstein A coefficient [s¯¹] & ,gamma_air_tmp ! air-broadened half-width [cm¯¹.atm¯¹] & ,gamma_self_tmp ! self-broadened half-width [cm¯¹.atm¯¹] & ,Elow_tmp ! lower-state energy [cm¯¹] & ,n_air_tmp ! temperature-dependance exponent for gamma_air [unitless] c & ,delta_air_tmp ! air pressure-induced line shift [cm¯¹.atm¯¹] c & ,Vup_tmp2 ! upper-state "global" quanta [unknown] c & ,Vlow_tmp2 ! lower-state "global" quanta [unknown] c & ,Qup_tmp2 ! upper-state "local" quanta [unknown] c & ,Qlow_tmp2 ! lower-state "local" quanta [unknown] c & ,(Ierr_tmp(j),j=1,3) ! uncertainty index [unknown] c & ,(Iref_tmp(j),j=1,3) ! reference index [unknown] n_self_tmp=n_air_tmp else if (fformat(ifile).eq.3) then ! CDSD file read(10,61) index_tmp ! molecule number [unitless] & ,isotope_tmp ! isotopologue number [unitless] & ,nu0_tmp ! vacuum wavenumber [cm¯¹] & ,Sref_tmp ! line intensity [cm¯¹/molecule.cm¯²] at standard Tref & ,A_tmp ! Einstein A coefficient [s¯¹] & ,gamma_air_tmp ! air-broadened half-width [cm¯¹.atm¯¹] & ,gamma_self_tmp ! self-broadened half-width [cm¯¹.atm¯¹] & ,Elow_tmp ! lower-state energy [cm¯¹] & ,n_air_tmp ! temperature-dependance exponent for gamma_air [unitless] & ,delta_air_tmp ! air pressure-induced line shift [cm¯¹.atm¯¹] & ,n_self_tmp ! temperature-dependance exponent for gamma_self [unitless] c Debug c n_self_tmp=n_air_tmp c Debug endif nlrif=nlrif+1 ! number of lines already read in current file isof=0 call indexes(Nmol,mol_index,iso_index, & mol_niso,index_tmp,isotope_tmp,molec,isotop,isof) c Debug c if (isof.eq.1) then c write(*,*) index_tmp c & ,nu0_tmp c & ,isotope_tmp c & ,molec,isotop c é ,iso_index(molec,isotop) c stop c endif c Debug c Debug c do iso=1,mol_niso(molec) c if (isotope_tmp.eq.iso_index(molec,isotop)) then c isof=1 c goto 124 c endif c enddo c 124 continue c Debug if (isof.eq.0) then goto 321 endif c c Modification - 2013-03-06 c Reference line intensities must be corrected to account for the fact that the reference intensity Sref c already takes into account the isotopic abundances: Sref is proportional to Ia, the abundance of c the isotopologue in terrestrial air. if ((Hiso_abundance(index_tmp,isotope_tmp).gt.0.0D+0).and. & (iso_abundance(molec,isotop).gt.0.0D+0)) then Sref_tmp=Sref_tmp & *iso_abundance(molec,isotop) & /Hiso_abundance(index_tmp,isotope_tmp) c Debug c Activate only for terrestrial air c tmp=iso_abundance(molec,isotop) c & /Hiso_abundance(index_tmp,isotope_tmp) c if (dabs(tmp-1.0D+0).gt.1.0D-8) then c call error(label) c write(*,*) 'tmp=',tmp c write(*,*) 'iso_abundance(',molec,',',isotop,')=', c & iso_abundance(molec,isotop) c write(*,*) 'Hiso_abundance(',index_tmp, c & ',',isotope_tmp,')=', c & Hiso_abundance(index_tmp,isotope_tmp) c stop c endif c Debug endif c Modification - 2013-03-06 c if (dr.eq.1) then ! dr=1 : on-the-fly data reorganization c If only 1 LBL database file has to be read, no need to sort: just append each line at the end of the array if (nfiles.eq.1) then l1=Nlines+1 c If there is more than 1 LBL database file to read, there is a need to sort lines else l1f=0 if (Nlines.eq.0) then l1=1 l1f=1 else do l=1,Nlines if (nu0_tmp.lt.nu0(l)) then l1=l l1f=1 goto 111 endif enddo ! l endif 111 continue if (l1f.eq.0) then l1=Nlines+1 endif if (l1.le.Nlines) then do l=Nlines,l1,-1 index(l+1)=index(l) isotope(l+1)=isotope(l) nu0(l+1)=nu0(l) Sref(l+1)=Sref(l) gamma_air(l+1)=gamma_air(l) gamma_self(l+1)=gamma_self(l) Elow(l+1)=Elow(l) n_air(l+1)=n_air(l) delta_air(l+1)=delta_air(l) n_self(l+1)=n_self(l) enddo ! l endif endif ! nfiles=1 or not index(l1)=index_tmp isotope(l1)=isotope_tmp nu0(l1)=nu0_tmp Sref(l1)=Sref_tmp gamma_air(l1)=gamma_air_tmp gamma_self(l1)=gamma_self_tmp Elow(l1)=Elow_tmp n_air(l1)=n_air_tmp delta_air(l1)=delta_air_tmp n_self(l1)=n_self_tmp Nlines=Nlines+1 else if (dr.eq.2) then ! dr=2 : simply reading LBL database, reorganization will occur later Nlines=Nlines+1 index(Nlines)=index_tmp isotope(Nlines)=isotope_tmp nu0(Nlines)=nu0_tmp Sref(Nlines)=Sref_tmp gamma_air(Nlines)=gamma_air_tmp gamma_self(Nlines)=gamma_self_tmp Elow(Nlines)=Elow_tmp n_air(Nlines)=n_air_tmp delta_air(Nlines)=delta_air_tmp n_self(Nlines)=n_self_tmp endif ! dr if (Nlines.eq.Nline_mx) then c write(*,*) 'Error from routine read_line_parameters:' c write(*,*) 'Nline_mx has been reached' c stop c Debug c write(*,*) '*********************' c write(*,*) 'ok' c write(*,*) '*********************' c Debug if (line.eq.nl) then if (ifile.eq.nfiles) then dbr=1 goto 321 else dbr=0 lastf=ifile+1 lastl=0 goto 322 endif else dbr=0 lastf=ifile lastl=nlrif goto 322 endif endif ! Nlines>Nline_mx 321 continue enddo ! line close(10) c Debug c write(*,*) 'Nlines=',Nlines c Debug enddo ! ifile dbr=1 ! once all files have been read, etc write(*,*) '...done' 322 continue c Debug c write(*,*) 'lastf=',lastf,' lastl=',lastl c write(*,*) 'Nlines=',Nlines c Debug c In the case LBL data was read without on-the-fly reorganization, c and if more than 1 LBL data file was read, c there is a need to sort lines with respect to increasing nu0 if ((dr.eq.2).and.(nfiles.gt.1)) then if (pct) then fichier='./date' call get_date(fichier,date1) c command="echo $(date +%s) > "//fichier(1:strlen(fichier)) c call exec(command) c open(11,file=fichier(1:strlen(fichier))) c read(11,*) date1 c close(11) endif write(*,*) 'Now sorting lines ...' c ------------------------------ c 05/09/2008 c ------------------------------ c c I now use "ssort" with array a(1:Nline_mx,2) instead of a(1:Nline_mx,9) c because of memory issues. One major inconvenience is that it is now impossible to sort c all data in one call of "ssort": it has to be called for every data that has to c be sorted, which makes things (a little bit) slower. c c Before I had this: c c call data2a(Nlines,nu0,index,isotope,Sref,gamma_air, c & gamma_self,Elow,n_air,delta_air,a) c nc=9 c call ssort(a,Nlines,nc) c call a2data(a,Nlines,nu0,index,isotope,Sref,gamma_air, c & gamma_self,Elow,n_air,delta_air) c c I had to write that: c c Debug c write(*,*) 'ok0' c Debug call sort_all_data1(Nlines) c Debug c write(*,*) 'ok1' c Debug c ------------------------------ c 05/09/2008 c ------------------------------ if (pct) then fichier='./date' call get_date(fichier,date2) write(*,*) 'line sorting took:',date2-date1,' s' endif endif c Debug c write(*,*) 'ok1' c Debug c call save_tmpdata(Nlines,index,isotope,nu0,Sref,gamma_air, c & gamma_self,Elow,n_air,delta_air,n_self) c$$$ c$$$c lines intensity analysis (and line dumping if required) c$$$ if (rl.eq.1) then ! if line rejection is required c$$$ call S_analysis(pr,Nlines,nu0,index,isotope,Sref, c$$$ & gamma_air,gamma_self,Elow,n_air,delta_air) c$$$ endif c Debug c write(*,*) 'ok2' c Debug c Debug c write(*,*) '--------------------------' c write(*,*) label(1:strlen(label)),' :out' c write(*,*) 'dbr=',dbr c write(*,*) 'lastf=',lastf c write(*,*) 'lastl=',lastl c write(*,*) 'Nlines=',Nlines c write(*,*) '--------------------------' c Debug return end subroutine S_analysis(pr,Nlines) implicit none include 'max.inc' include 'arrays1.inc' include 'arrays2.inc' c c Purpose: to analyse lines intensities, and dump some lines if required c c I/O double precision pr integer Nlines c temp integer Nlines2 integer line double precision S0 integer strlen character*(Nchar_mx) label label='subroutine S_analysis' c find S0 call Slimit(Nlines,Sref,pr,S0) c lines selection Nlines2=0 do line=1,Nlines if (Sref(line).gt.S0) then ! line is retained Nlines2=Nlines2+1 nu02(Nlines2)=nu0(line) index2(Nlines2)=index(line) isotope2(Nlines2)=isotope(line) Sref2(Nlines2)=Sref(line) gamma_air2(Nlines2)=gamma_air(line) gamma_self2(Nlines2)=gamma_self(line) Elow2(Nlines2)=Elow(line) n_air2(Nlines2)=n_air(line) delta_air2(Nlines2)=delta_air(line) n_self2(Nlines2)=n_self(line) endif enddo ! line c statistics write(*,*) '**************************' write(*,*) 'Number of lines rejected:',(Nlines-Nlines2), & ' /',Nlines write(*,*) 'Percentage of lines rejected:', & dble(Nlines-Nlines2)/dble(Nlines) write(*,*) '**************************' c outputs Nlines=Nlines2 do line=1,Nlines nu0(line)=nu02(line) index(line)=index2(line) isotope(line)=isotope2(line) Sref(line)=Sref2(line) gamma_air(line)=gamma_air2(line) gamma_self(line)=gamma_self2(line) Elow(line)=Elow2(line) n_air(line)=n_air2(line) delta_air(line)=delta_air2(line) n_self(line)=n_self2(line) enddo ! line c -------------------------------------- return end subroutine Slimit(Nlines,Sref,pr,S0) implicit none include 'max.inc' c c Purpose: to identify "S0", the maximum intensity of lines that will be kept c c I/O integer Nlines double precision Sref(1:Nline_mx),pr double precision S0 c temp integer line,Sminf,nl double precision Smin,Smax,S1,S2,dex,ex integer S0f double precision ex1,ex2 integer i,j,Nbins,bf,b parameter(Nbins=1000) double precision chart(1:Nbins+1,1:3) integer strlen character*(Nchar_mx) label label='subroutine Slimit' c Debug c write(*,*) 'Nlines=',Nlines c do line=1,Nlines c write(*,*) 'Sref(',line,')=',Sref(line) c enddo c write(*,*) 'pr=',pr c Debug c -------------------------------------- c Find Smin and Smax Sminf=0 do line=1,Nlines if (Sref(line).gt.0.0D+0) then Sminf=1 Smin=Sref(line) goto 321 endif enddo 321 continue if (Sminf.eq.0) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'No line with non-null intensity was found' c stop S0=0.0D+0 goto 666 endif Smax=0.0D+0 do line=1,Nlines if ((Sref(line).lt.Smin).and.(Sref(line).gt.0.0D+0)) then Smin=Sref(line) endif if (Sref(line).gt.Smax) then Smax=Sref(line) endif enddo ! line c Debug c write(*,*) 'Smin=',Smin c write(*,*) 'Smax=',Smax c Debug c -------------------------------------- c -------------------------------------- c chart c bins limits S2=Smax*2.0D+0 S1=Smin/2.0D+0 call exponent(S1,ex1) call exponent(S2,ex2) dex=(ex2-ex1)/dble(Nbins) c Debug c write(*,*) 'ex1=',ex1 c write(*,*) 'ex2=',ex2 c write(*,*) 'dex=',dex c Debug do i=0,Nbins ex=ex1+dex*i chart(i+1,1)=10.0D+0**ex c chart(i,1) and chart(i+1,1) are the intensity limits for bin index "i". enddo do i=1,Nbins chart(i,2)=0.0D+0 enddo c sorting lines into bins c write(*,*) 'Analysing the distribution of line intensities...' nl=0 do line=1,Nlines if (Sref(line).gt.0.0D+0) then call exponent(Sref(line),ex) b=int((ex-ex1)/dex)+1 bf=1 if (bf.eq.0) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'bin could not be found for:' write(*,*) 'Sref(',line,')=',Sref(line) stop else chart(b,2)=chart(b,2)+1.0D+0 c chart(i,2) counts the number of lines which intensity is greater that chart(i,1) c and lower than chart(i+1,1) endif nl=nl+1 endif enddo c write(*,*) '...done' c normalization do i=1,Nbins c chart(i,2) is now a distribution chart(i,2)=chart(i,2)/dble(nl) enddo c cumulative chart(1,3)=chart(1,2) do i=2,Nbins c chart(i,3) is the cumulative distribution chart(i,3)=chart(i-1,3)+chart(i,2) enddo c normalization if (chart(Nbins,3).ne.1.0D+0) then c write(*,*) 'warning: chart(3)=',chart(Nbins,3) do i=1,Nbins chart(i,3)=chart(i,3)/chart(Nbins,3) enddo ! i endif c -------------------------------------- c -------------------------------------- c reject weaker lines c look for mininum intensity S0 S0f=0 do i=1,Nbins if (chart(i,3).ge.pr) then S0f=1 S0=chart(i,1) goto 124 endif enddo 124 continue if (S0f.eq.0) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'minimum intensity could not be found for:' write(*,*) 'percentage of lines rejected=',pr write(*,*) 'chart(',Nbins,')=',(chart(Nbins,j),j=1,3) stop endif 666 continue return end subroutine exponent(val,ex) implicit none c c Purpose: to find x in val=10^x c double precision val,ex ex=dlog(val)/dlog(10.0D+0) return end