c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine write_kdist_results(i_min,i_max,Nb,Nq,nu_min,nu_max,
     &     w,kdist,delta_k,tr_exact,tr_quad,
     &     pres,temp,Nmol,molec_names,x,filelist)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c
c     Purpose: to write k-distribution output results
c     
c     Inputs:
c       + i_min,i_max: min and max interface indexes
c       + Nb: number of narrowband spectral intervals
c       + Nq: quadrature order
c       + nu_min, nu_max: minimum and maximum values of "nu"
c       + w: array of "Nq" quadrature weights
c       + kdist: output k-distribution data set
c       + delta_k: uncertainties over k-distribution results
c       + tr_exact: mean transmittivities computed from the initial discretized k-spectrum
c       + tr_quad: mean transmittivities computed from the output k-distribution data set
c       + pres: pressure for each level (atm)
c       + temp: temperature for each level (K)
c       + Nmol: number of molecular species for each level
c       + molec_names: name of each molecular species, for each level
c       + x: molar fraction of each species, for each level
c       + filelist: name of the each file that contains the relevant HR spectra
c
      integer Nb,Nq,band,quad,i_min,i_max,i
      double precision nu_min(1:Nbmx),nu_max(1:Nbmx)
      double precision w(Nqmx)
      double precision kdist(1:Nmax,1:Nbmx,1:Nqmx)
      double precision delta_k(1:Nmax,1:Nbmx,1:Nqmx)
      double precision tr_exact(1:Nmax,1:Nbmx)
      double precision tr_quad(1:Nmax,1:Nbmx)
      double precision err(1:Nmax,1:Nbmx)
      double precision pres(1:Nmax)
      double precision temp(1:Nmax)
      integer Nmol(1:Nmax)
      character*(Nchar_mx) molec_names(1:Nmax,1:Nmol_mx)
      double precision x(1:Nmax,1:Nmol_mx)
      character*(Nchar_mx) filelist(1:Nmax)
c     temp
      integer ifile,imol
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine write_kdist_results'
c     Debug
c      double precision ka_moy(1:Nmax,1:Nbmx)
c     Debug

      do i=i_min,i_max
         do band=1,Nb
            if (tr_exact(i,band).gt.1.0D-10) then
               err(i,band)=dabs(tr_exact(i,band)-tr_quad(i,band))
     &              /dabs(tr_exact(i,band))
            else
               err(i,band)=-1.0D+0
            endif
         enddo
      enddo
      
      open(12,file='results/kdist.txt')
      write(12,*) 'K-distribution result file'
      write(12,*) 'Number of HR spectra used:'
      write(12,*) i_max-i_min+1
      do ifile=1,i_max-i_min+1
         write(12,*) '----------------------------------------'
         write(12,*) 'HR spectrum index:',ifile
         write(12,*) filelist(ifile)(1:strlen(filelist(ifile)))
         write(12,*) 'Pressure (atm):'
         write(12,*) pres(ifile)
         write(12,*) 'Temperature (K):'
         write(12,*) temp(ifile)
         write(12,*) 'Number of molecular species:'
         write(12,*) Nmol(ifile)
         do imol=1,Nmol(ifile)
            write(12,*) 'molecular species index: ',imol
            write(12,*) 'name:'
            write(12,*) molec_names(ifile,imol)
     &           (1:strlen(molec_names(ifile,imol)))
            write(12,*) 'molar fraction:'
            write(12,*) x(ifile,imol)
         enddo ! imol
         if (ifile.eq.i_max-i_min+1) then
            write(12,*) '----------------------------------------'
         endif
      enddo ! ifile
      write(12,*)
      write(12,*) 'Number of spectral intervals:'
      write(12,*) Nb
      write(12,*) 'Quadrature order:'
      write(12,*) Nq
      write(12,*) 'Minimum HR spectrum index:'
      write(12,*) i_min
      write(12,*) 'Maximum HR spectrum index:'
      write(12,*) i_max
      write(12,*)
      write(12,*) 'Quadrature weights w(quad):'
      do quad=1,Nq
         write(12,*) w(quad)
      enddo
      write(12,*)
      write(12,*) 'nu_min / nu_max / k(quad),quad=1,Nq / err'
      do i=i_min,i_max
         write(12,*) 'HR spectrum index:',i
         do band=1,Nb
            write(12,50) nu_min(band),nu_max(band),
     &           (kdist(i,band,quad),quad=1,Nq),
     &           err(i,band)
         enddo                  ! band
         write(12,*)
      enddo
      close(12)

      open(14,file='results/delta_kdist.txt')
      write(14,*) 'K-distribution uncertainties result file'
      write(14,*) 'Number of HR spectra used:'
      write(14,*) i_max-i_min+1
      do ifile=1,i_max-i_min+1
         write(14,*) '----------------------------------------'
         write(14,*) 'HR spectrum index:',ifile
         write(14,*) filelist(ifile)(1:strlen(filelist(ifile)))
         write(14,*) 'Pressure (atm):'
         write(14,*) pres(ifile)
         write(14,*) 'Temperature (K):'
         write(14,*) temp(ifile)
         write(14,*) 'Number of molecular species:'
         write(14,*) Nmol(ifile)
         do imol=1,Nmol(ifile)
            write(14,*) 'molecular species index: ',imol
            write(14,*) 'name:'
            write(14,*) molec_names(ifile,imol)
     &           (1:strlen(molec_names(ifile,imol)))
            write(14,*) 'molar fraction:'
            write(14,*) x(ifile,imol)
         enddo ! imol
         if (ifile.eq.i_max-i_min+1) then
            write(14,*) '----------------------------------------'
         endif
      enddo ! ifile
      write(14,*)
      write(14,*) 'Number of spectral intervals:'
      write(14,*) Nb
      write(14,*) 'Quadrature order:'
      write(14,*) Nq
      write(14,*) 'Minimum interface index:'
      write(14,*) i_min
      write(14,*) 'Maximum interface index:'
      write(14,*) i_max
      write(14,*)
      write(14,*) 'Quadrature weights w(quad):'
      do quad=1,Nq
         write(14,*) w(quad)
      enddo
      write(14,*)
      write(14,*) 'nu_min / nu_max / delta(k)(quad),quad=1,Nq'
      do i=i_min,i_max
         write(14,*) 'HR spectrum index:',i
         do band=1,Nb
            write(14,50) nu_min(band),nu_max(band),
     &           (delta_k(i,band,quad),quad=1,Nq)
         enddo                  ! band
         write(14,*)
      enddo
      close(14)

      return
      end



      subroutine write_kdist_results_tmp(i,Nb,Nq,nu_min,nu_max,
     &     w,kdist0,delta_k0,tr_exact0,tr_quad0,
     &     pres,temp,Nmol,molec_names,x,kfile)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c
c     Purpose: to write temporary k-distribution output results
c     
c     Inputs:
c       + i: current (temporary) interface index
c       + Nb: number of narrowband spectral intervals
c       + Nq: quadrature order
c       + nu_min, nu_max: minimum and maximum values of "nu"
c       + w: array of "Nq" quadrature weights
c       + kdist0: output k-distribution data set
c       + delta_k0: uncertainties over k-distribution results
c       + tr_exact0: mean transmittivities computed from the initial discretized k-spectrum
c       + tr_quad0: mean transmittivities computed from the output k-distribution data set
c       + pres: pressure for each level (atm)
c       + temp: temperature for each level (K)
c       + Nmol: number of molecular species for each level
c       + molec_names: name of each molecular species, for each level
c       + x: molar fraction of each species, for each level
c       + kfile: name of the file that contains the relevant HR spectrum
c
c     I/O
      integer i,Nb,Nq
      double precision nu_min(1:Nbmx),nu_max(1:Nbmx)
      double precision w(Nqmx)
      double precision kdist0(1:Nbmx,1:Nqmx)
      double precision delta_k0(1:Nbmx,1:Nqmx)
      double precision tr_exact0(1:Nbmx)
      double precision tr_quad0(1:Nbmx)
      double precision pres(1:Nmax)
      double precision temp(1:Nmax)
      integer Nmol(1:Nmax)
      character*(Nchar_mx) molec_names(1:Nmax,1:Nmol_mx)
      double precision x(1:Nmax,1:Nmol_mx)
      character*(Nchar_mx) kfile
c     temp
      integer band,quad,imol
      double precision err(1:Nbmx)
      character*(Nchar_mx) res_file,kbar_file
      character*3 str3
      double precision kbar(1:Nbmx)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine write_kdist_results_tmp'

      call num2str3(i,str3)

      do band=1,Nb
         if (tr_exact0(band).gt.1.0D-10) then
            err(band)=dabs(tr_exact0(band)-tr_quad0(band))
     &           /dabs(tr_exact0(band))
         else
            err(band)=-1.0D+0
         endif
      enddo ! band
      
      res_file='./results/kdist_'
     &     //str3(1:strlen(str3))
     &     //'.txt'

      open(12,file=res_file(1:strlen(res_file)))
      write(12,*) 'K-distribution result file'
      write(12,*) 'Produced from HR spectrum file:'
      write(12,*) kfile(1:strlen(kfile))
      write(12,*) 'Pressure (atm):'
      write(12,*) pres(i)
      write(12,*) 'Temperature (K):'
      write(12,*) temp(i)
      write(12,*) 'Number of molecular species:'
      write(12,*) Nmol(i)
      do imol=1,Nmol(i)
         write(12,*) 'molecular species index: ',imol
         write(12,*) 'name:'
         write(12,*) molec_names(i,imol)(1:strlen(molec_names(i,imol)))
         write(12,*) 'molar fraction:'
         write(12,*) x(i,imol)
      enddo ! imol
      write(12,*)
      write(12,*) 'Number of spectral intervals:'
      write(12,*) Nb
      write(12,*) 'Quadrature order:'
      write(12,*) Nq
      write(12,*) 'HR spectrum index:'
      write(12,*) i
      write(12,*)
      write(12,*) 'Quadrature weights w(quad):'
      do quad=1,Nq
         write(12,*) w(quad)
      enddo
      write(12,*)
      write(12,*) 'nu_min / nu_max / k(quad),quad=1,Nq / err'
      do band=1,Nb
         write(12,50) nu_min(band),nu_max(band),
     &        (kdist0(band,quad),quad=1,Nq),
     &        err(band)
      enddo                     ! band
      close(12)

      do band=1,Nb
         kbar(band)=0.0D+0
         do quad=1,Nq
            kbar(band)=kbar(band)+kdist0(band,quad)*w(quad)
         enddo ! quad
      enddo ! band
c      kbar_file='./results/kbar_layer'
c     &     //str3(1:strlen(str3))
c     &     //'.txt'
c      open(10,file=kbar_file(1:strlen(kbar_file)))
c      do band=1,Nb
c         write(10,*) nu_min(band),kbar(band)
c         write(10,*) nu_max(band),kbar(band)
c      enddo ! band
c      close(10)

      res_file='./results/delta_kdist_'
     &     //str3(1:strlen(str3))
     &     //'.txt'

      open(14,file=res_file(1:strlen(res_file)))
      write(14,*) 'K-distribution uncertainties result file'
      write(14,*) 'Produced from HR spectrum file:'
      write(14,*) kfile(1:strlen(kfile))
      write(14,*) 'Pressure (atm):'
      write(14,*) pres(i)
      write(14,*) 'Temperature (K):'
      write(14,*) temp(i)
      write(14,*) 'Number of molecular species:'
      write(14,*) Nmol(i)
      do imol=1,Nmol(i)
         write(14,*) 'molecular species index: ',imol
         write(14,*) 'name:'
         write(14,*) molec_names(i,imol)(1:strlen(molec_names(i,imol)))
         write(14,*) 'molar fraction:'
         write(14,*) x(i,imol)
      enddo ! imol
      write(14,*)
      write(14,*) 'Number of spectral intervals:'
      write(14,*) Nb
      write(14,*) 'Quadrature order:'
      write(14,*) Nq
      write(14,*) 'HR spectrum index:'
      write(14,*) i
      write(14,*)
      write(14,*) 'Quadrature weights w(quad):'
      do quad=1,Nq
         write(14,*) w(quad)
      enddo
      write(14,*)
      write(14,*) 'nu_min / nu_max / delta(k)(quad),quad=1,Nq'
      do band=1,Nb
         write(14,50) nu_min(band),nu_max(band),
     &        (delta_k0(band,quad),quad=1,Nq)
      enddo                     ! band
      close(14)

      return
      end