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 max_error(pass,Nb,band,nu_low,nu_hi, & Sref_frac,Tref,T,P,Ps,m,i, & Nmol,mol_index,mol_niso,iso_index, & iso_mass,c2,Rgp,iso_Qref,iso_g, & k_mxep,sd_mxep,pct,sda,constant_dnu, & bslal,nlw_cz,np_lc,np_bl, & sens_P,sens_T,sens_x,pr_sd,large_db, & nk,nu,band_sigma,band_k, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) implicit none c c Purpose: to discretize the provided interval [nu_min,nu_max] into nk nu points c using variable meshing technique based on Voigt and Lorentz profiles tabulation c c Inputs: c + pass: pass index c + Nb: number of narrowbands c + band: index of the narrowband to discretize c + nu_low: lower limit of each narrowband [cm¯¹] c + nu_hi: upper limit of each narrowband [cm¯¹] c + Sref_frac: fraction of max. Sref used for line sorting c + Tref: array of reference temperature c + T: temperature, at each atmospheric level c + P: the total pressure, at each atmospheric level c + Ps: partial pressure of each species, at each level c + m : number of atmospheric levels c + i: index of the current atmospheric level c + Nmol: number of molecules used in the HITRAN database c + mol_index: indexes of the HITRAN molecules c + mol_niso: number of isotopes for each molecule c + iso_index: indexes of the HITRAN isotopes c + iso_mass(mol,iso) are molar masses [g/mol] c + c2: second radiation constant [K.cm] c + Rgp: perfect gas constant [atm.m³/(mol.K)] c + iso_Qref: total internal partition sum at reference temperature (Tref) c + iso_g(mol,iso) are state independent degeneracy factors c + k_mxep: maximum allowed error percentage of k c + sd_mxep: maximum allowed error percentage for spectral discretization c + pct: option "print computation time" c + sda: option for spectral discretization algorithm c + constant_dnu: constant spectral step (inv. cm) c + bslal : option "accuracy level for between-lines spectral meshing" c + nlw_cz: number of line widths to consider in the definition of lines central zones c + np_lc: number of points for line centers discretization c + np_bl: number of points for between-lines discretization c + sens_* : options for computation of sensitivities c + pr_sd: percentage of lines rejected for spectral discretization c + large_db: true if large LBL database is used (CDSD_4000) c c Output: c + nk: number of wavenumber values in the band c + nu: array of wavenumber values [cm¯¹] c + band_sigma & band_k: temporary results in the case of multipass computations c + band_dsigma* & band_dk*: temporary results in the case of multipass computations c include 'max.inc' include 'parameters.inc' include 'formats.inc' include 'arrays3.inc' include 'mpif.h' c Inputs integer pass,Nb,band double precision nu_low(1:Nbmx) double precision nu_hi(1:Nbmx) double precision Sref_frac double precision Tref(1:Nmol_max) double precision T(1:Nmax) double precision P(1:Nmax) double precision Ps(1:Nmol_max) integer m integer i double precision c2,Rgp double precision iso_Qref(1:Nmol_max,1:Niso_max) integer iso_g(1:Nmol_max,1:Niso_max) integer Nmol integer mol_index(1:Nmol_max) integer mol_niso(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) double precision iso_mass(1:Nmol_max,1:Niso_max) double precision k_mxep,sd_mxep logical pct integer sda,bslal double precision constant_dnu integer nlw_cz,np_lc,np_bl logical sens_P,sens_T,sens_x double precision pr_sd logical large_db c Outputs integer nk double precision nu(1:Nkmx) double precision band_sigma(1:Nkmx) double precision band_k(1:Nkmx) double precision band_dsigma_dP(1:Nkmx) double precision band_dsigma_dT(1:Nkmx) double precision band_dsigma_dx(1:Nkmx,1:Nmol_max) double precision band_dk_dP(1:Nkmx) double precision band_dk_dT(1:Nkmx) double precision band_dk_dx(1:Nkmx,1:Nmol_max) c temp integer nk_tmp,mol double precision nu_tmp(1:Nkmx) double precision nu_min,nu_max double precision numin(1:Nchunk_mx) double precision numax(1:Nchunk_mx) integer ik,ik_min integer nxtab,nl,nacc double precision xtab(1:Nxtab_mx,1:3) double precision xacc(1:Nacc_mx,1:2) c integer nk_func,nx(1:Nytab_mx) integer kindex(1:Nytab_mx,1:2) double precision ktab(1:Nktab_mx,1:3) double precision yvalues(1:Nytab_mx) integer*8 date1,date2,cputime character*(Nchar_mx) datep_file,smesh_file character*(Nchar_mx) nlinesp_file character*(Nchar_mx) kfile c double precision a(1:Nline_mx,2) double precision intw integer out c integer Nks double precision sgrid(1:Nkmx) c c MPI integer ierr,np,pindex integer stat(MPI_STATUS_SIZE) integer chunk,nchunks integer code0 integer proc,fin,fpf integer free(1:Nproc_mx) integer ch_idx(1:Nproc_mx) integer nsp c MPI c label integer strlen character*(Nchar_mx) label label='subroutine max_error' CALL MPI_COMM_SIZE(MPI_COMM_WORLD,np,ierr) CALL MPI_COMM_RANK(MPI_COMM_WORLD,pindex,ierr) if (pindex.eq.0) then ! root process only datep_file='./date_sd' call get_date(datep_file,date1) cputime=0 nu_min=nu_low(band) nu_max=nu_hi(band) if (sda.eq.0) then write(*,40) 'Reading grid for interval [',nu_min, & '-',nu_max,'] cm^{-1}, number of points=' else write(*,40) 'Discretization of interval [',nu_min, & '-',nu_max,'] cm^{-1} into...' endif if (pass.eq.1) then ! for the first pass, spectral meshing has to be computed if (sda.eq.0) then call read_sgrid(nu_min,nu_max,nk,nu) if (nk.gt.Nkmx) then call error(label) write(*,*) 'nk=',nk write(*,*) 'while Nkmx=',Nkmx stop endif fin=1 goto 222 endif call read_dx(sd_mxep,nxtab,xtab) call acceleration_dx(nxtab,xtab,nacc,xacc) call read_ktab(sd_mxep,ktab,nk_func,nx,yvalues,kindex) chunk=1 numin(chunk)=nu_min ! cm¯¹ call choose_intw(numin(chunk),nu_max,large_db,intw) numax(chunk)=nu_min+intw ! cm¯¹ do while(numax(chunk).lt.nu_max) chunk=chunk+1 numin(chunk)=numax(chunk-1) ! cm¯¹ call choose_intw(numin(chunk),nu_max,large_db,intw) numax(chunk)=numax(chunk-1)+intw ! cm¯¹ c Debug c write(*,*) chunk,numin(chunk),intw,numax(chunk) c Debug enddo nchunks=chunk code0=3 ! data transfer do proc=1,np-1 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) enddo c send integer data call MPI_BCAST(sda,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(bslal,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nmol,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(i,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(mol_index,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(mol_niso,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(iso_index,Nmol_max*Niso_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(iso_g,Nmol_max*Niso_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nacc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nxtab,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nk_func,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nx,Nytab_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(kindex,Nytab_mx*2,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nlw_cz,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(np_lc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(np_bl,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c send double precision data call MPI_BCAST(constant_dnu,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(c2,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Sref_frac,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(iso_mass,Nmol_max*Niso_max, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Tref,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(T,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(P,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Ps,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(xtab,Nxtab_mx*3,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(xacc,Nacc_mx*2,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ktab,Nktab_mx*3,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(yvalues,Nytab_mx,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(pr_sd,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) do proc=1,np-1 free(proc)=1 enddo nk=0 fin=0 chunk=0 do while (fin.eq.0) 444 continue call find_free_process(free,np-1,fpf,proc) if (fpf.eq.1) then ! a free child process found chunk=chunk+1 if (chunk.le.nchunks) then ch_idx(proc)=chunk code0=4 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(numin(chunk),1,MPI_DOUBLE_PRECISION, & proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(numax(chunk),1,MPI_DOUBLE_PRECISION, & proc,1,MPI_COMM_WORLD,ierr) c Debug c write(*,*) 'computing:',numin(chunk),numax(chunk) c Debug free(proc)=0 ! child process index 'proc' is busy else ! chunk>nchunks free(proc)=2 ! means child process index 'proc' has been stopped endif goto 444 else ! no free child process found call stopped_processes(free,np-1,nsp) if (nsp.eq.np-1) then fin=1 goto 111 endif c wait for results call MPI_RECV(proc,1,MPI_INTEGER & ,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) call MPI_RECV(nk_tmp,1,MPI_INTEGER & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) call MPI_RECV(nu_tmp,Nkmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c Debug c write(*,*) 'nk_tmp=',nk_tmp c Debug free(proc)=1 ! child process index 'proc' is free again if (ch_idx(proc).eq.1) then ik_min=1 else ik_min=2 endif do ik=ik_min,nk_tmp nk=nk+1 nu(nk)=nu_tmp(ik) if (nk.gt.Nkmx) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'Nkmx has been reached' stop endif c Debug call test_nan(nu_tmp(ik),out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu_tmp(',ik,')=',nu_tmp(ik) stop endif c Debug enddo endif ! fpf=0 or 1 111 continue enddo ! while fin=0 222 continue c sorting the 'nu' array: c c security if (nk.gt.Nline_mx) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'when sorting the "nu" array:' write(*,*) 'nk=',nk write(*,*) 'size of temporary "a" array is:' write(*,*) 'Nline_mx=',Nline_mx stop endif c 1 - transfering unsorted 'nu' to array 'a' do ik=1,nk a(ik,1)=nu(ik) a(ik,2)=nu(ik) enddo c 2 - sorting array 'a' call ssort(nk) c 3 - transfering sorted 'a' to array 'nu' do ik=1,nk nu(ik)=a(ik,1) enddo c Test do ik=2,nk if (nu(ik).lt.nu(ik-1)) then write(*,*) 'nu(',ik,')=',nu(ik), & ' < nu(',ik-1,')=',nu(ik-1) stop endif enddo c Test call get_date(datep_file,date2) cputime=cputime+(date2-date1) c Debug if (nu_min.lt.nu(1)) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu_min=',nu_min,' nu(1)=',nu(1) stop endif if (nu_max.gt.nu(nk)) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu_max=',nu_max,' nu(',nk,')=',nu(nk) stop endif c Debug call adjust_mesh(nu_min,nu_max,nu,nk) smesh_file='./results/nu.txt' nlinesp_file='./nlines' call get_nlines(smesh_file,nlinesp_file,nl) open(12,file=smesh_file(1:strlen(smesh_file))) do ik=1,nl read(12,*) enddo if (nl.eq.0) then ik_min=1 else ik_min=2 endif do ik=ik_min,nk write(12,*) ik,nu(ik) enddo close(12) do ik=1,Nkmx band_sigma(ik)=0.0D+0 band_k(ik)=0.0D+0 if (sens_P) then band_dsigma_dP(ik)=0.0D+0 band_dk_dP(ik)=0.0D+0 endif if (sens_T) then band_dsigma_dT(ik)=0.0D+0 band_dk_dT(ik)=0.0D+0 endif if (sens_x) then do mol=1,Nmol_max band_dsigma_dx(ik,mol)=0.0D+0 band_dk_dx(ik,mol)=0.0D+0 enddo ! mol endif enddo else ! for subsequent passes, spectral mesh has already been computed c temporary results have also to be read from temporary results file call result_file_name_tmp(i,band,kfile) nlinesp_file='./nlines' call get_nlines(kfile,nlinesp_file,nk) open(13,file=kfile(1:strlen(kfile))) do ik=1,nk read(13,52) nu(ik),band_sigma(ik),band_k(ik), & band_dsigma_dP(ik),band_dsigma_dT(ik), & (band_dsigma_dx(ik,mol),mol=1,Nmol), & band_dk_dP(ik),band_dk_dT(ik), & (band_dk_dx(ik,mol),mol=1,Nmol) enddo close(13) endif if (nk.gt.Nkmx) then write(*,*) 'Error in routine max_error:' write(*,*) 'Nkmx=',Nkmx,' is not sufficient' write(*,*) 'increase this parameter in max.inc' stop endif write(*,*) '...',nk,' points' if (pct) then write(*,*) 'Spectral discretization took :' call print_cputime(cputime) write(*,41) 'average time per nu point:', & dble(cputime)/dble(nk)*1.0D+3,' ms' endif endif ! pindex 666 continue return end subroutine max_error_1cm_alrl(pindex,numin,numax,c2, & sda,constant_dnu,bslal, & nlw_cz,np_lc,np_bl,Sref_frac,Nmol,iso_g, & mol_index,iso_index,mol_niso,Ps,Tref,T,P,i,iso_mass, & nacc,nxtab,nk_func,nx,kindex,xtab,xacc,ktab,yvalues,pr, & nk,nu) implicit none include 'max.inc' include 'parameters.inc' include 'formats.inc' c c Purpose: this routine is called by slave processes. Its purpose is to generate c a spectral mesh (the "nu" array, containing "nk" wavenumber values) c between two given wavenumber values "numin" and "numax". c The list of molecules that must be taken into account is provided, and the routine c looks for spectral lines within the [numin-numax] interval in the appropriate c LBL data files. c c Input integer pindex double precision numin double precision numax double precision c2 integer sda,bslal double precision constant_dnu integer nlw_cz,np_lc,np_bl double precision Sref_frac integer Nmol integer mol_index(1:Nmol_max) integer iso_g(1:Nmol_max,1:Niso_max) integer iso_index(1:Nmol_max,1:Niso_max) integer mol_niso(1:Nmol_max) double precision Ps(1:Nmol_max) double precision Tref(1:Nmol_max) double precision T(1:Nmax) double precision P(1:Nmax) integer i double precision iso_mass(1:Nmol_max,1:Niso_max) integer nacc integer nxtab integer nk_func integer nx(1:Nytab_mx) integer kindex(1:Nytab_mx,1:2) double precision xtab(1:Nxtab_mx,1:3) double precision xacc(1:Nacc_mx,1:2) double precision ktab(1:Nktab_mx,1:3) double precision yvalues(1:Nytab_mx) double precision pr c Output integer nk double precision nu(1:Nkmx) c tmp integer nk_slimits integer Nlines,ios,j,Nlines_prev double precision nu0(1:Nlpnb_mx) integer index(1:Nlpnb_mx) integer isotope(1:Nlpnb_mx) double precision Sref(1:Nlpnb_mx) double precision gamma_air(1:Nlpnb_mx) double precision gamma_self(1:Nlpnb_mx) double precision Elow(1:Nlpnb_mx) double precision n_air(1:Nlpnb_mx) double precision n_self(1:Nlpnb_mx) double precision delta_air(1:Nlpnb_mx) c double precision numoy c tmp integer l,lf,nextl,prevl integer nfiles integer nl(1:nfiles_mx) character*(Nchar_mx) filename(1:nfiles_mx),file_database integer fformat(1:nfiles_mx) integer line,lmin,Nl_tmp integer molec,isotop double precision dnu(1:Nmax_tab) double precision ppres double precision nuc,nu1,nuc1,nuc2 double precision nu2(1:Nlpnb_mx) double precision gamma_L(1:Nline_mx) double precision gamma_D(1:Nline_mx) double precision dnu_min,dnuf double precision eps integer ifile,index_tmp,isotope_tmp,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 character*1 star_tmp double precision Sref_max,Sref_lim character*(Nchar_mx) results_file double precision slimits1(1:Nti_mx,2) integer slimits2(1:Nti_mx,2) integer sint,sintf,ik,jk double precision nuprec,Qref,Q,S_tmp(1:Nline_mx),ST integer nl_lf,ilf integer*8 date integer*8 day,month,year integer*8 hour,min,sec character*(Nchar_mx) logfile,stringp_file,string,datep_file double precision dist_minl,dist_minr integer linel,liner integer linelf,linerf,lvt double precision dist,gammal,gammar,gamma1,gamma2 double precision nu0l,nu0r double precision dfrac,gfrac,dnufrac,dnufrac1,dnufrac2 double precision pmin,dfact,gfact integer isof,ik0,nk0 double precision nu_lim1,nu_lim2,nbwt double precision nu_llimit,nu_hlimit integer line1,line2 logical line1_found,line2_found character*(Nchar_mx) lbl_file c Debug integer case,out c Debug integer strlen character*(Nchar_mx) label label='subroutine max_error_1cm_alrl' c Debug c string='./results/smesh.log' c call stringp_file_name(pindex,string,logfile) c string='./nlines' c call stringp_file_name(pindex,string,stringp_file) c call get_nlines(logfile,stringp_file,nl_lf) c open(14,file=logfile(1:strlen(logfile))) c do ilf=1,nl_lf c read(14,*) c enddo c string='date' c call stringp_file_name(pindex,string,datep_file) c call get_date(datep_file,date) c call convert_date(date,pindex,day,month,year,hour,min,sec) c write(14,*) 'IN: [',numin,'-',numax,'] on:', c & day,'/',month,'/',year,' @',hour,':',min,':',sec c close(14) c Debug nk=1 nu(1)=numin if (sda.eq.3) then do while (nu(nk).lt.numax) nk=nk+1 if (nk.gt.Nkmx) then call error(label) write(*,*) 'nk=',nk write(*,*) 'while Nkmx=',Nkmx stop endif nu(nk)=nu(nk-1)+constant_dnu enddo if (nu(nk).ne.numax) then nu(nk)=numax endif goto 666 endif c Debug c write(*,*) 'numin=',numin,' numax=',numax c Debug c --------------------------------------------------------------- c Retrieving spectral lines in the [numin-numax] interval c --------------------------------------------------------------- lbl_file='./optimizations/LBL_files.txt' open(12,file=lbl_file(1:strlen(lbl_file)) & ,status='old',iostat=ios) if (ios.ne.0) then ! list_file not found call error(label) write(*,*) 'File could not be found:' write(*,*) lbl_file(1:strlen(lbl_file)) stop else read(12,*) nfiles do ifile=1,nfiles read(12,37) fformat(ifile),nl(ifile),filename(ifile) enddo close(12) endif ! ios c Sref_max=0.0D+0 c string='./nlines' c call stringp_file_name(pindex,string,nlinesp_file) c First, we look for the maximum value of Sref accross all files: Sref_max c Debug c write(*,*) 'looking for Sref_max in',numin,'-',numax c Debug Nl_tmp=0 Nlines_prev=0 do ifile=1,nfiles c reading slimits1 and slimits2 for the current LBL database file call slimits_filename(ifile,results_file) open(15,file=results_file(1:strlen(results_file)) & ,status='old',iostat=ios) if (ios.ne.0) then ! list_file not found call error(label) write(*,*) 'file could not be found:' write(*,*) results_file(1:strlen(results_file)) stop else read(15,*) nk_slimits if (nk_slimits.gt.Nti_mx) then call error(label) write(*,*) 'nk_slimits=',nk_slimits write(*,*) 'in file: ', & results_file(1:strlen(results_file)) write(*,*) 'while Nti_mx=',Nti_mx stop else do ik=1,nk_slimits read(15,*) (slimits1(ik,j),j=1,2), & (slimits2(ik,j),j=1,2) enddo close(15) endif ! nk_slimits>Nti_mx endif ! ios=0 c identification of the lines that correspond to the given spectral interval sintf=0 do ik=1,nk_slimits if ((numin+1.0D-3.gt.slimits1(ik,1)).and. & (numax-1.0D-3.lt.slimits1(ik,2))) then sint=ik sintf=1 c Debug c write(*,*) 'sint=',sint c write(*,*) slimits1(ik,1),slimits1(ik,2) c Debug goto 331 endif enddo 331 continue if (sintf.eq.0) then call error(label) write(*,*) 'Could not identify interval:',numin,'-',numax write(*,*) 'in file:',results_file(1:strlen(results_file)) do ik=1,nk_slimits write(*,*) 'slimits1(',ik,')=', & slimits1(ik,1),slimits1(ik,2) enddo ! ik c ik=2027 c write(*,*) 'slimits1(',ik,',1)=',slimits1(ik,1), c & ' slimits1(',ik,',2)=',slimits1(ik,2) c write(*,*) 'numoy=',numoy stop endif c If there are no lines in this interval, for the current file, this file can be skipped if ((slimits2(sint,1).eq.-1).and.(slimits2(sint,2).eq.-1)) then goto 334 endif file_database=filename(ifile) Nlines_prev=slimits2(sint,2)-slimits2(sint,1)+1 if (Nlines_prev.gt.Nline_mx) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'called by process: ',pindex write(*,*) 'slimits file:' write(*,*) results_file(1:strlen(results_file)) write(*,*) 'File to read:' write(*,*) file_database(1:strlen(file_database)) write(*,*) 'Interval: ',numin,',',numax write(*,*) 'Nb. of lines in LBL file, for the interval:' write(*,*) Nlines_prev write(*,*) 'while Nline_mx=',Nline_mx stop endif c Debug c write(*,*) file_database c Debug open(10,file=file_database(1:strlen(file_database)) & ,status='old',iostat=ios) c Debug c write(*,*) file_database(1:strlen(file_database)) c Debug 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 endif if (slimits2(sint,1).gt.1) then do line=1,slimits2(sint,1)-1 read(10,*) enddo endif c Debug c if ((slimits2(sint,1).ne.-1).and.(slimits2(sint,2).ne.-1)) then c write(*,*) 'slimits file: ',results_file c & (1:strlen(results_file)) c write(*,*) 'File: ',file_database(1:strlen(file_database)) c Nlines_db=Nlines_db+slimits2(sint,2)-slimits2(sint,1)+1 c write(*,*) slimits2(sint,1),slimits2(sint,2), c & slimits2(sint,2)-slimits2(sint,1)+1,Nlines_db c endif c Debug do line=slimits2(sint,1),slimits2(sint,2) 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¯¹] 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¯¹] else if (fformat(ifile).eq.3) then ! CDSD 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¯¹] endif c Debug c write(*,*) 'nu0_tmp(',line,')=',nu0_tmp,isotope_tmp c Debug call indexes(Nmol,mol_index,iso_index,mol_niso, & index_tmp,isotope_tmp,molec,isotop,isof) c Debug c write(*,*) 'isof(',line,')=',isof if (isof.eq.0) then goto 336 endif c Debug call TIPS(index_tmp,Tref(index_tmp), & isotope_tmp,iso_g(molec,isotop),Qref) call TIPS(index_tmp,T(i),isotope_tmp, & iso_g(molec,isotop),Q) c Debug c write(*,*) 'Q(',line,')=',Q c Debug c for some reason, the BD_TIPS_2003 routine can not compute c partition function Q for species index 34 (Q=0. is imposed). if (Q.ne.0.0D+0) then Nl_tmp=Nl_tmp+1 if (Nl_tmp.gt.Nline_mx) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'Nl_tmp=',Nl_tmp write(*,*) 'exceeds Nline_mx=',Nline_mx stop endif call line_intensity(Tref(index_tmp),T(i),nu0_tmp, & Sref_tmp,Elow_tmp,c2,Q,Qref,S_tmp(Nl_tmp)) endif c if (S_tmp.gt.Sref_max) then c Sref_max=S_tmp c endif c if (Sref_tmp.gt.Sref_max) then c Sref_max=Sref_tmp c endif 336 continue enddo ! line close(10) c Debug c write(*,*) 'lines read:',slimits2(sint,2)-slimits2(sint,1)+1 c Debug 334 continue enddo ! loop over files c Now the value of "Sref_max" is known c Debug c write(*,*) numin,numax,' Nl_tmp=',Nl_tmp c open(10,file='./S_tmp.txt') c do line=1,Nl_tmp c write(10,*) S_tmp(line) c enddo c close(10) c Debug c Sref_lim is the lower limit in line sorting call Slimit(Nl_tmp,S_tmp,pr,Sref_lim) c Sref_lim=Sref_max/Sref_frac c Debug c write(*,*) 'Sref_max is found' c write(*,*) 'numin=',numin c write(*,*) 'numax=',numax c write(*,*) 'Sref_max=',Sref_max c write(*,*) 'Sref_lim=',Sref_lim c write(*,*) 'Sref_frac=',Sref_frac c Debug c Debug c write(*,*) 'looking for lines in the interval' c Debug c Now we have to read each file again, and read line data that correspond to Sref>=Sref_lim Nlines=0 do ifile=1,nfiles c reading slimits1 and slimits2 for the current LBL database file call slimits_filename(ifile,results_file) open(15,file=results_file(1:strlen(results_file)) & ,status='old',iostat=ios) if (ios.ne.0) then ! list_file not found call error(label) write(*,*) 'file could not be found:' write(*,*) results_file(1:strlen(results_file)) stop else read(15,*) nk_slimits if (nk_slimits.gt.Nti_mx) then call error(label) write(*,*) 'nk_slimits=',nk_slimits write(*,*) 'in file: ', & results_file(1:strlen(results_file)) write(*,*) 'while Nti_mx=',Nti_mx stop else do ik=1,nk_slimits read(15,*) (slimits1(ik,j),j=1,2), & (slimits2(ik,j),j=1,2) enddo close(15) endif ! nk_slimits>Nti_mx endif ! ios=0 c identification of the lines that correspond to the given spectral interval sintf=0 do ik=1,nk_slimits if ((numin+1.0D-3.gt.slimits1(ik,1)).and. & (numax-1.0D-3.lt.slimits1(ik,2))) then sint=ik sintf=1 goto 332 endif enddo 332 continue if (sintf.eq.0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'Could not identify interval:',numin,'-',numax write(*,*) 'in file:',results_file(1:strlen(results_file)) c ik=2027 c write(*,*) 'slimits1(',ik,',1)=',slimits1(ik,1), c & ' slimits1(',ik,',2)=',slimits1(ik,2) c write(*,*) 'numoy=',numoy stop endif c If there are no lines in this interval, for the current file, this file can be skipped if ((slimits2(sint,1).eq.-1).and.(slimits2(sint,2).eq.-1)) then goto 335 endif file_database=filename(ifile) open(10,file=file_database(1:strlen(file_database)) & ,status='old',iostat=ios) c Debug c write(*,*) file_database(1:strlen(file_database)) c Debug 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 endif if (slimits2(sint,1).gt.1) then do line=1,slimits2(sint,1)-1 read(10,*) enddo endif do line=slimits2(sint,1),slimits2(sint,2) 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] & ,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 & ,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] endif c Take only lines for which Sref>=Sref_lim call indexes(Nmol,mol_index,iso_index,mol_niso, & index_tmp,isotope_tmp,molec,isotop,isof) ppres=Ps(molec) c Debug c if (Sref_tmp.eq.Sref_max) then c write(*,*) numin,numax,Sref_tmp,nu0_tmp, c & nu0_tmp+delta_air_tmp*P(i) c endif c write(*,*) 'Sref_tmp=',Sref_tmp c Debug c computation of S(T) c Debug if (isof.eq.0) then goto 338 endif c Debug call TIPS(index_tmp,Tref(index_tmp),isotope_tmp, & iso_g(molec,isotop),Qref) call TIPS(index_tmp,T(i),isotope_tmp, & iso_g(molec,isotop),Q) c for some reason, the BD_TIPS_2003 routine can not compute c partition function Q for species index 34 (Q=0. is imposed). if (Q.ne.0.0D+0) then call line_intensity(Tref(index_tmp),T(i),nu0_tmp, & Sref_tmp,Elow_tmp,c2,Q,Qref,ST) endif c if ((ST.ge.Sref_lim) & .and.(ppres.ne.0.0D+0) & .and.(isof.eq.1)) then Nlines=Nlines+1 if (Nlines.gt.Nlpnb_mx) then write(*,*) 'Error in routine max_error_1cm_alrl :' write(*,*) 'Nlpnb_mx=',Nlpnb_mx,' was reached' stop endif c Debug c write(*,*) 'retained, Nlines=',Nlines c Debug 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 c Debug Elow(Nlines)=Elow_tmp n_air(Nlines)=n_air_tmp delta_air(Nlines)=delta_air_tmp n_self(Nlines)=n_self_tmp c Debug c else c write(*,*) 'rejected' c Debug endif 338 continue enddo ! line close(10) c Debug c write(*,*) 'closing:',file_database(1:strlen(file_database)) c Debug 335 continue enddo ! loop over files c Debug c write(*,*) 'lines identified: Nlines=',Nlines,' Now sorting' c write(*,*) numin,numax,' Nlines=',Nlines c stop c Debug c This is (probably ?) a bug from the HITEMP database: c some values of the gamma_self are null ! Which is a problem because c in the case the gamma_self is null for a given line, the Lorentz c width will be null for this line, the main consequence is that the spectral c discretization will fail for this line (and crash the program). c When a null gamma_self is detected, I set its value to the c first non-null gamma_self encountered. do line=1,Nlines if (gamma_self(line).eq.0.0D+0) then if (line.lt.Nlines/2) then lf=0 do l=line+1,Nlines if (gamma_self(l).ne.0.0D+0) then lf=1 nextl=l goto 127 endif enddo 127 continue if (lf.eq.0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nextl could not be found' write(*,*) 'Nlines=',Nlines write(*,*) 'line=',line write(*,*) 'Are you reading HITEMP data file ?' stop else gamma_self(line)=gamma_self(nextl) endif else lf=0 do l=line-1,1,-1 if (gamma_self(l).ne.0.0D+0) then lf=1 prevl=l goto 128 endif enddo 128 continue if (lf.eq.0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'prevl could not be found' write(*,*) 'Nlines=',Nlines write(*,*) 'line=',line write(*,*) 'Are you reading HITEMP data file ?' write(*,*) 'gamma_self(',line,')=',gamma_self(line) write(*,*) 'nu0(',line,')=',nu0(line) stop else gamma_self(line)=gamma_self(prevl) endif endif endif enddo c Debug do line=1,Nlines if (gamma_self(line).eq.0.0D+0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'gamma_self(',line,')=',gamma_self(line) write(*,*) 'error 0' stop endif enddo c Debug c Sorting data (not really necessary, but who cares) c Debug c write(*,*) 'calling sort_all_data2' c Debug call sort_all_data2(Nlines,nu0,index,isotope,Sref, & gamma_air,gamma_self,Elow,n_air,delta_air,n_self) c Debug c write(*,*) 'sorting done. Now discretizing interval' c Debug c --------------------------------------------------------------- c Spectral lines have been identified. c Now discretizing the [numin-numax] interval. c --------------------------------------------------------------- if (sda.eq.4) then line1_found=.false. do line=1,Nlines if ((.not.line1_found).and.(nu0(line).gt.numin)) then line1_found=.true. line1=line ! is the first line after 'numin' goto 661 endif 661 continue if (.not.line1_found) then call error(label) write(*,*) 'line 1 was not found' write(*,*) 'numin=',numin write(*,*) 'nu0(1)=',nu0(1) write(*,*) 'nu0(',Nlines,')=',nu0(Nlines) stop endif enddo ! line line2_found=.false. do line=1,Nlines if ((.not.line2_found).and.(nu0(line).gt.numax)) then line2_found=.true. line2=line ! is the first line after 'numax' goto 662 endif 662 continue if (line2_found) then line2=line2-1 ! use last line before 'numax' else if (nu0(Nlines).lt.numax) then line2_found=.true. line2=Nlines ! use last line endif endif if (.not.line2_found) then call error(label) write(*,*) 'line 2 was not found' write(*,*) 'numax=',numax write(*,*) 'nu0(1)=',nu0(1) write(*,*) 'nu0(',Nlines,')=',nu0(Nlines) stop endif enddo ! line c take spectral discretization as line centers do line=line1,line2 nk=nk+1 if (nk.gt.Nkmx) then call error(label) write(*,*) 'nk=',nk write(*,*) 'while Nkmx=',Nkmx stop endif nu(nk)=nu0(line) enddo ! line nk=nk+1 if (nk.gt.Nkmx) then call error(label) write(*,*) 'nk=',nk write(*,*) 'while Nkmx=',Nkmx stop endif nu(nk)=numax goto 666 endif if (Nlines.eq.0) then call choose_nbw(nu(nk),nbwt) if (sda.eq.1) then c nk0=10*int(nbwt) nk0=int(nbwt) else c nk0=int(nbwt) nk0=int(nbwt)/2 endif c write(*,*) c write(*,*) 'WARNING from ',label(1:strlen(label)),' :' c write(*,*) 'No line found in interval:',numin,'-',numax c write(*,*) 'number of points added:',nk0 dnuf=(numax-numin)/(nk0-1) c uniform discretization between numin and numax using nk0 points call fast_discretization(numin,numax,dnuf,nk,nu) goto 666 endif c Debug c write(*,*) 'Nlines=',Nlines c Debug do line=1,Nlines nuc=nu0(line)+delta_air(line)*P(i) c Debug c write(*,*) 'nuc(',line,')=',nuc c Debug call indexes(Nmol,mol_index,iso_index,mol_niso, & index(line),isotope(line),molec,isotop,isof) c Debug if (isof.eq.0) then write(*,*) 'error from routine max_error_1cm_alrl:' write(*,*) 'isof=',isof stop endif c Debug call Lorentz_linewidth(Tref(index(line)),T(i),P(i), & Ps(molec),n_air(line),n_self(line), & gamma_air(line),gamma_self(line),gamma_L(line)) c Debug if (gamma_L(line).eq.0.0D+0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'gamma_L(',line,')=',gamma_L(line) write(*,*) 'Tref(',index(line),')=',Tref(index(line)) write(*,*) 'T(',i,')=',T(i) write(*,*) 'P(',i,')=',P(i) write(*,*) 'Ps(',molec,')=',Ps(molec) write(*,*) 'n_air(',line,')=',n_air(line) write(*,*) 'n_self(',line,')=',n_self(line) write(*,*) 'gamma_air(',line,')=',gamma_air(line) write(*,*) 'gamma_self(',line,')=',gamma_self(line) stop endif c Debug call Doppler_linewidth(T(i),iso_mass(molec,isotop), & nu0(line),gamma_D(line)) enddo c "gfrac" is used to define line's central zone: if "nu0" is the line central wavenumber, c the central zone for this line will be considered as the [nu0-gfrac*gamma,nu0+gfrac*gamma] c with "gamma" the line width. gfrac=dble(nlw_cz)/2.0D+0 c "dnufrac1" is the number of wavenumber points that have to be taken within c the line central zone. dnufrac1=dble(np_lc)/2.0D+0 ! line centers c "dnufrac2" is the number of wavenumber points that have to be taken between c two consecutive line central zones. dnufrac2=dble(np_bl) ! between lines if (sda.eq.2) then c ************************* c This is when degraded resolution is required c ************************* c discretization of the interval [nu(nk)-nuc1] call param1(nu0(1),P(i),delta_air(1),gamma_L(1),gamma_D(1), & nuc1,gamma1) c Debug c write(*,*) 'nuc1=',nuc1 c Debug nu_llimit=nuc1-gfrac*gamma1 ! lower limit of the line's central zone c At this stage, nk=1 and nu(1)=numin if (nu(nk).le.nu_llimit) then nu_lim1=nu(nk) ! numin nu_lim2=nu_llimit ! lower limit of the first line's central zone dnuf=(nu_lim2-nu_lim1)/(dnufrac2/2.0D+0) ! spectral resolution c uniform discretization between nu_min and the lower limit of the first line's central zone call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c This routine will update the "nu" array and "nk", starting from nk=1 and nu(1)=numin endif nu_lim1=nu(nk) ! lower limit of the first line's central zone nu_lim2=nuc1 ! first line's center dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c uniform discretization between the beginning of the first line's central zone c and the first line's center call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c "nk" and "nu" have been updated accordingly c In the case there is at most one line within the [numin-numax] interval, the c routine can jump to the last step: the discretization of the spectral interval c between the line center and numax if (Nlines.le.1) then goto 436 endif c Otherwise, it has to discretize the interval between the first line's center c and the last line's center do line=1,Nlines-1 call param1(nu0(line),P(i),delta_air(line),gamma_L(line), & gamma_D(line),nuc1,gamma1) c "nuc1" and "gamma1" are the line center and line width of the line index "line" call param1(nu0(line+1),P(i),delta_air(line+1), & gamma_L(line+1),gamma_D(line+1),nuc2,gamma2) c "nuc2" and "gamma2" are the line center and line width of the line index "line+1" c Debug c write(*,*) 'nuc1=',nuc1,' nuc2=',nuc2 c Debug c dist=nu0(line+1)-nu0(line) dist=nuc2-nuc1 c "dist" is the spectral distance between these two lines centers. if ((gamma1+gamma2)*gfrac.ge.dist) then c In the case lines central zones do overlap c "nu(nk)" is always the center of the previous line's center nu_lim1=nu(nk) ! line 1 center wavenumber nu_lim2=nuc2 ! line 2 center wavenumber dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c uniform discretization between two consecutive line centers call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c Debug c write(*,*) 'no overlap, nu(',nk,')=',nu(nk) c Debug else c In the case line central zones do not overlap (lines are well separated) nu_lim1=nu(nk) ! line's center wavenumber nu_lim2=nuc1+gamma1*gfrac ! upper limit of the first line's central zone dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c uniform discretization between the center of line index "line" and c the upper limit of the first line's central zone call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) nu_lim1=nu(nk) ! upper limit of the first line's central zone nu_lim2=nuc2-gamma2*gfrac ! lower limit of the second line's central zone dnuf=(nu_lim2-nu_lim1)/dnufrac2 ! spectral resolution c uniform discretization between the upper limit of the first line's central zone c and the lower limit of the second line's central zone call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) nu_lim1=nu(nk) ! lower limit of the second line's central zone c nu_lim2=nu0(line+1) ! second line's center wavenumber nu_lim2=nuc2 ! second line's center wavenumber dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c uniform discretization between the lower limit of the second line's central zone c and the second line's center call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c Debug c write(*,*) 'overlap, nu(',nk,')=',nu(nk) c Debug endif enddo ! do it for every line except the last one 436 continue if (nu(nk).gt.numax) then nu(nk)=numax goto 339 endif call param1(nu0(Nlines),P(i),delta_air(Nlines), & gamma_L(Nlines),gamma_D(Nlines),nuc2,gamma2) c "nuc2" and "gamma2" are the center wavenumber and the line width c of the last line within the [numin,numax] spectral interval nu_hlimit=nuc2+gamma2*gfrac ! higher limit of the last line's central zone c Debug c write(*,*) 'nuc2=',nuc2 c write(*,*) 'nu_hlimit=',nu_hlimit c do i=1,nk c write(*,*) 'nu(',i,')=',nu(i) c enddo if (nu(nk).gt.numax) then write(*,*) 'Error:' write(*,*) 'nu(',nk,')=',nu(nk) write(*,*) 'numax=',numax endif c Debug if (nu(nk)+gamma2*gfrac.le.numax) then c if the last line's central zone is within the [numin,numax] spectral interval nu_lim1=nu(nk) ! last line's center nu_lim2=nu_hlimit ! higher limit of the last line's central zone dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c uniform discretization between the last line's center c and the higher limit of the last line's central zone call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) nu_lim1=nu(nk) ! higher limit of the last line's central zone nu_lim2=numax ! end of the [numin,numax] spectral interval dnuf=(nu_lim2-nu_lim1)/(dnufrac2/2.0D+0) ! spectral resolution c uniform discretization between the higher limit of the last line's central zone c and numax call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c Debug if (nu(nk).ne.numax) then write(*,*) 'error 10:' write(*,*) 'nu(',nk,')=',nu(nk),' numax=',numax write(*,*) 'nu_lim1=',nu_lim1,' nu_lim2=',nu_lim2 write(*,*) 'dnuf=',dnuf stop endif c Debug else c if the upper limit of the last line's central zone is beyond "numax" nu_lim1=nu(nk) ! last line's center nu_lim2=numax ! end of the [numin,numax] spectral interval dnuf=(nu_lim2-nu_lim1)/dnufrac1 ! spectral resolution c Debug if (nu_lim2.le.nu_lim1) then write(*,*) 'error:' write(*,*) 'nu_lim1=',nu_lim1 write(*,*) 'nu_lim2=',nu_lim2 endif c Debug c uniform discretization between the last line's center and "numax" call fast_discretization(nu_lim1,nu_lim2,dnuf,nk,nu) c Debug if (nu(nk).ne.numax) then write(*,*) 'error 11:' write(*,*) 'nu(',nk,')=',nu(nk),' numax=',numax write(*,*) 'nu_lim1=',nu_lim1,' nu_lim2=',nu_lim2 write(*,*) 'dnuf=',dnuf write(*,*) 'npt=',int((nu_lim2-nu_lim1)/dnuf)+1 stop endif c Debug endif c Debug 116 continue do ik=2,nk if (nu(ik).le.nu(ik-1)) then c remove nu(ik) do jk=ik,nk-1 nu(jk)=nu(jk+1) enddo ! jk nk=nk-1 goto 116 endif enddo c write(*,*) 'nu(',nk,')=',nu(nk),' numax=',numax c stop c Debug c do not execute what is after then "endif": go immediately to c the end of the routine goto 666 endif ! sda=2 339 continue c ************************* c This is when reference discretization is required c ************************* c nk=1 and nu(1)=numin because when the routine starts from there, c it's because the reference discretization algorithm has been required do while (nu(nk).lt.numax) ! discretization of the [numin-numax] interval nu1=nu(nk) if ((sda.eq.1).and.(bslal.eq.0)) then ! disable between-line spectral meshing algorithm lvt=1 c reference discretization algorithm with no special feature goto 551 endif lvt=1 ! standard c --------------------------------------------------------------------- c Special modification for points between two very close lines: c the spectral discretization may not be correct using standard method c c 1 - Identify lines at the left and at the right of current nu1 linelf=0 do line=1,Nlines if ((nu1-nu0(line)).gt.0.0D+0) then linelf=1 linel=line dist_minl=nu1-nu0(line) goto 441 endif enddo 441 continue c Debug c if ((nu1.gt.4215.98).and.(nu1.lt.4216.03)) then c write(*,*) 'nu1=',nu1,' linelf=',linelf,' linel=',linel c endif c Debug if (linelf.eq.1) then do line=1,Nlines if (((nu1-nu0(line)).gt.0.0D+0).and. & ((nu1-nu0(line)).lt.dist_minl)) then linel=line dist_minl=nu1-nu0(line) endif enddo endif linerf=0 do line=1,Nlines if ((nu0(line)-nu1).gt.0.0D+0) then linerf=1 liner=line dist_minr=nu0(line)-nu1 goto 442 endif enddo 442 continue c Debug c if ((nu1.gt.4215.98).and.(nu1.lt.4216.03)) then c write(*,*) 'nu1=',nu1,' linerf=',linerf,' liner=',liner c endif c Debug if (linerf.eq.1) then do line=1,Nlines if (((nu0(line)-nu1).gt.0.0D+0).and. & ((nu0(line)-nu1).lt.dist_minr)) then liner=line dist_minr=nu0(line)-nu1 endif enddo endif c 2 - compute distance between nearest lines if (linelf.eq.1) then nu0l=nu0(linel) else nu0l=numin endif if (linerf.eq.1) then nu0r=nu0(liner) else nu0r=numax endif dist=nu0r-nu0l c Debug c if ((nu1.gt.4215.98).and.(nu1.lt.4216.03)) then c write(*,*) 'nu1=',nu1,' dist=',dist c endif c Debug c 3 - if this distance is too short, c and nu1 is in the middle of the two nearest lines, c dnu_min is very likely wrong if (linelf.eq.1) then gammal=gamma_L(linel) if (gamma_D(linel).gt.gammal) then gammal=gamma_D(linel) endif else gammal=0.0D+0 endif if (linerf.eq.1) then gammar=gamma_L(liner) if (gamma_D(liner).gt.gammar) then gammar=gamma_D(liner) endif else gammar=0.0D+0 endif c Debug c if ((nu1.gt.4215.98).and.(nu1.lt.4216.03)) then c write(*,*) 'nu1=',nu1,' gammal=',gammal,' gammar=',gammar c endif c Debug c set variables dfrac,gfrac,dnufrac call choice_bslal(bslal,dfrac,gfrac,dnufrac,pmin,dfact,gfact) c In the case the reference spectral discretization algorithm is required, c and special algorithm is required between strong lines, c and frequency is between strong close lines if ( & (sda.eq.1).and.(bslal.gt.0).and. & (P(i).ge.pmin).and. & (dist.lt.dfrac*(gammal+gammar)).and. & (nu1.ge.nu0l+gammal*gfrac).and. & (nu1.le.nu0r-gammar*0*gfrac) & ) then c Special case of close lines with (potentially) strong overlap -> standard algorithm c based on Lorentz and Voigt functions tabulation may be incorrect dnu_min=(gammal+gammar)/dnufrac ! should be fine lvt=0 case=1 c In the case the degraded spectral discretization algorithm is required else if ( & (sda.eq.2) & ) then c gfrac=1.0D+0 c dnufrac1=3.0D+0 c dnufrac2=1.5D+0 c in the case nu1 is close to left line center if (dabs(nu1-nu0l).lt.gfrac*gammal) then dnu_min=gammal/dnufrac1 c in the case nu1 is close to right line center else if (dabs(nu1-nu0r).lt.gfrac*gammar) then dnu_min=gammar/dnufrac1 c in all other cases, nu1 is between left and right lines centers else dnu_min=(gammal+gammar)/dnufrac2 endif lvt=0 case=3 endif c for every spectral discretization algorithm, in the case c pressure is too low if ( & (P(i).lt.pmin).and. & (dist.gt.dfact*(gammal+gammar)).and. & ((nu1.ge.nu0l+gammal*gfact).and. & (nu1.le.nu0r-gammar*0*gfact)).and. & (nu0r-gammar*gfact-nu1.gt.1.0D-12) & ) then c Special case of very low pressure levels, with very distant strong lines -> nothing c will have to be computed until the beginning of the next strong spectral line dnu_min=nu0r-gammar*gfact-nu1 c Debug c write(*,*) nu1,dnu_min,nu1+dnu_min,nu0l,nu0r,gammal,gammar c Debug lvt=0 case=2 endif c Debug c if ((nu1.gt.4215.98).and.(nu1.lt.4216.03)) then c write(*,*) 'linelf=',linelf,' linerf=',linerf c write(*,*) 'linel=',linel,' nu0=',nu0l, c & ' gammal=',gammal c write(*,*) 'liner=',linel,' nu0=',nu0r, c & ' gammar=',gammar c write(*,*) 'dist=',dist,dfrac*(gammal+gammar) c write(*,*) 'nu1=',nu1,nu0l+gammal*gfrac, c & nu0r-gammar*0*gfrac c endif c Debug c --------------------------------------------------------------------- 551 continue if (lvt.eq.1) then if (nk.gt.1) then nuprec=nu(nk-1) else nuprec=nu(nk) endif do line=1,Nlines c Debug if (gamma_L(line).eq.0.0D+0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'gamma_L(',line,')=',gamma_L(line) stop endif c Debug call find_nu2_tab(pindex,nuprec,nu1,nu0(line), & delta_air(line),P(i),gamma_L(line),gamma_D(line), & xacc,nacc,xtab,nxtab,ktab,nk_func,nx,yvalues,kindex, & nu2(line)) dnu(line)=dabs(nu2(line)-nu1) c Debug call test_nan(dnu(line),out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'error 2' write(*,*) 'dnu(',line,')=',dnu(line) write(*,*) 'nu2(',line,')=',nu2(line) write(*,*) 'nu1=',nu1 stop endif c Debug enddo call find_min(dnu,Nlines,dnu_min,lmin) endif if (dnu_min.lt.0.0D+0) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'dnu_min=',dnu_min stop endif nk=nk+1 if (nk.gt.Nkmx) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'Nkmx has been exceeded' stop endif nu(nk)=nu1+dnu_min c Debug call test_nan(nu(nk),out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'error 1' write(*,*) 'nu(',nk,')=',nu(nk) write(*,*) 'nu1=',nu1 write(*,*) 'dnu_min=',dnu_min stop endif c Debug eps=dabs(dabs(nu(nk)/nu(nk-1))-1.0D+0) if (eps.lt.1.0D-12) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu(',nk-1,')=',nu(nk-1) write(*,*) 'nu(',nk,')=',nu(nk) if (lvt.eq.0) then write(*,*) 'lvt=',lvt,' case=',case write(*,*) 'linelf=',linelf,' linerf=',linerf write(*,*) 'nu0l=',nu0l,' nu0r=',nu0r write(*,*) 'gammal=',gammal,' gammar=',gammar, & ' dist=',dist write(*,*) 'nu1=',nu1,' dnu_min=',dnu_min if (case.eq.2) then write(*,*) 'P(',i,')=',P(i) write(*,*) 'nu0r=',nu0r write(*,*) 'gammar*gfact=',gammar*gfact write(*,*) 'nu0r-gammar*gfact=',nu0r-gammar*gfact write(*,*) 'dnu_min=',nu0r-gammar*gfact-nu1 endif endif stop endif c Debug call test_nan(nu(nk),out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'error 0' write(*,*) 'nu(',nk,')=',nu(nk) stop endif c Debug enddo ! while (nu(nk).lt.numax) 666 continue c Debug c string='./nlines' c call stringp_file_name(pindex,string,stringp_file) c call get_nlines(logfile,stringp_file,nl_lf) c open(14,file=logfile(1:strlen(logfile))) c do ilf=1,nl_lf c read(14,*) c enddo c string='./date' c call stringp_file_name(pindex,string,datep_file) c call get_date(datep_file,date) c call convert_date(date,pindex,day,month,year,hour,min,sec) c write(14,*) 'OUT: [',numin,'-',numax,'] on:', c & day,'/',month,'/',year,' @',hour,':',min,':',sec c close(14) c Debug return end subroutine calc_k_interp(k1,k2,nu1,nu2,nu,k_interp) implicit none double precision k1,k2,nu1,nu2,nu double precision k_interp k_interp=k1+(k2-k1)/(nu2-nu1)*(nu-nu1) return end subroutine find_nu2_tab(pindex,nuprec,nu1,nu0,delta,P, & gamma_L,gamma_D,xacc,nacc,xtab,nxtab,ktab,nk_func, & nx,yvalues,kindex,nu2) implicit none c Purpose: to compute frequency nu2 given frequency nu1, Lorentz line width gamma_L c and Doppler line width gamma_D c c Inputs: c + pindex: index of the process that calls this routine c + nuprec: previous value of nu in the nu grid c + nu1: starting vavenumber [cm¯¹] c + nu0: vacuum wavenumber [cm¯¹] c + delta: air pressure-induced line shift [cm¯¹.atm¯¹] c + P: the total pressure [atm] c + gamma_L: Lorentz line half-width [cm¯¹] c + gamma_D: Doppler line half-width [cm¯¹] c + xacc: acceleration array to retrieve information in "xtab" faster c + nacc: dimension of array "xacc" c + xtab: tabulation array of Lorentz profile c + nxtab: dimension of array xtab c + ktab: values of [x1,x1p,x2] for each k function tabulated c + nk_func: number of k functions tabulated c + nx: number of [x1,x1p,x2] triplets for each k function tabulated c + yvalues: value of y for every k function tabulated c + kindex: c - kindex(i,1) is the index (in ktab) that locates the beginning of [x1,x1p,x2] sequence for yvalue 'i' c - kindex(i,2) is the index (in ktab) that locates the end of [x1,x1p,x2] sequence for yvalue 'i' c c Output: c + nu2: ending wavenumber [cm¯¹]; this is the wavenumber at which c the maximum relative error has been reached c include 'max.inc' integer pindex double precision nuprec double precision nu1,nu0,delta,P,gamma_L,gamma_D,nu2,nu2d,nu2l double precision nuc,x1,x2,x2d,y1,y2,y double precision xtab(1:Nxtab_mx,1:3) double precision xacc(1:Nacc_mx,1:2) integer nxtab,nacc,i,j,x2f,i0f,i0,x2df double precision a,b,sign,eps_max c integer nk_func,nx(1:Nytab_mx) integer kindex(1:Nytab_mx,1:2) double precision ktab(1:Nktab_mx,1:3) double precision yvalues(1:Nytab_mx) c integer kfunc,kfunc0,kfunc0f,kfunc1,kfunc2 integer nx1,nx2,out double precision ktab1(1:Nktab_mx,1:3) double precision ktab2(1:Nktab_mx,1:3) double precision ktab2p(1:Nktab_mx,1:3) double precision ktab_yi(1:Nktab_mx,1:3) double precision kmax1,kmax2,kmax_min integer doppler,lorentz integer i1,i2,i2f,k1,k2,k,x2dts double precision dnu_min,dnu_fail parameter(dnu_min=5.0D-7) ! cm¯¹ c Debug double precision eps c Debug integer strlen character*(Nchar_mx) label label='subroutine find_nu2_tab' c Debug c write(*,*) 'This is find_nu2_tab' c write(*,*) 'nu1=',nu1 c write(*,*) 'gamma_L=',gamma_L,' gamma_D=',gamma_D c Debug eps_max=1.0D-2 ! 1% dnu_fail=nu1-nuprec c Debug if (dnu_fail.lt.0.0D+0) then write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'nuprec=',nuprec write(*,*) 'nu1=',nu1 write(*,*) endif c Debug if (dnu_fail.eq.0.0D+0) then dnu_fail=dnu_min endif c Test do kfunc=1,nk_func k1=kindex(kfunc,1) k2=kindex(kfunc,2) c write(*,*) k2-k1+1,nx(kfunc) do k=k1+1,k2 if (ktab(k-1,2).ne.ktab(k,1)) then write(*,*) 'ktab(',k-1,',2)=',ktab(k-1,2) write(*,*) 'ktab(',k,',1)=',ktab(k,1) stop endif enddo ! k enddo ! kfunc c Test nuc=nu0+delta*P doppler=0 lorentz=0 x2dts=1 c First, try to locate (x,y) in the ktab array x1=(nu1-nuc)/gamma_D*dsqrt(dlog(2.0D+0)) y=gamma_L/gamma_D*dsqrt(dlog(2.0D+0)) c Debug c write(*,*) 'x1=',x1,' y=',y c write(*,*) 'yvalues(1)=',yvalues(1), c & ' yvalues(',nk_func,')=',yvalues(nk_func) c Debug kfunc0f=0 do kfunc=2,nk_func if ((yvalues(kfunc-1).le.y).and.(yvalues(kfunc).ge.y)) then kfunc0=kfunc kfunc0f=1 c Debug c write(*,*) 'kfunc0=',kfunc0 c write(*,*) yvalues(kfunc0-1),y,yvalues(kfunc0) c Debug goto 321 endif enddo ! i 321 continue if (kfunc0f.eq.0) then c write(*,*) 'Error from routine find_nu2_tab:' c write(*,*) 'y=',y,' was not found among yvalues' c write(*,*) 'y min=',yvalues(1) c write(*,*) 'y max=',yvalues(nk_func) c write(*,*) if (y.lt.yvalues(1)) then c if y is lower than the first value of y that is tabulated, then interpolation will c be performed between the first and second values of y interpolated kfunc0=2 kfunc0f=1 c Debug c write(*,*) 'y=',y,' < yvalues(1)=',yvalues(1) c Debug else if (y.gt.yvalues(nk_func)) then c if y is greater that the last value of y that is tabulated, then switch to lorentz mode: c use the interpolation over the Lorentz line profile doppler=0 lorentz=1 c Debug c write(*,*) 'y=',y,' > yvalues(',nk_func,')=', c & yvalues(nk_func) c Debug goto 333 else write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'y=',y write(*,*) endif ! value of y endif ! kfunc0f.eq.0 c kmax1=ktab(kindex(kfunc0-1,2),2) ! limit value of the x-range covered by k-function index kfunc0-1 c kmax2=ktab(kindex(kfunc0,2),2) ! limit value of the x-range covered by k-function index kfunc0 kmax1=ktab(kindex(kfunc0-1,2),2) ! limit value of the x-range covered by k-function index kfunc0-1 kmax2=ktab(kindex(kfunc0,2),2) ! limit value of the x-range covered by k-function index kfunc0 c Debug c write(*,*) 'kmax1=',kmax1,' kmax2=',kmax2 c Debug kmax_min=kmax1 ! kmax_min is the minium value between kmax1 and kmax2 kfunc1=kfunc0-1 kfunc2=kfunc0 nx1=nx(kfunc0-1) nx2=nx(kfunc0) y1=yvalues(kfunc0-1) y2=yvalues(kfunc0) if (kmax2.lt.kmax_min) then kmax_min=kmax2 kfunc1=kfunc0 kfunc2=kfunc0-1 nx1=nx(kfunc0) nx2=nx(kfunc0-1) y1=yvalues(kfunc0) y2=yvalues(kfunc0-1) endif c Debug c write(*,*) 'kmax_min=',kmax_min c Debug if (dabs(x1).le.kmax_min) then doppler=1 lorentz=0 else doppler=0 lorentz=1 endif 333 continue c Debug c write(*,*) 'doppler=',doppler,' lorentz=',lorentz c Debug if ((doppler.eq.1).and.(lorentz.eq.0)) then c Debug c write(*,*) 'kfunc0=',kfunc0,' y=',y c Debug c In Doppler mode x2df=1 do i=1,nx1 do j=1,3 ktab1(i,j)=ktab(kindex(kfunc1,1)+i-1,j) enddo ! j enddo ! i c Test k1=kindex(kfunc1,1) k2=kindex(kfunc1,2) c continuity test do k=k1+1,k2 if (ktab1(k-1,2).ne.ktab1(k,1)) then write(*,*) 'ktab1(',k-1,',2)=',ktab1(k-1,2) write(*,*) 'ktab1(',k,',1)=',ktab1(k,1) do i=k1,k2 write(*,*) 'ktab1(',i,')=',(ktab1(i,j),j=1,3) enddo c stop x2df=0 goto 521 endif enddo ! k c Test do i=1,nx2 do j=1,3 ktab2(i,j)=ktab(kindex(kfunc2,1)+i-1,j) enddo ! j c continuity test if (i.gt.1) then if ((ktab2(i-1,2).ne.ktab2(i,1)).and.(i.ne.1)) then write(*,*) 'ktab2(',i-1,',2)=',ktab2(i-1,2) write(*,*) 'ktab2(',i,',1)=',ktab2(i,1) c stop x2df=0 goto 521 endif endif c Debug c write(*,*) 'ktab2(',i,')=',(ktab2(i,j),j=1,3) c Debug enddo ! i c Debug c write(*,*) 'nx1=',nx1 c write(*,*) 'ktab1(1,1)=',ktab1(1,1), c & ' ktab1(',nx1,',2)=',ktab1(nx1,2) c write(*,*) 'nx2=',nx2 c write(*,*) 'ktab2(1,1)=',ktab2(1,1), c & ' ktab2(',nx2,',2)=',ktab2(nx2,2) c Debug do i1=1,nx1 i2f=0 c do i2=2,nx2 c if ((ktab2(i2-1,1).le.ktab1(i1,1)).and. c & (ktab2(i2,1).ge.ktab1(i1,1))) then c i2f=1 c goto 322 c endif c enddo ! i2 do i2=1,nx2 if ((ktab2(i2,1).le.ktab1(i1,1)).and. & (ktab2(i2,2).ge.ktab1(i1,1))) then i2f=1 goto 322 c do i=k1,k2 c write(*,*) 'ktab1(',i,')=',(ktab1(i,j),j=1,3) c enddo endif enddo ! i2 322 continue if (i2f.eq.0) then write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'ktab1(',i1,',1)=',ktab1(i1,1), & ' could not be interpolated' write(*,*) 'nu1=',nu1 c Debug write(*,*) ktab2(nx2,1),ktab2(nx2,2) write(*,*) ((ktab2(nx2,1).le.ktab1(i1,1)).and. & (ktab2(nx2,2).ge.ktab1(i1,1))) write(*,*) 'kmax1=',kmax1,' kmax2=',kmax2, & ' kmax_min=',kmax_min write(*,*) 'nx1=',nx1,' nx2=',nx2 write(*,*) 'ktab1(',nx1,',1)=',ktab1(nx1,1), & ' ktab2(',nx2,',2)=',ktab2(nx2,2) do i=1,nx1 write(*,*) 'ktab1(',i,')=',(ktab1(i,j),j=1,3) enddo write(*,*) c write(*,*) 'nx2=',nx2 c do i2=1,nx2 c write(*,*) 'ktab2(',i2,',1)=',ktab2(i2,1),ktab2(i2,2) c call test_nan(ktab2(i2,1),out) c if (out.eq.1) then c stop c endif c enddo c stop x2df=0 goto 521 c Debug endif c First interpolation to determine values of ktab2p(i1,j) ktab2p(i1,1)=ktab1(i1,1) ktab2p(i1,2)=ktab1(i1,2) do j=3,3 c call lin_interp(ktab2(i2-1,1),ktab2(i2,1), c & ktab2p(i1,1),ktab2(i2-1,j),ktab2(i2,j), c & ktab2p(i1,j)) if (i2.le.nx2) then call lin_interp(ktab2(i2,1),ktab2(i2,2),ktab2p(i1,1), & ktab2(i2,j),ktab2(i2+1,j),ktab2p(i1,j)) else ktab2p(i1,j)=ktab2(i2,j) endif enddo ! j c Test if (ktab2p(i1,3).lt.ktab2p(i1,2)) then write(*,*) 'Error in first interpolation:' write(*,*) 'ktab1(',i1,')=',(ktab1(i1,j),j=1,3) write(*,*) 'ktab2(',i1,')=',(ktab2(i1,j),j=1,3) write(*,*) 'ktab2p(',i1,')=',(ktab2p(i1,j),j=1,3) write(*,*) 'nu1=',nu1 write(*,*) x2df=0 goto 521 endif c Test c write(*,*) 'ktab2p(',i1,')=',(ktab2p(i1,j),j=1,3) c Debug c Second interpolation to determine values of ktab_yi(i1,j) ktab_yi(i1,1)=ktab1(i1,1) ktab_yi(i1,2)=ktab1(i1,2) do j=3,3 call lin_interp(y1,y2,y,ktab1(i1,j),ktab2p(i1,j), & ktab_yi(i1,j)) if ((ktab_yi(i1,1).lt.0.0D+0). & and.(ktab_yi(i1,2).ge.0.0D+0)) then ktab_yi(i1,3)=0.0D+0 endif enddo ! j c Debug if (ktab_yi(i1,3).lt.ktab_yi(i1,2)) then write(*,*) 'Error in second interpolation:' write(*,*) 'ktab_yi(',i1,',2)=',ktab_yi(i1,2),' >', & ' ktab_yi(',i1,',3)=',ktab_yi(i1,3) write(*,*) 'y1=',y1,' y2=',y2,' y=',y write(*,*) 'ktab1(',i1,',3)=',ktab1(i1,3), & ' ktab2p(',i1,',3)=',ktab2p(i1,3) write(*,*) 'nu1=',nu1 write(*,*) x2df=0 goto 521 endif c Debug enddo ! i1 c Debug do i=2,nx1 if (ktab2p(i-1,2).ne.ktab2p(i,1)) then write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'ktab2p(',i-1,',2)=',ktab2p(i-1,2) write(*,*) 'ktab2p(',i,',1)=',ktab2p(i,1) write(*,*) 'nu1=',nu1 write(*,*) x2df=0 goto 521 endif if (ktab_yi(i-1,2).ne.ktab_yi(i,1)) then write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'ktab_yi(',i-1,',2)=',ktab_yi(i-1,2) write(*,*) 'ktab_yi(',i,',1)=',ktab_yi(i,1) write(*,*) 'nu1=',nu1 write(*,*) x2df=0 goto 521 endif enddo c Debug c Now, x1 has to be identified in ktab_yi x2f=0 do i=1,nx1 if ((x1.ge.ktab_yi(i,1)).and.(x1.lt.ktab_yi(i,2))) then x2f=1 x2d=ktab_yi(i,3) c Debug c write(*,*) 'i=',i c write(*,*) ktab1(i,3),ktab2p(i,3),ktab_yi(i,3) c Debug goto 323 endif enddo ! i 323 continue if (x2f.eq.0) then write(*,*) 'Error from routine find_nu2_tab' write(*,*) 'x2 could not be found in doppler mode' write(*,*) 'nu1=',nu1 write(*,*) 'x1=',x1 write(*,*) 'ktab_yi(1,1)=',ktab_yi(1,1), & ' ktab_yi(',nx1,',2)=',ktab_yi(nx1,2) do i=1,nx1 write(*,*) 'ktab_yi(',i,')=',(ktab_yi(i,j),j=1,3) enddo write(*,*) endif nu2d=x2d*gamma_D/dsqrt(dlog(2.0D+0))+nuc c Debug eps=dabs(dabs(nu2d/nu1)-1.0D+0) if ((eps.le.1.0D-10).or.(nu2d.le.nu1)) then x2dts=1 else x2dts=0 endif c Debug endif ! lorentz or doppler 521 continue c Also compute nu2 criterion based on Lorentz lineshape tabulation x1=(nu1-nuc)/gamma_L c Debug call test_inf(x1,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'x1=',x1 write(*,*) 'nu1=',nu1 write(*,*) 'nuc=',nuc write(*,*) 'gamma_L=',gamma_L stop endif c Debug if (dabs(x1).lt.1.0D-6) then x1=0.0D+0 endif x2f=0 i0f=0 do i=1,nacc-1 if ((x1.gt.xacc(i,1)).and.(x1.lt.xacc(i+1,1))) then i0=i i0f=1 goto 123 endif enddo 123 continue if (i0f.eq.0) then a=1.0D+0/(1.0D+0+x1**2)+2*x1**2/(1.0D+0+x1**2)**2 b=-2*x1/(1.0D+0+x1**2)**2 if (x1.gt.xtab(nxtab,2)) then sign=-1.0D+0 else if (x1.lt.xtab(1,1)) then sign=1.0D+0 endif x2=x1*(1+sign*eps_max)+sign*eps_max*a/b x2f=1 c Debug call test_nan(x2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'x2=',x2 write(*,*) 'x1=',x1 write(*,*) 'sign=',sign write(*,*) 'eps_max=',eps_max write(*,*) 'a=',a write(*,*) 'b=',b stop endif c Debug c Debug if (x2.lt.x1) then write(*,*) 'error in find_nu2_tab:' write(*,*) 'x1=',x1,' x2=',x2 endif c Debug else ! i0f=1 x2f=0 do i=int(xacc(i0,2)),int(xacc(i0+1,2)) if ((x1.ge.xtab(i,1)).and.(x1.lt.xtab(i,2))) then x2=xtab(i,3) x2f=1 c Debug call test_nan(x2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'x2=',x2 write(*,*) 'error2' stop endif c Debug goto 124 endif enddo 124 continue c Debug if (x2f.eq.0) then write(*,*) 'Error from routine find_nu2_tab' write(*,*) 'x2 could not be found in lorentz mode' write(*,*) 'nu1=',nu1 write(*,*) 'x1=',x1 write(*,*) 'i0=',i0 write(*,*) 'xacc(',i0,',1)=',xacc(i0,1), & ' xacc(',i0,',2)=',xacc(i0,2) write(*,*) 'xacc(',i0+1,',1)=',xacc(i0+1,1), & ' xacc(',i0+1,',2)=',xacc(i0+1,2) write(*,*) 'xtab(',xacc(i0,2),',1)=', & xtab(int(xacc(i0,2)),1) write(*,*) 'xtab(',xacc(i0+1,2),',2)=', & xtab(int(xacc(i0+1,2)),2) write(*,*) endif c Debug endif if (dabs(xtab(i,3)/x1-1.0D+0).lt.1.0D-5) then i=i+1 x2=xtab(i,3) c Debug call test_nan(x2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'x2=',x2 write(*,*) 'error3' stop endif c Debug endif nu2l=x2*gamma_L+nuc c Debug call test_nan(nu2l,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2l=',nu2l write(*,*) 'x2=',x2 write(*,*) 'gamma_L=',gamma_L write(*,*) 'nuc=',nuc stop endif c Debug c keep nu2 min between Doppler and Lorentz tabulation results if ((doppler.eq.1).and.(lorentz.eq.0).and.(nu2d.lt.nu2l) & .and.(x2dts.eq.0).and.(x2df.eq.1)) then c Debug c write something on screen in case it ever happens... c write(*,*) 'nu2 doppler has been retained !' c write(*,*) 'nu1=',nu1,' x1=',x1,' x2=',x2,' nu2=',nu2 c write(*,*) 'x2d=',x2d,' nu2d=',nu2d c Debug nu2=nu2d c Debug call test_nan(nu2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2=',nu2 write(*,*) 'nu2d=',nu2d stop endif c Debug else if (x2f.eq.1) then nu2=nu2l c Debug call test_nan(nu2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2=',nu2 write(*,*) 'nu2l=',nu2l stop endif c Debug else write(*,*) 'Error from routine find_nu2_tab:' write(*,*) 'x2 could not be computed at all !' write(*,*) nu2=nu1+dnu_fail c Debug call test_nan(nu2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2=',nu2 write(*,*) 'nu1=',nu1 write(*,*) 'dnu_fail=',dnu_fail stop endif c Debug endif c Debug eps=dabs(dabs(nu2/nu1)-1.0D+0) c In spite of all tests, it happens sometimes that eps<1.0D-10, which c is really really a huge issue because it means that nu2 is not greater c than nu1 (@ machine precision). c In this case, I can't find a solution but to add a very small amount to c nu1, so that the algorithm does not stay stuck at nu1. if (eps.lt.1.0D-10) then nu2=nu1+dnu_min ! cm¯¹, which is really small. c Debug call test_nan(nu2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2=',nu2 write(*,*) 'dnu_min=',dnu_min write(*,*) 'eps=',eps stop endif c Debug endif c Debug call test_nan(nu2,out) if (out.eq.1) then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'nu2=',nu2 stop endif c Debug return end subroutine read_dx(sd_mxep,nxtab,xtab) implicit none include 'max.inc' include 'formats.inc' double precision sd_mxep character*(Nchar_mx) file_dnu,nlines_file integer nxtab,ik,ios,strlen double precision xtab(1:Nxtab_mx,1:3) call lorentztab_filename(sd_mxep,file_dnu) c Debug c write(*,*) 'read_dx; file_dnu=',file_dnu(1:strlen(file_dnu)) c Debug c file_dnu='./data/dx_optimal.txt' nlines_file='./nlines' call get_nlines(file_dnu,nlines_file,nxtab) open(11,file=file_dnu(1:strlen(file_dnu)), & status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from routine read_dx' write(*,*) 'File could not be found:' write(*,*) file_dnu(1:strlen(file_dnu)) stop endif do ik=1,nxtab read(11,50) xtab(ik,1),xtab(ik,2),xtab(ik,3) enddo close(11) return end subroutine read_ktab(sd_mxep,ktab,nk,nx,yvalues,kindex) implicit none c c Purpose: to read the k(x,y) function tabulation file c c Output: c + ktab: values of [x1,x1p,x2] for each k function tabulated c + nk: number of k functions tabulated c + nx: number of [x1,x1p,x2] triplets for each k function tabulated c + yvalues: value of y for every k function tabulated c + kindex: c - kindex(i,1) is the index (in ktab) that locates the beginning of [x1,x1p,x2] sequence for yvalue 'i' c - kindex(i,2) is the index (in ktab) that locates the end of [x1,x1p,x2] sequence for yvalue 'i' c include 'max.inc' include 'formats.inc' integer strlen,i,j,ios character*(Nchar_mx) file_ktab,strtmp1,strtmp2,nlines_file integer nk,nx(1:Nytab_mx),nl,nlr,fin,nkr integer kindex(1:Nytab_mx,1:2) double precision sd_mxep double precision ktab(1:Nktab_mx,1:3) double precision yvalues(1:Nytab_mx) c file_ktab='./data/tabulation_voigt.txt' call voigttab_filename(sd_mxep,file_ktab) c Debug c write(*,*) 'file_ktab=',file_ktab(1:strlen(file_ktab)) c Debug nlines_file='./nlines' call get_nlines(file_ktab,nlines_file,nl) nlr=0 fin=0 nk=0 nkr=0 open(11,file=file_ktab(1:strlen(file_ktab)), & status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from routine read_ktab' write(*,*) 'File could not be found:' write(*,*) file_ktab(1:strlen(file_ktab)) stop endif do while (fin.eq.0) nk=nk+1 read(11,36) strtmp1,yvalues(nk),strtmp2,nx(nk) c Debug c write(*,*) 'yvalues(',nk,')=',yvalues(nk) c Debug kindex(nk,1)=nkr+1 do i=1,nx(nk) read(11,50) (ktab(nkr+i,j),j=1,3) c Debug c write(*,*) 'ktab(',nkr+i,')=',(ktab(nkr+i,j),j=1,3) c Debug enddo ! i read(11,*) nkr=nkr+nx(nk) kindex(nk,2)=nkr nlr=nlr+nx(nk)+2 if (nlr.ge.nl) then fin=1 endif enddo close(11) c Debug c write(*,*) 'nk=',nk c do i=1,nk c write(*,*) 'nx(',i,')=',nx(i),' kindex(',i,',1)=',kindex(i,1), c & ' kindex(',i,',2)=',kindex(i,2) c enddo c stop c Debug return end subroutine acceleration_dx(nxtab,xtab,nacc,xacc) implicit none include 'max.inc' integer nxtab,i,nint,nacc,ninc double precision xtab(1:Nxtab_mx,1:3) double precision xacc(1:Nacc_mx,1:2) nint=100 nacc=nint+1 ninc=nxtab/(nacc-1) xacc(1,1)=xtab(1,1) xacc(1,2)=1 do i=2,nacc-1 xacc(i,1)=xtab(ninc*(i-1),2) xacc(i,2)=ninc*(i-1) enddo xacc(nacc,1)=xtab(nxtab,2) xacc(nacc,2)=nxtab c Debug c do i=1,nacc c write(*,*) 'xacc(',i,',1)=',xacc(i,1), c & ' xacc(',i,',2)=',xacc(i,2) c enddo c stop c Debug return end subroutine find_nu2(nu1_in,nu0,delta,P,gamma,eps_max,nu2_out) implicit none c Purpose: to compute frequency nu2 given frequency nu1 and line width gamma c c Inputs: c + nu1_in: starting vavenumber [cm¯¹] c + nu0: vacuum wavenumber [cm¯¹] c + delta: air pressure-induced line shift [cm¯¹.atm¯¹] c + P: the total pressure [atm] c + gamma: line half-width [cm¯¹.atm¯¹] c + eps_max: maximum allowed relative error [0-1] between ka given by the Lorentz c profile and ka given by linear interpolation between nu1 and nu2 c c Output: c + nu2_out: ending wavenumber [cm¯¹]; this is the wavenumber at which c the maximum relative error "eps_max" has been reached c double precision nu1_in,nu1,nu0,delta,P,gamma,eps_max,nu2,nu2_out double precision nuc,x1,x2 double precision cb,cc,cd,ce,cp,cq,cr,cs,ct,cu double precision y0,a0,b0 integer nr,i double precision z0,z1,z2 integer ns,s,s0,nwn2,iwn2(1:4) double precision yr(1:3),xr(1:3) double precision z(1:4),x(1:4),wn2(1:4),dwn(1:4) double precision dwn_min integer sign,ntries double precision xm,val,eps,eps_xr double precision beta1,beta2,alpha double precision ai,bi,ci,di c Debug c write(*,*) 'coucou find_nu2' c Debug eps=1.0D-4 eps_xr=1.0D-6 nuc=nu0+delta*P nu1=nu1_in x1=(nu1-nuc)/gamma if (x1.lt.1.0D-5) then nu2=nu1+gamma*5.0D-2 goto 666 endif x2=x1+gamma/1.0D+2 xm=(x1+x2)/2.0D+0 val=(1.0D+0+xm**2)/(x2-x1)* & ((x2-xm)/(1.0D+0+x1**2)+ & (xm-x1)/(1.0D+0+x2**2)) if (val.gt.1) then sign=1 else sign=-1 endif c Debug c write(*,*) 'val=',val,' sign=',sign c Debug ntries=0 123 continue ntries=ntries+1 cb=2.0D+0*x1 cc=(2.0D+0-8.0D+0*(1.0D+0+sign*eps_max))*x1**2 & +6.0D+0-8.0D+0*(1.0D+0+sign*eps_max) cd=2.0D+0*x1*(2.0D+0+x1**2) ce=(x1**2+4.0D+0)*(2.0D+0+x1**2) & -8.0D+0*(1.0D+0+sign*eps_max)*(1.0D+0+x1**2) c Debug c write(*,*) 'cb=',cb c write(*,*) 'cc=',cc c write(*,*) 'cd=',cd c write(*,*) 'ce=',ce c Debug if (dabs(x1).lt.1.0D-2) then ns=1 x2=x1+1.0D-2 wn2(1)=x2*gamma+nuc else cp=-3.0D+0/8.0D+0*cb**2+cc cq=1.0D+0/8.0D+0*cb**3-1.0D+0/2.0D+0*cb*cc+cd cr=-3.0D+0*(cb/4.0D+0)**4+cc*(cb/4.0D+0)**2-cb*cd/4.0D+0+ce cs=-cp/2.0D+0 ct=-cr cu=cp*cr/2.0D+0-1.0D+0/8.0D+0*cq**2 call findroots_3rd(1.0D+0,cs,ct,cu,nr,yr) y0=yr(1) ! we only need one solution a0=dsqrt(-cp+2.0D+0*y0) b0=-cq/(2.0D+0*a0) ns=0 call findroots_2nd(1.0D+0,-a0,y0-b0,nr,z0,z1,z2) if (nr.eq.1) then ns=ns+1 z(ns)=z0 else if (nr.eq.2) then ns=ns+1 z(ns)=z1 ns=ns+1 z(ns)=z2 endif call findroots_2nd(1.0D+0,a0,y0+b0,nr,z0,z1,z2) if (nr.eq.1) then ns=ns+1 z(ns)=z0 else if (nr.eq.2) then ns=ns+1 z(ns)=z1 ns=ns+1 z(ns)=z2 endif do s=1,ns x(s)=z(s)-cb/4.0D+0 wn2(s)=gamma*x(s)+nuc enddo c Debug c write(*,*) 'ns=',ns c do s=1,ns c write(*,*) 'x(',s,')=',x(s),' wn2(',s,')=',wn2(s) c enddo c Debug endif ! x1=0 nwn2=0 do s=1,ns if ((wn2(s).gt.0.0D+0).and. & (x(s).gt.x1).and. & ((x(s)/x1).gt.0.0D+0)) then nwn2=nwn2+1 iwn2(nwn2)=s endif enddo if (ntries.eq.2) then goto 321 else if (nwn2.eq.0) then sign=-sign goto 123 endif endif 321 continue c Debug c write(*,*) 'ntries=',ntries,' nwn2=',nwn2 c do i=1,nwn2 c write(*,*) 'wn2(',iwn2(i),')=',wn2(iwn2(i)) c enddo c Debug c if you end up with exactly one solution if (nwn2.eq.1) then nu2=wn2(iwn2(nwn2)) endif c if ntries=2, there is an inversion between the given line profile c and the linear profile that uses the solution that has been found c => this solution is WRONG if ((ntries.eq.2).and.(nwn2.eq.1)) then x2=x(iwn2(1)) c Debug c write(*,*) 'nu2 found=',x2*gamma+nuc c Debug beta1=1.0D+0/(1.0D+0+x1**2) beta2=1.0D+0/(1.0D+0+x2**2) alpha=(beta2-beta1)/(x2-x1) ai=alpha bi=beta1-alpha*x1 ci=alpha di=beta1-alpha*x1-1.0D+0 c find xi, abscissa of the intersection between the real line profile c and the linear profile that uses the (wrong) solution call findroots_3rd(ai,bi,ci,di,nr,xr) c Debug c write(*,*) 'nr=',nr c do i=1,nr c write(*,*) 'xr(',i,')=',xr(i),' nur=',xr(i)*gamma+nuc c enddo c Debug s=0 do i=1,3 if ((dabs(xr(i)-x1).gt.eps_xr).and. & (dabs(xr(i)-x2).gt.eps_xr).and. & (xr(i).gt.x1).and.(xr(i).lt.x2)) then s=i goto 421 endif enddo 421 continue if (s.eq.0) then c write(*,*) 'no intersection found' nu2=(x1+x(iwn2(1)))/2.0D+0*gamma+nuc else nu2=xr(s)*gamma+nuc c write(*,*) 'intersection wn=',nu2 endif c Debug c write(*,*) 'nu2=',nu2 c stop c write(*,*) 'wn:',nu2 c write(*,*) 's=',s c write(*,*) 'xr(',s,')=',xr(s),' x1=',x1,' x2=',x2 c write(*,*) 'x(iwn2(1))=',x(iwn2(1)) c write(*,*) dabs(xr(s)-x1),eps_xr c Debug goto 666 endif c you might end up with zero solution with x1<0 if ((ntries.eq.2).and.(nwn2.eq.0) & .and.(dabs(x1).lt.0.3D+0)) then nu2=nu1+gamma*5.0D-2 if ((x1.lt.0.0D+0).and.(nu2.gt.nuc)) then nu2=nuc endif goto 666 endif c you might end up with two very close solutions if ((nwn2.eq.2).and. & (dabs(wn2(iwn2(1))/wn2(iwn2(2))-1.0D+0).lt.eps)) then nwn2=1 ! only one is retained endif c if even after that the number of solutions is not exactly 1 c if (nwn2.ne.1) then c write(*,*) 'Number of solutions:',nwn2 ! print the number of solutions c endif c and if there is still more than one solution c Debug c if (nwn2.gt.1) then c do s=1,nwn2 c write(*,*) 'nu2=',wn2(iwn2(s)) ! print every solution found c enddo c endif c Debug if (nwn2.gt.1) then do s=1,nwn2 dwn(iwn2(s))=wn2(iwn2(s))-nu1 c Debug c write(*,*) 'dwn(',s,')=',dwn(iwn2(s)) c Debug enddo do s=1,nwn2 if (dwn(iwn2(s)).gt.0.0D+0) then s0=s dwn_min=dwn(iwn2(s0)) goto 323 endif enddo 323 continue c Debug c write(*,*) 's0=',s0,' dwn_min=',dwn_min c Debug do s=1,nwn2 if ((dwn(iwn2(s)).lt.dwn_min).and. & (dwn(iwn2(s)).gt.0.0D+0)) then dwn_min=dwn(iwn2(s)) s0=s ! select solution that is the closest from nu1 c Debug c write(*,*) 's=',s,' s0=',s0,' dwn_min=',dwn_min c Debug endif enddo nu2=wn2(iwn2(s0)) nwn2=1 c Debug c write(*,*) 'selected solution:',wn2(iwn2(s0)) c Debug endif c in the case there is still not exactly 1 solution, if (nwn2.ne.1) then write(*,*) 'warning in find_nu2: nwn2=',nwn2 write(*,*) 'nu1=',nu1_in write(*,*) 'ntries=',ntries stop ! you have to stop endif if (nwn2.eq.1) then nu2=wn2(iwn2(nwn2)) endif 666 continue nu2_out=nu2 c Debug if ((dabs(nu2_out/nu1_in-1.0D+0).lt.1.0D-10) & .or.(nu2_out.lt.nu1_in)) then write(*,*) 'error in find_nu2:' write(*,*) 'nu1=',nu1_in write(*,*) 'nu2=',nu2_out write(*,*) 'nuc=',nuc,' nu1=',nu1,' x1=',x1 write(*,*) 'ntries=',ntries,' ns=',ns stop endif c write(*,*) 'ok find_nu2: nu1=',nu1_in,' nu2=',nu2_out c Debug return end subroutine findroots_3rd(a,b,c,d,nr,x) implicit none include './parameters.inc' c Purpose: to find ONLY ONE root of a specified c third-order polynomia double precision a,b,c,d,x(1:3),z(1:3) integer nr,i double precision p,q,delta double precision u,v,y(1:3) double precision epsilon c Debug c write(*,*) 'coucou findroots_3rd' c Debug epsilon=1.0D-2 c Debug c write(*,*) 'a=',a c write(*,*) 'b=',b c write(*,*) 'c=',c c write(*,*) 'd=',d c Debug p=-1.0D+0/3.0D+0*(b/a)**2+c/a q=1.0D+0/27.0D+0*(b/a)*(2.0D+0*(b/a)**2-9.0D+0*c/a)+d/a delta=q**2+4.0D+0/27.0D+0*p**3 c Debug c write(*,*) 'p=',p c write(*,*) 'q=',q c write(*,*) 'delta=',delta c Debug if (delta.gt.0.0D+0) then nr=1 call powalpha((-q+dsqrt(delta))/2.0D+0,1.0D+0/3.0D+0,u) call powalpha((-q-dsqrt(delta))/2.0D+0,1.0D+0/3.0D+0,v) z(1)=u+v c Debug c write(*,*) 'u=',u,u**3,(-q+dsqrt(delta))/2.0D+0 c write(*,*) 'v=',v,v**3,(-q-dsqrt(delta))/2.0D+0 c write(*,*) 'z0=',z0 c Debug else if (delta.eq.0.0D+0) then nr=2 z(1)=3.0D+0*q/p z(2)=-3.0D+0*q/(2.0D+0*p) else if (delta.lt.0.0D+0) then nr=3 do i=1,3 z(i)=2.0D+0*dsqrt(-p/3.0D+0)* & dcos(1.0D+0/3.0D+0* & dacos(-q/2.0D+0*dsqrt(27.0D+0/(-p**3))) & +2.0D+0*i*pi/3.0D+0) enddo endif do i=1,nr x(i)=z(i)-b/(3.0D+0*a) y(i)=a*x(i)**3+b*x(i)**2+c*x(i)+d c if (dabs(y(i)/x(i)).gt.epsilon) then c write(*,*) 'resolution polyn. 3e degré:' c write(*,*) 'x(',i,')=',x(i),' a*x³+b*x²+c*x+d=',y(i) c stop c endif enddo c Debug c write(*,*) 'ok findroots_3rd' c Debug return end subroutine findroots_2nd(a,b,c,nr,x0,x1,x2) implicit none c Purpose: to find ALL roots of a specified second-order polynomia integer nr double precision a,b,c,x0,x1,x2,delta c Debug c write(*,*) 'coucou findroots_2nd' c Debug delta=b**2-4.0D+0*a*c c Debug c write(*,*) 'in find_roots_2nd: a=',a c write(*,*) 'in find_roots_2nd: b=',b c write(*,*) 'in find_roots_2nd: c=',c c write(*,*) 'in find_roots_2nd: delta=',delta c Debug if (delta.eq.0) then nr=1 x0=-b/(2.0D+0*a) c Debug c write(*,*) 'in find_roots_2nd: nr=',nr c Debug else if (delta.gt.0.0D+0) then nr=2 x1=(-b-dsqrt(delta))/(2.0D+0*a) x2=(-b+dsqrt(delta))/(2.0D+0*a) c Debug c write(*,*) 'in find_roots_2nd: nr=',nr c Debug else if (delta.lt.0.0D+0) then nr=0 endif c Debug c write(*,*) 'ok findroots_2nd' c Debug return end subroutine powalpha(x,alpha,y) implicit none double precision x,alpha,y c Debug c write(*,*) 'coucou powalpha' c Debug if (x.ge.0.0D+0) then y=x**alpha else y=-(-x)**alpha endif c Debug c write(*,*) 'ok powalpha' c Debug return end subroutine find_fllines(Nlines,Nmol,mol_index,iso_index,mol_niso, & index,isotope,Ps,i,nu0,clines,nu_min,nu_max, & fline,lline) implicit none include 'max.inc' include 'parameters.inc' integer Nlines,line integer Nmol integer mol_index(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) integer mol_niso(1:Nmol_max) integer index(1:Nline_mx) integer isotope(1:Nline_mx) double precision Ps(1:Nmol_max) integer i,mol double precision nu0(1:Nline_mx) integer clines(1:Nline_mx) double precision nu_min,nu_max,dnu,nu1,nu2 integer fline,lline,line1,line2,nint integer l1f,l2f,isof integer molec,isotop nint=0 dnu=(nu_max-nu_min)*nint if (nu_min.gt.dnu) then nu1=nu_min-dnu else nu1=0.0D+0 endif nu2=nu_max+dnu l1f=0 c Debug c write(*,*) 'This is find_fllines' c write(*,*) 'Nlines=',Nlines c write(*,*) 'Nmol=',Nmol c do mol=1,Nmol c write(*,*) 'Ps(',mol,')=',Ps(mol) c enddo c do line=1,500 c write(*,*) 'nu0(',line,')=',nu0(line) c enddo c write(*,*) c stop c Debug do line=1,Nlines call indexes(Nmol,mol_index,iso_index,mol_niso, & index(line),isotope(line),molec,isotop,isof) c Debug c if (Ps(molec).eq.0.0D+0) then c write(*,*) 'Ps(',molec,')=',Ps(molec) c stop c endif c if ((nu0(line).gt.nu_min).and.(nu0(line).lt.nu_max)) then c write(*,*) 'nu0(',line,')=',nu0(line) c endif c Debug if ((nu0(line).gt.nu1) & .and.(l1f.eq.0) & .and.(Ps(molec).ne.0.0D+0)) then l1f=1 line1=line goto 121 endif enddo if (l1f.eq.0) then write(*,*) 'Error from routine find_fllines' write(*,*) 'No line found for nu>',nu1 write(*,*) 'Nlines=',Nlines write(*,*) 'Nmol=',Nmol do mol=1,Nmol write(*,*) 'Ps(',mol,')=',Ps(mol) enddo stop endif 121 continue l2f=0 do line=line1,Nlines call indexes(Nmol,mol_index,iso_index,mol_niso, & index(line),isotope(line),molec,isotop,isof) if ((nu0(line).gt.nu2) & .and.(l2f.eq.0) & .and.(Ps(molec).ne.0.0D+0)) then l2f=1 line2=line-1 goto 122 endif enddo if (l2f.eq.0) then write(*,*) 'Warning from routine find_fllines' write(*,*) 'No line found for nu<',nu2 c stop endif 122 continue if ((l1f.eq.0).and.(l2f.eq.0)) then line1=0 line2=0 endif fline=line1 lline=line2 c Debug c write(*,*) 'line1=',line1,' nu0(',line1,')=',nu0(line1) c write(*,*) 'fline=',fline,' nu0(',fline,')=',nu0(fline) c write(*,*) 'line2=',line2,' nu0(',line2,')=',nu0(line2) c write(*,*) 'lline=',lline,' nu0(',lline,')=',nu0(lline) c stop c Debug return end subroutine select_lines(nl,S,ppres,ind,nl2,ind2) implicit none include 'max.inc' integer nl,l,lmax,nl2 integer ind(1:Nlpnb_mx) integer ind2(1:Nlpnb_mx) double precision S(1:Nmax_tab) double precision ppres(1:Nlpnb_mx) double precision Smax call find_max(S,nl,Smax,lmax) c Debug c if (nl.eq.75) then c write(*,*) 'Smax=',Smax c endif c Debug nl2=0 do l=1,nl if ((S(l).ge.Smax/1.0D+2).and. & (ppres(l).ne.0.0D+0)) then nl2=nl2+1 ind2(nl2)=ind(l) c Debug c if (nl.eq.75) then c write(*,*) 'nl2=',nl2,' l=',l,' S(',l,')=',S(l) c endif c Debug endif enddo c Debug c write(*,*) 'Smax=',Smax,' lmax=',lmax c write(*,*) 'nl=',nl,' nl2=',nl2 c Debug return end subroutine adjust_mesh(nu_min,nu_max,nu,nk) implicit none include 'max.inc' integer i,n,nk,k1,k2,k1f,k2f double precision nu_min,nu_max double precision nu(1:Nkmx) k1f=0 do i=2,nk if ((nu_min.ge.nu(i-1)).and.(nu_min.lt.nu(i))) then k1=i k1f=1 goto 111 endif enddo 111 continue if (k1f.eq.0) then write(*,*) 'Error from routine adjust_mesh' write(*,*) 'nu_min=',nu_min,' could not be found' write(*,*) 'nk=',nk write(*,*) 'nu(1)=',nu(1) write(*,*) 'nu(',nk,')=',nu(nk) c Debug do i=1,nk write(*,*) 'nu(',i,')=',nu(i) enddo c Debug stop endif k2f=0 do i=k1,nk if ((nu_max.gt.nu(i-1)).and.(nu_max.le.nu(i))) then k2=i k2f=1 goto 112 endif enddo 112 continue if (k2f.eq.0) then write(*,*) 'Error from routine adjust_mesh' write(*,*) 'nu_max=',nu_max,' could not be found' write(*,*) 'nk=',nk write(*,*) 'nu(1)=',nu(1) write(*,*) 'nu(',nk,')=',nu(nk) stop endif nu(1)=nu_min n=1 do i=k1,k2-1 n=n+1 nu(n)=nu(i) enddo n=n+1 nu(n)=nu_max do i=1,n nu(i)=nu(i) enddo nk=n return end