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 sort_all_data1(Nlines) implicit none include 'max.inc' include 'arrays1.inc' include 'arrays3.inc' c c Purpose: to sort the 9 input arrays in the order of increasing nu0 c integer Nlines integer l integer strlen character*(Nchar_mx) label label='subroutine sort_all_data1' c sort index(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=dble(index(l)) enddo call ssort(Nlines) do l=1,Nlines index(l)=int(a(l,2)) enddo c sort isotope(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=dble(isotope(l)) enddo call ssort(Nlines) do l=1,Nlines isotope(l)=int(a(l,2)) enddo c sort Sref(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=Sref(l) enddo call ssort(Nlines) do l=1,Nlines Sref(l)=a(l,2) enddo c sort gamma_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=gamma_air(l) enddo call ssort(Nlines) do l=1,Nlines gamma_air(l)=a(l,2) enddo c sort gamma_self(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=gamma_self(l) enddo call ssort(Nlines) do l=1,Nlines gamma_self(l)=a(l,2) enddo c sort Elow(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=Elow(l) enddo call ssort(Nlines) do l=1,Nlines Elow(l)=a(l,2) enddo c sort n_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=n_air(l) enddo call ssort(Nlines) do l=1,Nlines n_air(l)=a(l,2) enddo c sort delta_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=delta_air(l) enddo call ssort(Nlines) do l=1,Nlines delta_air(l)=a(l,2) enddo c sort n_self(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=n_self(l) enddo call ssort(Nlines) do l=1,Nlines n_self(l)=a(l,2) enddo c sort nu0(l) do l=1,Nlines nu0(l)=a(l,1) enddo return end subroutine sort_all_data2(Nlines,nu0,index,isotope,Sref, & gamma_air,gamma_self,Elow,n_air,delta_air,n_self) implicit none include 'max.inc' include 'arrays3.inc' c c Purpose: to sort the 9 input arrays in the order of increasing nu0 c integer Nlines 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 delta_air(1:Nlpnb_mx) double precision n_self(1:Nlpnb_mx) integer l integer strlen character*(Nchar_mx) label label='subroutine sort_all_data2' c sort index(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=dble(index(l)) enddo call ssort(Nlines) do l=1,Nlines index(l)=int(a(l,2)) enddo c sort isotope(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=dble(isotope(l)) enddo call ssort(Nlines) do l=1,Nlines isotope(l)=int(a(l,2)) enddo c sort Sref(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=Sref(l) enddo call ssort(Nlines) do l=1,Nlines Sref(l)=a(l,2) enddo c sort gamma_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=gamma_air(l) enddo call ssort(Nlines) do l=1,Nlines gamma_air(l)=a(l,2) enddo c sort gamma_self(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=gamma_self(l) enddo call ssort(Nlines) do l=1,Nlines gamma_self(l)=a(l,2) enddo c sort Elow(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=Elow(l) enddo call ssort(Nlines) do l=1,Nlines Elow(l)=a(l,2) enddo c sort n_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=n_air(l) enddo call ssort(Nlines) do l=1,Nlines n_air(l)=a(l,2) enddo c sort delta_air(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=delta_air(l) enddo call ssort(Nlines) do l=1,Nlines delta_air(l)=a(l,2) enddo c sort n_self(l) do l=1,Nlines a(l,1)=nu0(l) a(l,2)=n_self(l) enddo call ssort(Nlines) do l=1,Nlines n_self(l)=a(l,2) enddo c sort nu0(l) do l=1,Nlines nu0(l)=a(l,1) enddo return end subroutine sort_precalc_data(Nlines,index, & sorted_index,nuc,dnuc_dP, & gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx, & gamma_D,dgammaD_dT, & S,dS_dT,densities, & drho_dP,drho_dT,drho_dx) implicit none include 'max.inc' include 'arrays3.inc' c c Purpose: to sort the 14 output arrays of "precalc" in the order of decreasing line intensity c c I/O integer Nlines integer index(1:Nline_mx) integer sorted_index(1:Nline_mx) double precision nuc(1:Nline_mx) double precision dnuc_dP(1:Nline_mx) double precision gamma_L(1:Nline_mx) double precision dgammaL_dP(1:Nline_mx) double precision dgammaL_dT(1:Nline_mx) double precision dgammaL_dx(1:Nline_mx) double precision gamma_D(1:Nline_mx) double precision dgammaD_dT(1:Nline_mx) double precision S(1:Nline_mx) double precision dS_dT(1:Nline_mx) double precision densities(1:Nline_mx) double precision drho_dP(1:Nline_mx) double precision drho_dT(1:Nline_mx) double precision drho_dx(1:Nline_mx) c temp integer l c label integer strlen character*(Nchar_mx) label label='subroutine sort_precalc_data' if (Nlines.gt.Nline_mx) then write(*,*) 'Error from: ',label(1:strlen(label)) write(*,*) 'Nlines=',Nlines write(*,*) 'while Nline_mx=',Nline_mx stop endif c sort index do l=1,Nlines a(l,1)=S(l) a(l,2)=dble(index(l)) enddo call ssort(Nlines) do l=Nlines,1,-1 sorted_index(Nlines-l+1)=int(a(l,2)) enddo c sort nuc do l=1,Nlines a(l,1)=S(l) a(l,2)=nuc(l) enddo call ssort(Nlines) do l=Nlines,1,-1 nuc(Nlines-l+1)=a(l,2) enddo c sort dnuc_dP do l=1,Nlines a(l,1)=S(l) a(l,2)=dnuc_dP(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dnuc_dP(Nlines-l+1)=a(l,2) enddo c sort gamma_L do l=1,Nlines a(l,1)=S(l) a(l,2)=gamma_L(l) enddo call ssort(Nlines) do l=Nlines,1,-1 gamma_L(Nlines-l+1)=a(l,2) enddo c sort dgammaL_dP do l=1,Nlines a(l,1)=S(l) a(l,2)=dgammaL_dP(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dgammaL_dP(Nlines-l+1)=a(l,2) enddo c sort dgammaL_dT do l=1,Nlines a(l,1)=S(l) a(l,2)=dgammaL_dT(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dgammaL_dT(Nlines-l+1)=a(l,2) enddo c sort dgammaL_dx do l=1,Nlines a(l,1)=S(l) a(l,2)=dgammaL_dx(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dgammaL_dx(Nlines-l+1)=a(l,2) enddo c sort gamma_D do l=1,Nlines a(l,1)=S(l) a(l,2)=gamma_D(l) enddo call ssort(Nlines) do l=Nlines,1,-1 gamma_D(Nlines-l+1)=a(l,2) enddo c sort dgammaD_dT do l=1,Nlines a(l,1)=S(l) a(l,2)=dgammaD_dT(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dgammaD_dT(Nlines-l+1)=a(l,2) enddo c sort dS_dT do l=1,Nlines a(l,1)=S(l) a(l,2)=dS_dT(l) enddo call ssort(Nlines) do l=Nlines,1,-1 dS_dT(Nlines-l+1)=a(l,2) enddo c sort densities do l=1,Nlines a(l,1)=S(l) a(l,2)=densities(l) enddo call ssort(Nlines) do l=Nlines,1,-1 densities(Nlines-l+1)=a(l,2) enddo c sort drho_dP do l=1,Nlines a(l,1)=S(l) a(l,2)=drho_dP(l) enddo call ssort(Nlines) do l=Nlines,1,-1 drho_dP(Nlines-l+1)=a(l,2) enddo c sort drho_dT do l=1,Nlines a(l,1)=S(l) a(l,2)=drho_dT(l) enddo call ssort(Nlines) do l=Nlines,1,-1 drho_dT(Nlines-l+1)=a(l,2) enddo c sort drho_dx do l=1,Nlines a(l,1)=S(l) a(l,2)=drho_dx(l) enddo call ssort(Nlines) do l=Nlines,1,-1 drho_dx(Nlines-l+1)=a(l,2) enddo c sort S(l) do l=Nlines,1,-1 S(Nlines-l+1)=a(l,1) enddo return end subroutine intensity_distribution(Nlines,S, & Nbin,bin,cdf_indexes, & cdf,first_bin,last_bin) implicit none include 'max.inc' c c Purpose: to compute the cumulated distribution of line intensities c c Input: c + Nlines: number of lines c + S: intensities of the "Nlines" lines c c Output: c + Nbin: number of S bins c + bin: values of the "Nbin" S bins c + cdf_indexes: index of the first and last line for each bin c + cdf: cumulated distribution of line intensities c + first_bin: index of the first bin that must be used c + lasst_bin: index of the last bin that must be used c c I/O integer Nlines double precision S(1:Nline_mx) integer Nbin double precision bin(1:Nbin_mx,1:2) integer cdf_indexes(1:Nbin_mx,1:2) double precision cdf(1:Nbin_mx) integer first_bin integer last_bin c temp integer npd integer i,j,line integer Nl(1:Nbin_mx) integer out_code logical binf c label integer strlen character*(Nchar_mx) label label='subroutine intensity_distribution' c Number of intervals per decade npd=4 out_code=1 do while (out_code.eq.1) call initialize_bin(Nlines,S,npd,Nbin,bin) do i=1,Nbin Nl(i)=0 enddo ! i do line=1,Nlines call binify(Nbin,bin,S(line),npd,Nl) enddo ! line call analyze_number_distribution(Nbin,bin,Nl,npd,out_code) enddo ! while (out_code=1) c remove empty bins call remove_empty_bins(Nbin,bin,Nl) c invert order call reverse_bins(Nbin,bin,Nl) c Debug c do i=1,Nbin c write(*,*) 'Nl(',i,')=',Nl(i), c & '[',bin(i,1),';',bin(i,2),']' c enddo ! i c stop c Debug c produce cdf call make_cdf(Nbin,bin,Nl,cdf_indexes,cdf) c identify first and last bins binf=.false. do i=1,Nbin if (cdf(i).gt.0.0D+0) then binf=.true. first_bin=i goto 111 endif enddo ! i 111 continue if (.not.binf) then call error(label) write(*,*) 'first_bin was not found:' do i=1,Nbin write(*,*) 'bin(',i,')=',(bin(i,j),j=1,2) enddo ! i stop endif binf=.false. if (Nbin.eq.1) then if (cdf(Nbin).eq.1.0D+0) then binf=.true. last_bin=Nbin goto 112 endif else ! Nbin.ne.1 do i=Nbin,2,-1 if ((cdf(i).eq.1.0D+0).and.(cdf(i-1).lt.1.0D+0)) then binf=.true. last_bin=i goto 112 endif enddo ! i endif 112 continue if (.not.binf) then call error(label) write(*,*) 'last_bin was not found:' do i=1,Nbin write(*,*) 'bin(',i,')=',(bin(i,j),j=1,2) enddo ! i do i=1,Nbin write(*,*) 'cdf(',i,')=',cdf(i) enddo ! i stop endif c Debug c open(11,file='./bins.txt') c write(11,*) 'Nbin=',Nbin c write(11,*) 'first_bin=',first_bin c write(11,*) 'last_bin=',last_bin c write(11,*) c write(11,*) 'bin:' c do i=1,Nbin c write(11,*) 'i=',i,' bin=',(bin(i,j),j=1,2) c enddo ! i c write(11,*) 'cdf_indexes:' c do i=1,Nbin c write(11,*) 'i=',i,' cdf_indexes=',(cdf_indexes(i,j),j=1,2) c enddo ! i c write(11,*) 'cdf:' c do i=1,Nbin c write(11,*) 'i=',i,' cdf=',cdf(i) c enddo ! i c close(11) c Debug return end subroutine initialize_bin(Nlines,S,npd,Nbin,bin) implicit none include 'max.inc' c c Purpose: to initialize the "bin" array c c Input: c + Nlines: number of lines c + S: line intensities c + npd: numher of bins per decade c c Output: c + Nbin: number of bins c + bin: limits of the "Nbin" bins c c I/O integer Nlines double precision S(1:Nline_mx) integer npd integer Nbin double precision bin(1:Nbin_mx,1:2) c temp integer line,i,j double precision Smin,Smax integer n0 double precision x(1:Npd_mx,1:2) c label integer strlen character*(Nchar_mx) label label='subroutine initialize_bin' Smax=0.0D+0 do line=1,Nlines if (S(line).gt.Smax) then Smax=S(line) endif enddo ! line Smin=Smax do line=1,Nlines if (S(line).lt.Smin) then Smin=S(line) endif enddo ! line c Debug c write(*,*) 'Smin=',Smin c write(*,*) 'Smax=',Smax c Debug if (Smax.eq.0.0D+0) then call error(label) write(*,*) 'Smax=',Smax write(*,*) 'Nlines=',Nlines c Debug c open(12,file='./optimizations/intensities.txt') c do line=1,Nlines c write(12,*) line,S(line) c enddo ! line c close(12) c Debug stop else call identify_decade(Smax,n0) endif c Debug c write(*,*) 'n0=',n0 c Debug call decade(n0,npd,x) Nbin=npd do i=1,Nbin do j=1,2 bin(i,j)=x(i,j) enddo ! j enddo ! i return end subroutine identify_decade(x,n0) implicit none include 'max.inc' c c Purpose: to identify the decade 'x' belongs to c c Input: c + x: value of x c c Output: c + n0: first exponent of the decade 'x' belongs to; 'x' belongs to the decade [10^n0;10^(n0+1)] c c I/O double precision x integer n0 c temp integer n1 double precision x1,x2 c label integer strlen character*(Nchar_mx) label label='subroutine identify_decade' if (x.eq.0.0D+0) then call error(label) write(*,*) 'Bad input argument:' write(*,*) 'x=',x stop else if (x.ge.1.0D+0) then n1=0 else n1=1 endif n0=int(dlog(x)/dlog(10.0D+0))-n1 x1=(10.0D+0)**n0 x2=(10.0D+0)**(n0+1) if ((x.lt.x1).and.(dabs(x1-x).lt.1.0D-10)) then x=x1 else if ((x.gt.x2).and.(dabs(x-x2).lt.1.0D-10)) then x=x2 endif if ((x.lt.x1).or.(x.gt.x2)) then call error(label) write(*,*) 'x=',x write(*,*) 'n0=',n0 write(*,*) 'decade=[',x1,';',x2,']' write(*,*) 'x < x1=',x.lt.x1 write(*,*) 'x > x2=',x.gt.x2 stop c Debug c else c write(*,*) 'x=',x c write(*,*) 'decade=[',x1,';',x2,'] is ok' c Debug endif endif return end subroutine decade(n0,npd,x) implicit none include 'max.inc' c c Purpose: to discretize a decade c c Input: c + n0: x0=10^n0 is the lower x-limit of the decade; the decade limits are therefore [10^n0;10^(n0+1)] c + npd: number of intervals per decade c c Output: c + x: limits of the 'npd' intervals c c I/O integer n0 integer npd double precision x(1:Npd_mx,1:2) c temp integer int double precision alpha double precision x0 integer j c label integer strlen character*(Nchar_mx) label label='subroutine decade' if (npd.le.0) then call error(label) write(*,*) 'Bad input argument:' write(*,*) 'npd=',npd stop endif if (npd.gt.Npd_mx) then call error(label) write(*,*) 'npd=',npd write(*,*) 'while Npd_mx=',Npd_mx stop endif x0=(10.0D+0)**n0 alpha=dexp(dlog(10.0D+0)/dble(npd)) x(1,1)=x0 if (npd.gt.1) then do int=1,npd-1 x(int,2)=x0*alpha**int x(int+1,1)=x(int,2) enddo ! int endif x(npd,2)=x0*alpha**npd c Debug c do int=1,npd c write(*,*) (x(int,j),j=1,2) c enddo ! int c Debug return end subroutine binify(Nbin,bin,S,npd,Nl) implicit none include 'max.inc' c c Purpose: to classify a given S into one of the S bins c c Input: c + Nbin: number of S bins c + bin: values of S for each S bin c + S: current value of S c + npd: number of intervals per decade c + Nl: number of lines in each bin c c Output: c + Nl: updated array c c I/O integer Nbin double precision bin(1:Nbin_mx,1:2) double precision S integer npd integer Nl(1:Nbin_mx) c temp integer i,j,ibin logical bin_found integer n0 c label integer strlen character*(Nchar_mx) label label='subroutine binify' if (((S.lt.bin(1,1)).or.(S.gt.bin(Nbin,2))) & .and.(S.gt.0.0D+0)) then c Debug c write(*,*) 'calling identify_decade with S=',S c Debug call identify_decade(S,n0) c Debug c write(*,*) 'n0=',n0 c Debug call extend_bins(Nbin,bin,npd,n0,Nl) c Debug c write(*,*) 'bin was extended to decade:',n0 c Debug endif bin_found=.false. do i=1,Nbin if ((S.ge.bin(i,1)) & .and.(S.le.bin(i,2))) then bin_found=.true. ibin=i goto 111 endif enddo ! i 111 continue if (.not.bin_found) then call error(label) write(*,*) 'bin was not found:' write(*,*) 'S=',S do i=1,Nbin write(*,*) 'i:',i, & ' bin=',(bin(i,j),j=1,2) enddo ! i stop else c Add one to the number of lines in bin Nl(ibin)=Nl(ibin)+1 endif return end subroutine extend_bins(Nbin,bin,npd,n0,Nl) implicit none include 'max.inc' c c Purpose: to extend the S mesh up to a given decade c c Input: c + Nbin: number of bins in the initial S mesh c + bin: initial mesh c + npd: number of interals per decade c + n0: index of the decade [10^n0;10^(n0+1)] representing the new limit of the mesh c + Nl: number of lines in each S bin c c Output: c + Nbin, bin, Nl: updated c c I/O integer Nbin double precision bin(1:Nbin_mx,1:2) integer npd integer n0 integer Nl(1:Nbin_mx) c temp double precision x(1:Npd_mx,1:2) double precision x0 logical keep_going integer direction,index double precision x1,x2 integer first_index,i,j integer n1,n2 c label integer strlen character*(Nchar_mx) label label='subroutine extend_bins' x1=(bin(1,1)+bin(1,2))/2.0D+0 if (x1.eq.0.0D+0) then call error(label) write(*,*) 'x1=',x1 write(*,*) 'bin(1,1)=',bin(1,1) write(*,*) 'bin(1,2)=',bin(1,2) stop endif call identify_decade(x1,n1) x2=(bin(Nbin,1)+bin(Nbin,2))/2.0D+0 if (x2.eq.0.0D+0) then call error(label) write(*,*) 'x2=',x2 write(*,*) 'bin(',Nbin,',1)=',bin(Nbin,1) write(*,*) 'bin(',Nbin,',2)=',bin(Nbin,2) stop endif call identify_decade(x2,n2) if (n0.lt.n1) then direction=-1 first_index=n1 else if (n0.gt.n2) then direction=1 first_index=n2 else c call error(label) c write(*,*) 'n0=',n0 c write(*,*) 'x0=',(10.0D+0)**n0 c write(*,*) 'bin=[',bin(1,1),';',bin(Nbin,2),']' c write(*,*) 'No need to add any bin' c stop goto 666 endif c 'direction' is the direction decades must be added c 'first_index' is the exponent of the first or last decade in 'bin' index=first_index keep_going=.true. do while (keep_going) index=index+direction call decade(index,npd,x) if (Nbin+npd.gt.Nbin_mx) then call error(label) write(*,*) 'Nbin=',Nbin write(*,*) 'can not add ',npd,'bins' write(*,*) 'with Nbin_mx=',Nbin_mx stop endif if (direction.eq.-1) then do i=Nbin+npd,npd+1,-1 Nl(i)=Nl(i-npd) do j=1,2 bin(i,j)=bin(i-npd,j) enddo ! j enddo ! i do i=1,npd Nl(i)=0 do j=1,2 bin(i,j)=x(i,j) enddo ! j enddo ! i else ! direction=1 do i=1,npd Nl(Nbin+i)=0 do j=1,2 bin(Nbin+i,j)=x(i,j) enddo ! j enddo ! i endif ! direction Nbin=Nbin+npd if (index.eq.n0) then keep_going=.false. endif enddo ! while (keep_going) 666 continue return end subroutine analyze_number_distribution(Nbin,bin,Nl,npd,out_code) implicit none include 'max.inc' c c Purpose: to analyze the number distribution c in order to take appropriate actions when required c c Input: c + Nbin: number of bins c + bin: S limits of the 'Nbin' bins c + Nl: number of lines in each bin c + npd: number of intervals per decade c c Output: c + npd: updated, when needed c + out_code: output code (integer value) c c I/O integer Nbin double precision bin(1:Nbin_mx,1:2) integer Nl(1:Nbin_mx) integer npd integer out_code c temp integer i c label integer strlen character*(Nchar_mx) label label='subroutine analyze_number_distribution' c analyze number distribution out_code=0 ! 0: default value (no problem) do i=1,Nbin if (Nl(i).gt.Nline_mx) then call warning(label) write(*,*) 'Nl(',i,')=',Nl(i) out_code=1 ! 1: too many lines for some bins endif enddo ! bin if (out_code.eq.1) then npd=npd+1 write(*,*) 'Now trying to reduce number of lines by using:' write(*,*) 'Number of intervals per decade=',npd if (npd.gt.Npd_mx) then call error(label) write(*,*) 'npd=',npd write(*,*) 'while Npd_mx=',Npd_mx stop endif endif return end subroutine reverse_bins(Nbin,bin,Nl) implicit none include 'max.inc' c c Purpose: to reverse bins order c c I/O c + Nbin: number of bins c + bin: S limits of the 'Nbin' bins c + Nl: number of lines in each bin c integer Nbin double precision bin(1:Nbin_mx,1:2) integer Nl(1:Nbin_mx) c temp integer i,j double precision tbin(1:Nbin_mx,1:2) integer tNl(1:Nbin_mx) c label integer strlen character*(Nchar_mx) label label='subroutine reverse_bins' c copy orignial arrays do i=1,Nbin tNl(i)=Nl(i) do j=1,2 tbin(i,j)=bin(i,j) enddo ! j enddo ! i c reverse original arrays do i=1,Nbin Nl(i)=tNl(Nbin-i+1) do j=1,2 bin(i,j)=tbin(Nbin-i+1,j) enddo ! j enddo ! i return end subroutine remove_empty_bins(Nbin,bin,Nl) implicit none include 'max.inc' c c Purpose: to remove empty bins c c I/O c + Nbin: number of bins c + bin: S limits of the 'Nbin' bins c + Nl: number of lines in each bin c integer Nbin double precision bin(1:Nbin_mx,1:2) integer Nl(1:Nbin_mx) c temp integer i,j,ibin logical ebf,keep_going c label integer strlen character*(Nchar_mx) label label='subroutine remove_empty_bins' c Debug c write(*,*) label(1:strlen(label)) c write(*,*) 'in Nbin=',Nbin c Debug keep_going=.true. do while (keep_going) ebf=.false. do i=1,Nbin if (Nl(i).eq.0) then ebf=.true. ibin=i goto 111 endif enddo ! i 111 continue c Debug c write(*,*) 'ebf=',ebf c if (ebf) then c write(*,*) 'ibin=',ibin c endif c Debug if (ebf) then if (ibin.eq.1) then do i=1,Nbin-1 Nl(i)=Nl(i+1) do j=1,2 bin(i,j)=bin(i+1,j) enddo ! j enddo ! i else if (ibin.lt.Nbin) then do i=ibin,Nbin-1 Nl(i)=Nl(i+1) if (i.gt.ibin) then bin(i,1)=bin(i+1,1) endif bin(i,2)=bin(i+1,2) enddo ! i endif Nbin=Nbin-1 else keep_going=.false. endif enddo ! while (keep_going) c Debug c write(*,*) 'out Nbin=',Nbin c Debug return end subroutine make_cdf(Nbin,bin,Nl,cdf_indexes,cdf) implicit none include 'max.inc' c c Purpose: to produce the final cdf histogram c c Input: c + Nbin: number of bins c + bin: S limits of the 'Nbin' bins c + Nl: number of lines in each bin c c Output: c + cdf_indexes: indexes of the first and last line, for each bin c + cdf: cdf histogram c c I/O integer Nbin double precision bin(1:Nbin_mx,1:2) integer Nl(1:Nbin_mx) integer cdf_indexes(1:Nbin_mx,1:2) double precision cdf(1:Nbin_mx) c temp integer i,j,Nc,Nlines,Nc_prev character*(Nchar_mx) Nl_file character*(Nchar_mx) hfile character*(Nchar_mx) hist_file c label integer strlen character*(Nchar_mx) label label='subroutine make_histogram' Nlines=0 do i=1,Nbin Nlines=Nlines+Nl(i) enddo ! i Nc=0 cdf_indexes(1,1)=0 do i=1,Nbin Nc_prev=Nc Nc=Nc+Nl(i) if (Nc.gt.Nc_prev) then cdf_indexes(i,1)=Nc_prev+1 cdf_indexes(i,2)=Nc else cdf_indexes(i,1)=Nc_prev cdf_indexes(i,2)=Nc endif cdf(i)=dble(Nc)/dble(Nlines) enddo ! i c record number distribution for current species Nl_file='./optimizations/' & //'Number_distribution.txt' open(15,file=Nl_file(1:strlen(Nl_file))) do i=1,Nbin write(15,*) 'Nl(',i,')=',Nl(i), & '[',bin(i,1),';',bin(i,2),']' enddo ! i close(15) c record cdf hfile='intensity_cdf.txt' hist_file='./results/' & //hfile(1:strlen(hfile)) open(11,file=hist_file(1:strlen(hist_file))) do i=1,Nbin write(11,*) bin(i,1),cdf(i) write(11,*) bin(i,2),cdf(i) enddo close(11) c record cdf indexes open(11,file='./optimizations/cdf_indexes.txt') do i=1,Nbin write(11,*) i,(cdf_indexes(i,j),j=1,2) enddo ! i close(11) return end