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 find_min(tab,n,min,imin)
      implicit none
      include 'max.inc'
      
      double precision tab(1:Nmax_tab)
      double precision min
      integer n,imin,i

      min=tab(1)
      imin=1
      do i=1,n
         if (tab(i).lt.min) then
            min=tab(i)
            imin=i
         endif
      enddo

      return
      end



      subroutine find_max(tab,n,max,imax)
      implicit none
      include 'max.inc'
      
      double precision tab(1:Nmax_tab)
      double precision max
      integer n,imax,i

      max=tab(1)
      imax=1
      do i=1,n
         if (tab(i).gt.max) then
            max=tab(i)
            imax=i
         endif
      enddo

      return
      end

      subroutine find_max_vt(tab,n,max,imax)
      implicit none
      include 'max.inc'
      
      double precision tab(1:Nmax_vt)
      double precision max
      integer n,imax,i

      max=tab(1)
      imax=1
      do i=1,n
         if (tab(i).gt.max) then
            max=tab(i)
            imax=i
         endif
      enddo

      return
      end



      subroutine test_nan(value,out)
      implicit none
      include 'max.inc'
      integer t1,t2,out
      double precision value
c
c     Purpose: to test whether "value" is a NaN or not
c
c     Input: value (double precision)
c     
c     Output: out
c     if out=0, "value" has a real value
c     if out=1, "value" is a NaN
c
      out=0 ! default
      if (value.le.0.0D+0) then
         t1=0
      else
         t1=1
      endif
      if (value.ge.0.0D+0) then
         t2=0
      else
         t2=1
      endif
      if ((t1.eq.1).and.(t2.eq.1)) then
         out=1
      endif

      return
      end


      subroutine test_inf(value,out)
      implicit none
      include 'max.inc'
      double precision value,huge,inf
      double precision d1mach
      integer out
c
c     Purpose: to test whether "value" is INF or not
c
c     Input: value (double precision)
c
c     Output: out
c       + out=0 means "value" has a finite value
c       + out=1 means "value" is infinite
c

      huge=d1mach(2) ! is the maximum amplitude
      out=0 ! default=no problem
      if (dabs(value).ge.huge) then
         out=1 ! value=Inf
      endif

      return
      end


      function frac(x)
      implicit none
      include 'max.inc'
      double precision frac,x
      
      frac=x-int(x)
      
      return
      end


      subroutine find_free_process(free,nproc,fpf,proc)
      implicit none
c
c     Purpose: to find 'proc', the index of a free child process
c
c     Input:
c       + free: 0/1 value for each child process (0: busy process, 1: free process)
c       + nproc: number of child processes
c
c     Output:
c       + fpf: 0 if no free child process found; 1 otherwise
c       + proc: index of the free child process found, if fpf=1
c
      include 'max.inc'
      integer free(1:Nproc_mx)
      integer proc,nproc,fpf

      fpf=0
      do proc=1,nproc
         if (free(proc).eq.1) then
            fpf=1
            goto 123
         endif
      enddo
 123  continue

      return
      end



      subroutine stopped_processes(free,nproc,nsp)
      implicit none
c
c     Purpose: to identify the number of child processes that have been stopped
c
c     Input:
c       + free: array of free processes (a value of 2 means that the corresponding child process has been stopped)
c       + nproc: number of child processes
c     
c     Output:
c       + nsp: number of child processes that have been stopped
c     
      include 'max.inc'
      integer free(1:Nproc_mx)
      integer proc,nproc,nsp
      nsp=0
      do proc=1,nproc
         if (free(proc).eq.2) then
            nsp=nsp+1
         endif
      enddo

      return
      end


      subroutine stringp_file_name(pindex,string,stringp_file)
      implicit none
      include 'max.inc'
      include 'formats.inc'
      integer strlen,pindex
      character*(Nchar_mx) string
      character*3 pindex_ch3
      character*(Nchar_mx) stringp_file

      call num2str3(pindex,pindex_ch3)
      stringp_file=string(1:strlen(string))//
     &     '_'//
     &     pindex_ch3(1:strlen(pindex_ch3))

      return
      end


      subroutine set_chunks(nchunks,np,nk,ik1,ik2,Nkrun)
      implicit none
      include 'max.inc'
      integer nchunks,np,nk
      integer Nkpp,Nklp,chunk
      integer Nkrun(1:Nchunk_mx)
      integer ik1(1:Nchunk_mx)
      integer ik2(1:Nchunk_mx)


      if (nchunks.eq.-1) then
         nchunks=np-1
      endif
      if (nk.ge.nchunks) then
         Nkpp=nk/nchunks
         Nklp=nk-Nkpp*(nchunks-1)
         
         do chunk=1,nchunks
            ik1(chunk)=(chunk-1)*Nkpp+1
            if (chunk.lt.nchunks) then
               ik2(chunk)=chunk*Nkpp
               Nkrun(chunk)=Nkpp
            else
               ik2(chunk)=nk
               Nkrun(chunk)=Nklp
            endif
         enddo                  ! chunk

      else
         do chunk=1,nk
            ik1(chunk)=chunk
            ik2(chunk)=ik1(chunk)
            Nkrun(chunk)=1
         enddo
         do chunk=nk+1,nchunks
            ik1(chunk)=1
            ik2(chunk)=ik1(chunk)
            Nkrun(chunk)=0
         enddo
      endif

      return
      end


      subroutine get_date(datep_file,date)
      implicit none
      include 'max.inc'
      integer strlen
      integer*8 date
      character*(Nchar_mx) datep_file
      character*(Nchar_mx) command

      command="echo $(date +%s) > "//
     &     datep_file(1:strlen(datep_file))
      call exec(command)
      open(11,file=datep_file(1:strlen(datep_file)))
      read(11,*) date
      close(11)

      return
      end


      subroutine set_limits(over,rid,Nb,m,i_inf,i_sup,band_inf,
     &     band_sup,uaal,uanb,last_i,last_band,i_min,i_max,
     &     band_min,band_max,dbr,lastf,lastl,pass,raz)
      implicit none
      include 'max.inc'
      integer over
      logical rid
      integer Nb,m,last_i,last_band,raz
      integer i_min,i_max,band_min,band_max
      integer uaal,uanb
      integer i_inf
      integer i_sup
      integer band_inf
      integer band_sup
      integer dbr,lastf,lastl,pass

      if ((over.eq.0).and.(rid)) then
c     special case of an interruption before anything could actually be done
         if ((last_band.eq.0).and.
     &        (last_i.eq.0)) then
            raz=1
         else
            raz=0
         endif
      else
         raz=1
      endif
           
      if (raz.eq.0) then ! computation has to be resumed
         if (last_band.eq.Nb) then
            i_min=last_i+1
            band_min=1
         else
            i_min=last_i
            band_min=last_band+1
         endif
         if (uaal.eq.1) then
            i_max=m
         else
            i_max=i_sup
         endif
         if (uanb.eq.1) then
            band_max=Nb
         else
            band_max=band_sup
         endif
c         raz=0
      else ! start a new computation
         if (uaal.eq.1) then
            i_min=1
            i_max=m
         else
            i_min=i_inf
            i_max=i_sup
         endif
         if (uanb.eq.1) then
            band_min=1
            band_max=Nb
         else
            band_min=band_inf
            band_max=band_sup
         endif
         lastf=1
         lastl=0
         pass=1
c         raz=1
      endif

      return
      end



      subroutine write_over(ind)
      implicit none
      include 'max.inc'
      integer ind

      open(10,file='./over')
      write(10,*) ind
      close(10)

      return
      end



      subroutine write_backup(i,band,dbr,lastf,lastl,pass)
      implicit none
      include 'max.inc'
      integer i,band,dbr,lastf,lastl,pass

      open(10,file='./backup.txt')
      write(10,*) i
      write(10,*) band
      write(10,*) dbr
      write(10,*) lastf
      write(10,*) lastl
      write(10,*) pass
      close(10)

      return
      end



      subroutine print_cputime(cputime)
      implicit none
      include 'max.inc'
      include 'formats.inc'
      integer*8 cputime
      integer day,hour,min,sec

      call convert_time(cputime,day,hour,min,sec)
      if (day.ge.1) then
         write(*,47) day,'days',hour,'hours',min,'min',sec,'sec'
      else if (hour.ge.1) then
         write(*,46) hour,'hour',min,'min',sec,'sec'
      else if (min.ge.1) then
         write(*,45) min,'min',sec,'sec'
      else
         write(*,44) sec,'sec'
      endif
         

      return
      end




      subroutine convert_time(cputime,day,hour,min,sec)
      implicit none
      include 'max.inc'
      integer*8 cputime
      integer day,hour,min,sec
      integer*8 s
      integer nspd,nsph,nspm
      day=0
      hour=0
      min=0
      sec=0
      s=cputime
      nspm=60
      nsph=60*nspm
      nspd=24*nsph
      
      do while (s.ge.nspd)
         day=day+1
         s=s-nspd
      enddo
      do while (s.ge.nsph)
         hour=hour+1
         s=s-nsph
      enddo
      do while (s.ge.nspm)
         min=min+1
         s=s-nspm
      enddo
      sec=s
      
      return
      end



      subroutine rmfile(tfile)
      implicit none
      include 'max.inc'
c     Purpose: to erase the specified file
      integer strlen
      character*(Nchar_mx) tfile,path
      character*(Nchar_mx) command

      path='./'
      command='rm -rf '//path(1:strlen(path))//tfile(1:strlen(tfile))
      call exec(command)

      return
      end


      subroutine exec(command)
      implicit none
      include 'max.inc'
      integer strlen
      character*(Nchar_mx) command

      call system(command(1:strlen(command)))
c     Debug
c      write(*,*) 'command was executed:'
c      write(*,*) command(1:strlen(command))
c     Debug

      return
      end



      subroutine set_Nk_total(raz,i,pass,uanb,Nb,band_inf,
     &     band_min,band_sup,Nk_total,Nk_remain)
      implicit none
      include 'max.inc'
      integer i,pass,uanb,band_inf,band_min,band_sup,band,bs,be
      integer raz,Nk_total,nl,Nb,Nk_remain
      character*(Nchar_mx) kfile,linesp_file
      
      linesp_file='./nlines'

      if (raz.eq.1) then
         Nk_total=0
      else
         if (uanb.eq.1) then
            bs=1
         else
            bs=band_inf
         endif
         if (pass.eq.1) then
            be=band_min-1
         else
            if (uanb.eq.1) then
               be=Nb
            else
               be=band_sup
            endif
         endif
         if (be.ge.bs) then
            do band=bs,be
               call result_file_name_tmp(i,band,kfile)
               call get_nlines(kfile,linesp_file,nl)
               Nk_total=Nk_total+nl
            enddo               ! band
         endif
      endif
      Nk_remain=0

      return
      end



      subroutine calc_Nk_remain(i,pass,Npass,band_min,band_max,
     &     Nk_total,Nk_remain)
      implicit none
      include 'max.inc'
      integer i,pass,Npass,band_min,band_max,Nk_total
      integer nl,Nk_remain,band,bs,be
      character*(Nchar_mx) kfile,linesp_file

      
      linesp_file='./nlines'
      Nk_remain=0

      if (pass.gt.1) then
         bs=band_min
         be=band_max
         
         if (be.ge.bs) then
            do band=bs,be
               call result_file_name_tmp(i,band,kfile)
               call get_nlines(kfile,linesp_file,nl)
               Nk_remain=Nk_remain+nl
            enddo               ! band
         endif
         Nk_remain=Nk_remain+(Npass-pass)*Nk_total
      endif

      return
      end



      subroutine set_Sref_frac(sda,p,Sref_frac)
      implicit none
      include 'max.inc'
      double precision p,Sref_frac
      integer sda
c     input:
c       + sda: spectral discretization algorithm choice (1/2)
c       + p [atm]
c     output:
c       + Sref_frac

c      if (p.gt.1.0D-1) then
c         Sref_frac=200.0D+0
c      else if (p.gt.1.0D-2) then
c         Sref_frac=130.0D+0
c      else if (p.gt.1.0D-3) then
c         Sref_frac=100.0D+0
c      else if (p.gt.1.0D-4) then
c         Sref_frac=80.0D+0
c      else if (p.gt.1.0D-5) then
c         Sref_frac=40.0D+0
c      else if (p.gt.1.0D-6) then
c         Sref_frac=30.0D+0
c     else
      if (sda.eq.1) then
         Sref_frac=2.0D+3
      else
         Sref_frac=1.0D+8
      endif
c      endif

      return
      end




      subroutine choice_bslal(bslal,dfrac,gfrac,dnufrac,
     &     pmin,dfact,gfact)
      implicit none
      include 'max.inc'
      integer bslal
      double precision dfrac,gfrac,dnufrac
      double precision pmin,dfact,gfact
      
      if (bslal.eq.1) then
         dfrac=3.0D+0
         gfrac=0.65D+0
         dnufrac=10.0D+0
      else if (bslal.eq.2) then
         dfrac=4.0D+0
         gfrac=0.6D+0
         dnufrac=20.0D+0
      else if (bslal.eq.3) then
         dfrac=5.0D+0
         gfrac=0.55D+0
         dnufrac=30.0D+0
      else if (bslal.eq.4) then
         dfrac=6.0D+0
         gfrac=0.5D+0
         dnufrac=40.0D+0
      endif
      pmin=1.0D-3 ! atm
      dfact=20.0D+0
      gfact=10.0D+0

      return
      end



      subroutine fast_discretization(nu1,nu2,dnu,nk,nu)
      implicit none
      include 'max.inc'
c
c     Purpose: to generate a uniform spectral discretization between
c     two values of wavenumber
c
c     Inputs:
c       + nu1: lower wavenumber (cm¯¹)
c       + nu2: upper wavenumber (cm¯¹)
c       + dnu: (fixed) spectral resolution (cm¯¹)
c     
c     Outputs:
c       + nk: number of spectral points
c       + nu: spectral mesh
c

c     I/O
      double precision nu1,nu2,dnu
      integer nk
      double precision nu(1:Nkmx)
c     temp
      integer npt,i
      integer strlen
      character*(Nchar_mx) label
      label='routine fast_discretization'

c     Tests

c     special case of lines that have the same central wavenumbers
      if (nu2.eq.nu1) then
         goto 666 ! do nothing
      endif
c     special case of lines which center moved so much due to pressure
c     that the spectral domain to discretize is negative
      if ((nu2.le.nu1).or.(dnu.le.0.0D+0)) then
         goto 666 ! do nothing
      endif
c     Tests

c     (nu2-nu1) may not be a multiple of dnu; this is why I have to compute npt,
c     the number of points (as an integer):
      npt=int((nu2-nu1)/dnu)+1
c     Debug
c      write(*,*) 'nu1=',nu1
c      write(*,*) 'nu2=',nu2
c      write(*,*) 'dnu=',dnu
c      write(*,*) 'npt=',npt
c     Debug

c     Tests
      if (npt.le.0) then
         write(*,*) 'Error from ',label(1:strlen(label)),' :'
         write(*,*) 'npt=',npt
         write(*,*) 'nu1=',nu1
         write(*,*) 'nu2=',nu2
         write(*,*) 'nu2-nu1=',nu2-nu1
         write(*,*) 'dnu=',dnu,(nu2-nu1)/dnu
         stop
      endif
      if (npt.eq.1) then
         npt=2
      endif
c     Tests

c     uniform spectral discretization
      do i=2,npt
         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+(nu2-nu1)/dble(npt-1)*dble(i-1)
      enddo

 666  continue
      return
      end



      subroutine param1(nu0,P,delta,gamma_L,gamma_D,nuc,gamma)
      implicit none
      include 'max.inc'
      double precision nu0,P,delta,gamma_L,gamma_D,nuc,gamma

      nuc=nu0+delta*P
c     "nuc" is the line's central wavenumber
      gamma=gamma_L
      if (gamma_D.gt.gamma) then
         gamma=gamma_D
      endif
c     "gamma" is the highest value between "gamma_L" and "gamma_D"
      
      return
      end


      subroutine choose_nbw(nu,nbw)
      implicit none
      include 'max.inc'
      double precision nu,nbw

      if (nu.lt.5.0D+3) then
         nbw=20.0D+0
      else if (nu.lt.10.0D+3) then
         nbw=100.0D+0
      else if (nu.lt.12.5D+3) then
         nbw=150.0D+0
      else if (nu.lt.25.0D+3) then
         nbw=200.0D+0
      else
         nbw=500.0D+0
      endif
c     Debug
c      write(*,*) 'nu=',nu,' nbw=',nbw
c     Debug
      
      return
      end



      subroutine choose_intw(nu_min,nu_max,large_db,intw)
      implicit none
      include 'max.inc'
c
c     Purpose: to return "intw", the wavenumber width of the small spectral interval
c     that has to be discretized by a given process.
c
c     I/O
      double precision nu_min,nu_max,intw
      logical large_db
c     temp
      double precision nu,dnu
      integer nlimits
      double precision limits(1:100)
      integer strlen
      character*(Nchar_mx) label
      label='routine choose_intw'

      call provide_nulimits(nlimits,limits)

      nu=(nu_min+nu_max)/2.0D+0
      dnu=nu_max-nu_min
      if (dnu.le.0.0D+0) then
         write(*,*) 'Error from ',label(1:strlen(label)),' :'
         write(*,*) 'nu_max=',nu_max
         write(*,*) 'is lower or equal to:'
         write(*,*) 'nu_min=',nu_min
         stop
      endif
         
      if (large_db) then
         intw=1.0D-1
      else
         if (nu.lt.limits(1)) then
            intw=1.0D+0
         else if (nu.lt.limits(2)) then
            intw=10.0D+0
         else
            intw=100.0D+0
         endif
      endif


      if (intw.gt.dnu) then
         intw=dnu
      endif

      return
      end



      subroutine provide_nulimits(nlimits,limits)
      implicit none
      include 'max.inc'
c
c     Purpose: to provide the wavenumber limits where narrowband width changes
c     These limits are used by routine "choose_intw" that returns the spectral width
c     of the interval that will be discretized by a given process;
c     but it also has to be used by routine "band_limits" that computes the spectral
c     narrowband widths, so that small intervals returned by "choose_intw" fit
c     into narrowband intervals.
c
      integer nlimits
      double precision limits(1:100)

      nlimits=2
      limits(1)=6.0D+3
      limits(2)=2.0D+4

      return
      end



      subroutine line_intervals(db,line_index,Nbin,bin,cdf_indexes,
     &     k,max_ep,epsilon,
     &     current_bin,trigger,
     &     k_history,sum_lines_contributions)
      implicit none
      include 'max.inc'
c
c     Purpose: optimize the convergent calculation of k(nu)
c
c     Input:
c       + Nbin: number of S bins
c       + line_index: index of the current line
c       + bin: limits of the "Nbin" S bins
c       + cdf_indexes: line indexes for each bin
c       + k: current value of k
c       + max_ep: maximum error percentage (user-defined) over k
c       + epsilon: used to switch algorithms
c       + current_bin: index of the bin the current line is located in
c       + trigger: line index that triggers the call to this routine
c       + k_history: values of k at the end of each line bin
c
c     Output:
c       + current_bin, trigger, k_history: updated
c       + sum_lines_contributions: switch that stops the line contribution adding process
c
c     I/O
      logical db
      integer line_index
      integer Nbin
      double precision bin(1:Nbin_mx,1:2)
      integer cdf_indexes(1:Nbin_mx,1:2)
      double precision k
      double precision max_ep
      double precision epsilon
      integer current_bin
      integer trigger
      double precision k_history(0:Nbin_mx)
      logical sum_lines_contributions
c     temp
      integer i
      double precision evolution(1:Nbin_mx)
      double precision last_pc,max_pc
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine line_intervals'

c     update history of k for that nu
      k_history(current_bin)=k

c     compute progression of k-computation
      if (k_history(current_bin).gt.0.0D+0) then
         do i=1,current_bin
c     percentage of current value of k reached at each step
            evolution(i)=k_history(i)/k_history(current_bin)
         enddo ! i
         max_pc=0.0D+0
         do i=2,current_bin
            if ((evolution(i)-evolution(i-1)).gt.max_pc) then
c     percentage of variation between two successive steps -> find highest value
               max_pc=evolution(i)-evolution(i-1)
            endif
         enddo ! i
c     last variation
         if (current_bin.gt.1) then
            last_pc=evolution(current_bin)-evolution(current_bin-1)
         else
            last_pc=1.0D+0
         endif
c     stop computation if last variation is negligible compared to the
c     maximum variation observed from the beginning
         if (last_pc.lt.epsilon*max_ep*max_pc) then
            sum_lines_contributions=.false.
         endif
c     not really necessary (default behavior)
      else
c     If k is still null, keep running !
         sum_lines_contributions=.true.
      endif

c     also stop computation when last line has been reached
      if (current_bin.eq.Nbin) then
         sum_lines_contributions=.false.
      else
         current_bin=current_bin+1
         trigger=cdf_indexes(current_bin,2)
      endif

c     Debug
c      sum_lines_contributions=.false.
c     Debug

      return
      end