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 band_limits(nu_min,nu_max,T,nlev,B_mxep,sdis, & Nb,nu_low,nu_hi) implicit none c c Purpose: to discretize the spectral region and provide narrowband limits c c Inputs: c + nu_min: minimum wavenumber of the spectral region to discretize [cm¯¹] c + nu_max: maximum wavenumber of the spectral region to discretize [cm¯¹] c + T: temperature for each atmospheric layer [K] c + nlev : number of atmospheric levels c + B_mxep: maximum error percentage over B for discretization (0-1) c c Outputs: c + Nb: number of narrowbands c + nu_low: lower limit of each narrowband [cm¯¹] c + nu_hi: upper limit of each narrowband [cm¯¹] c include 'max.inc' include 'formats.inc' integer sdis,i,nlev,strlen double precision nu_min,nu_max,B_mxep double precision T(1:Nmax) integer Nb,Nit_max,it double precision nu_low(1:Nbmx) double precision nu_hi(1:Nbmx) double precision nu2tab(1:Nmax_tab) double precision Bep(1:Nmax_tab),Bepm integer imin,ios,band,imax double precision nu,nu1,nu2,epsm,dnu,dnu0,nus,nbw character*(Nchar_mx) file_Berror,file_nu double precision B,lambda1,lambda2,B1,B2,Bmin integer Nsc,case double precision sc(1:100) call provide_nulimits(Nsc,sc) c in the case the spectral discretization has to be computed if (sdis.eq.1) then Nit_max=100 dnu0=1.0D+1 ! cm¯¹ Nb=0 write(*,*) 'Now discretizing [',nu_min,'-',nu_max, & ']cm¯¹ into...' c Debug goto 421 c Debug nu2=nu_min do while (nu2.lt.nu_max) Nb=Nb+1 if (Nb.gt.Nbmx) then write(*,*) 'Error from routine band_limits:' write(*,*) 'Nbmx=',Nbmx,' is not sufficient' write(*,*) 'increase this parameter in max.inc' stop endif nu_low(Nb)=nu2 do i=1,nlev nu1=nu_low(Nb) nu=nu1+dnu0 nus=nu1 do it=1,Nit_max call find_epsm(nu1,nu,T(i),epsm) c Debug c write(*,*) nu1,nu,epsm,B_mxep c Debug if (epsm.lt.0.99D+0*B_mxep) then dnu=nu-nus nus=nu nu=nu+dnu endif if (epsm.gt.B_mxep) then nu=(nu+nus)/2.0D+0 endif if ((epsm.gt.0.99D+0*B_mxep).and.(epsm.lt.B_mxep)) & then nu2tab(i)=nu goto 123 endif if ((epsm.lt.B_mxep).and.(nu.ge.nu_max)) then nu2tab(i)=nu goto 123 endif enddo if (it.ge.Nit_max) then write(*,*) 'Error from routine band_limits:' write(*,*) 'Nit_max reached' stop endif 123 continue c Debug c write(*,*) i,nu2tab(i) c Debug enddo ! i call find_min(nu2tab,nlev,nu2,imin) if (nu2.gt.nu_max) then nu_hi(Nb)=nu_max else nu_hi(Nb)=nu2 endif dnu0=1.1D+0*(nu2-nu1) c Debug c write(*,*) 'nu_low(',Nb,')=',nu_low(Nb), c & ' nu_hi(',Nb,')=',nu_hi(Nb) c Debug enddo c Debug 421 continue Nb=1 nu_low(1)=nu_min call choose_nbw(nu_low(1),nbw) nu_hi(1)=nu_min+nbw do while (nu_hi(Nb).lt.nu_max) Nb=Nb+1 c Tests if (Nb.gt.Nbmx) then write(*,*) 'Error from band_limits' write(*,*) 'Nbmx has been exceeded' stop endif c Tests nu_low(Nb)=nu_hi(Nb-1) call choose_nbw(nu_low(Nb),nbw) nu_hi(Nb)=nu_low(Nb)+nbw c Debug c write(*,*) nu_low(Nb),nbw,nu_hi(Nb) c Debug c special cases do case=1,Nsc if ((nu_low(Nb).lt.sc(case)).and. & (nu_hi(Nb).gt.sc(case))) then nu_hi(Nb)=sc(case) endif enddo ! case enddo nu_hi(Nb)=nu_max c Nb=1 c nu_low(1)=nu_min c nu_hi(1)=nu_max open(10,file='./bands.txt') write(10,*) Nb do band=1,Nb write(10,*) nu_low(band),nu_hi(band) enddo close(10) c Debug if (Nb.eq.1) then write(*,*) '...',Nb,' narrowband interval' else write(*,*) '...',Nb,' narrowband intervals' endif c in the case a predefined spectral discretization has to be used else ! sdis=2 file_nu='./data/narrowbands.in' open(10,file=file_nu(1:strlen(file_nu)), & status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from routine band_limits:' write(*,*) 'File could not ne found:' write(*,*) file_nu(1:strlen(file_nu)) stop endif read(10,*) Nb if (Nb.gt.Nbmx) then write(*,*) 'Error from routine band_limits:' write(*,*) 'Nb=',Nb write(*,*) 'while Nbmx=',Nbmx stop endif do band=1,Nb read(10,*) nu_low(band),nu_hi(band) enddo close(10) c small verification if (Nb.gt.1) then do band=2,Nb if (nu_low(band).ne.nu_hi(band-1)) then write(*,*) 'Error from routine band_limits' write(*,*) 'Data error in file:' write(*,*) file_nu(1:strlen(file_nu)) write(*,*) 'nu_low(',band,')=',nu_low(band), & ' while nu_hi(',band-1,')=',nu_hi(band-1) stop endif enddo endif ! Nb > 1 write(*,*) 'Narrowband discretization read from file:' write(*,*) file_nu(1:strlen(file_nu)) nu_min=nu_low(1) nu_max=nu_hi(Nb) file_Berror='./B_error.txt' open(11,file=file_Berror(1:strlen(file_Berror))) do band=1,Nb ! narrowbands do i=1,nlev ! atmospheric levels call inversecm2micrometer(nu_low(band),lambda1) call inversecm2micrometer(nu_hi(band),lambda2) B1=B(T(i),lambda1) B2=B(T(i),lambda2) Bmin=B1 if (B2.lt.B1) then Bmin=B2 endif Bep(i)=dabs(B2-B1)/Bmin enddo ! i call find_max(Bep,nlev,Bepm,imax) write(11,*) band,nu_low(band),nu_hi(band),Bepm, & ' T(',imax,')=',T(imax) enddo ! band close(11) write(*,*) 'Errors associated to the narrowband', & ' discretization have been recorded in:' write(*,*) file_Berror(1:strlen(file_Berror)) endif ! sdis return end subroutine find_epsm(nu1,nu,T,epsm) implicit none include 'max.inc' double precision B,T double precision nu,nu1,lambda1,B1 double precision epsm,nul,lambdal,Bl,Bmin double precision eps(1:Nmax_tab) integer N,i,imax c Debug integer out c Debug c Debug c write(*,*) 'nu1=',nu1,' nu=',nu,' T=',T if (nu.lt.nu1) then write(*,*) 'This is routine find_epsm' write(*,*) 'nu1=',nu1,' nu=',nu endif c Debug N=nint((nu-nu1)*1.0D+3) if (N.gt.Nmax_tab) then do while (N.gt.Nmax) N=N/10 enddo endif call inversecm2micrometer(nu1,lambda1) B1=B(T,lambda1) c Debug c write(*,*) 'nu1=',nu1,' lambda1=',lambda1,' B1=',B1 c stop c Debug do i=1,N nul=nu1+(nu-nu1)/(N-1)*(i-1) call inversecm2micrometer(nul,lambdal) Bl=B(T,lambdal) Bmin=B1 if (Bl.lt.Bmin) then Bmin=Bl endif eps(i)=dabs(Bl-B1)/Bmin c Debug call test_inf(eps(i),out) if (out.eq.1) then write(*,*) 'B1=',B1,' Bl=',Bl,' Bmin=',Bmin, & ' eps(',i,')=',eps(i) write(*,*) 'nu1',nu1,' B1=',B1,' nul=',nul,' Bl=',Bl, & ' nu=',nu,' B=',B(T,1.0D+4/nu) stop endif c Debug enddo call find_max(eps,N,epsm,imax) return end