program make_composition
      implicit none
c     
c     Purpose: to create input data file "composition_venus.in"
c
      include 'max.inc'
      include 'formats.inc'
      integer m,i,j,nlev_mf,mol,imf,ind,iso
      logical indf
      integer band
      double precision alt(1:Nmax)
      double precision pres(1:Nmax)
      double precision adia(1:Nmax)
      double precision temp(1:Nmax)
      double precision Tsol
      double precision Tesp
      double precision dens_mix(1:Nmax)
      double precision grav(1:Nmax)
      double precision zeta(1:Nmax)
      double precision capa(1:Nmax)
      double precision P_mf(1:Nmax)
      double precision x(1:Nmax,1:Nmol_max)
      double precision nu_lo(1:Nbmx)
      double precision nu_hi(1:Nbmx)
      double precision nu_inf,nu_tr,nu_sup
      double precision dnu_lw,dnu_sw
      integer Nband,Nband_lw,Nband_sw
      character*(Nchar_mx) composition_file,model_file,nu_file
      character*(Nchar_mx) input_file
      character*(Nchar_mx) iso_comp_file
      character*(Nchar_mx) planet_descriptor
c     (P,T) levels
      double precision DeltaT
      integer iprange,itrange
      integer NPrange,NTrange(1:Nrange_mx)
      double precision Prange(1:Nrange_mx)
      double precision Trange_limits(1:Nrange_mx,1:2)
      double precision Trange(1:Nrange_mx,1:Nrange_mx)
      integer Nlev,ilev
      double precision Plev(1:Nmax)
      double precision Tlev(1:Nmax)
      double precision zlev(1:Nmax)
      double precision xlev(1:Nmax,1:Nmol_max)
      double precision x_HDO(1:Nmax)
      double precision abundance_HDO(1:Nmax)
      logical hdo_found,h2o_found
      integer ihdo,ih2o
      logical imol_found
      integer imol
      logical iso161_found,iso162_found
      integer i161,i162
      integer conv_err
c     HITRAN data
      integer HNmol
      integer Hmol_index(1:Nmol_max)
      integer Hmol_niso(1:Nmol_max)
      character*10 Hmol_name(1:Nmol_max)
      integer Hiso_index(1:Nmol_max,1:Niso_max)
      double precision Hiso_abundance(1:Nmol_max,1:Niso_max)
      double precision Hiso_Qref(1:Nmol_max,1:Niso_max)
      integer Hiso_g(1:Nmol_max,1:Niso_max)
      double precision Hiso_mass(1:Nmol_max,1:Niso_max)
c     abundances of isotopologues
      integer Nmol
      integer mol_index(1:Nmol_max)
      integer mol_niso(1:Nmol_max)
      character*10 mol_name(1:Nmol_max)
      integer iso_index(1:Nmol_max,1:Niso_max)
      double precision iso_abundance(1:Nmax,1:Nmol_max,1:Niso_max)
      double precision iso_Qref(1:Nmol_max,1:Niso_max)
      integer iso_g(1:Nmol_max,1:Niso_max)
      double precision iso_mass(1:Nmol_max,1:Niso_max)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='program make_composition'

      planet_descriptor='VENUS 140km'

c     get molecules names and indexes from HITRAN
      input_file='./molparam_hitran2008.txt'
      call read_molparam(input_file,
     &     HNmol,Hmol_index,Hmol_niso,Hmol_name,
     &     Hiso_index,Hiso_abundance,Hiso_Qref,Hiso_g,Hiso_mass)

c     Prange: range of pressure levels for spectra generation
      NPrange=14
      if (NPrange.gt.Nrange_mx) then
         call error(label)
         write(*,*) 'NPrange=',NPrange
         write(*,*) 'while Nrange_mx=',Nrange_mx
         stop
      endif
      Prange(1)=1.0D+7 ! Pa
      do iprange=2,NPrange
         Prange(iprange)=Prange(iprange-1)/10.0D+0 ! Pa
      enddo ! iprange
c     Trange: range of temperature levels for spectra generation
      DeltaT=50.0D+0            ! K
c     Prange(1)=10^7 Pa
      Trange_limits(1,1)=450.0D+0 ! K
      Trange_limits(1,2)=750.0D+0 ! K
c     Prange(2)=10^6 Pa
      Trange_limits(2,1)=300.0D+0 ! K
      Trange_limits(2,2)=750.0D+0 ! K
c     Prange(3)=10^5 Pa
      Trange_limits(3,1)=200.0D+0 ! K
      Trange_limits(3,2)=550.0D+0 ! K
c     Prange(4)=10^4 Pa
      Trange_limits(4,1)=150.0D+0 ! K
      Trange_limits(4,2)=400.0D+0 ! K
c     Prange(5)=10^3 Pa
      Trange_limits(5,1)=150.0D+0 ! K
      Trange_limits(5,2)=300.0D+0 ! K
c     Prange(6)=10^2 Pa
      Trange_limits(6,1)=100.0D+0 ! K
      Trange_limits(6,2)=300.0D+0 ! K
c     Prange(7)=10^1 Pa
      Trange_limits(7,1)=100.0D+0 ! K
      Trange_limits(7,2)=300.0D+0 ! K
c     Prange(8)=10^0 Pa
      Trange_limits(8,1)=100.0D+0 ! K
      Trange_limits(8,2)=300.0D+0 ! K
c     Prange(9)=10^-1 Pa
      Trange_limits(9,1)=100.0D+0 ! K
      Trange_limits(9,2)=300.0D+0 ! K
c     Prange(10)=10^-2 Pa
      Trange_limits(10,1)=100.0D+0 ! K
      Trange_limits(10,2)=300.0D+0 ! K
c     Prange(11)=10^-3 Pa
      Trange_limits(11,1)=100.0D+0 ! K
      Trange_limits(11,2)=300.0D+0 ! K
c     Prange(12)=10^-4 Pa
      Trange_limits(12,1)=100.0D+0 ! K
      Trange_limits(12,2)=300.0D+0 ! K
c     Prange(13)=10^-5 Pa
      Trange_limits(13,1)=100.0D+0 ! K
      Trange_limits(13,2)=300.0D+0 ! K
c     Prange(14)=10^-6 Pa
      Trange_limits(14,1)=100.0D+0 ! K
      Trange_limits(14,2)=300.0D+0 ! K

      do iprange=1,NPrange
         NTrange(iprange)=int((Trange_limits(iprange,2)
     &        -Trange_limits(iprange,1))/DeltaT+1)
         if (NTrange(iprange).gt.Nrange_mx) then
            call error(label)
            write(*,*) 'NTrange(',iprange,')=',NTrange(iprange)
            write(*,*) 'while Nrange_mx=',Nrange_mx
            stop
         endif
         do itrange=1,NTrange(iprange)
            Trange(iprange,itrange)=Trange_limits(iprange,1)
     &           +(itrange-1)*DeltaT
         enddo ! itrange
c     Debug
c         write(*,*) 'Prange(',iprange,')=',Prange(iprange)
c         write(*,*) 'Tmin=',Trange_limits(iprange,1),
c     &        ' Tmax=',Trange_limits(iprange,2),
c     &        ' NTrange(',iprange,')=',NTrange(iprange)
c         write(*,*) (Trange(iprange,itrange),
c     &        itrange=1,NTrange(iprange))
c     Debug
      enddo ! iprange

      Nlev=0
      do iprange=1,NPrange
         Nlev=Nlev+NTrange(iprange)
      enddo ! iprange
      if (Nlev.gt.Nmax) then
         call error(label)
         write(*,*) 'Nlev=',Nlev
         write(*,*) 'while Nmax=',Nmax
         stop
      endif
      ilev=0
      do iprange=1,NPrange
         do itrange=1,NTrange(iprange)
            ilev=ilev+1
            Plev(ilev)=Prange(iprange) ! Pa
            Tlev(ilev)=Trange(iprange,itrange) ! K
            zlev(ilev)=0.0D+0 ! m
         enddo ! itrange
      enddo ! iprange

c     Reading physical model of the atmosphere
      call read_model(m,alt,pres,temp,Tsol,Tesp,
     &     dens_mix,adia,grav,zeta,capa)

      call read_mixratios(Nmol,nlev_mf,P_mf,x,mol_name)
c     2013-10-09

c     Identify HDO molecule index
      hdo_found=.false.
      do mol=1,Nmol
         if (mol_name(mol).eq.'HDO') then
            hdo_found=.true.
            ihdo=mol
            goto 112
         endif
      enddo ! mol
 112  continue
      if (.not.hdo_found) then
         call error(label)
         write(*,*) 'HDO could not be identified among:'
         do mol=1,Nmol
            write(*,*) mol,mol_name(mol)(1:strlen(mol_name(mol)))
         enddo ! mol
      endif

c     Identify H2O molecule index
      h2o_found=.false.
      do mol=1,Nmol
         if (mol_name(mol).eq.'H2O') then
            h2o_found=.true.
            ih2o=mol
            goto 113
         endif
      enddo ! mol
 113  continue
      if (.not.h2o_found) then
         call error(label)
         write(*,*) 'H2O could not be identified among:'
         do mol=1,Nmol
            write(*,*) mol,mol_name(mol)(1:strlen(mol_name(mol)))
         enddo ! mol
      endif

c     Interpolation of species concentrations for every (P,T) level
      do ilev=1,Nlev
         indf=.false.
         do imf=1,nlev_mf-1
            if ((Plev(ilev).gt.P_mf(imf)).and.
     &           (Plev(ilev).lt.P_mf(imf+1))) then
               indf=.true.
               ind=imf
               goto 321
            endif
         enddo ! imf
 321     continue
c     
         if (indf) then
            do mol=1,Nmol
               xlev(ilev,mol)=x(ind,mol)
     &              +(x(ind+1,mol)-x(ind,mol))
     &              /(dlog(P_mf(ind+1))-dlog(P_mf(ind)))
     &              *(dlog(Plev(ilev))-dlog(P_mf(ind)))
            enddo               ! mol
         else                   ! indf=.false.
            if (Plev(ilev).gt.P_mf(nlev_mf)) then
               do mol=1,Nmol
                  xlev(ilev,mol)=x(nlev_mf,mol)
               enddo
            endif
            if (Plev(ilev).lt.P_mf(1)) then
               do mol=1,Nmol
                  xlev(ilev,mol)=x(1,mol)
               enddo            ! mol
            endif               ! Plev(ilev)<P_mf(1)
         endif                  ! indf
      enddo ! ilev

c     Add HDO content to H2O
      do ilev=1,Nlev
         abundance_HDO(ilev)=xlev(ilev,ihdo)
     &        /(xlev(ilev,ih2o)+xlev(ilev,ihdo))
         xlev(ilev,ih2o)=xlev(ilev,ih2o)+xlev(ilev,ihdo)
      enddo ! ilev
c     Remove HDO from the list of main molecules and concentration profiles
      do mol=ihdo,Nmol-1
         mol_name(mol)=mol_name(mol+1)
         do ilev=1,Nlev
            xlev(ilev,mol)=xlev(ilev,mol+1)
         enddo                  ! ilev
      enddo                     ! mol
      Nmol=Nmol-1

      do ilev=1,Nlev
c     Modification of values obtained from the physical model
         do mol=1,Nmol
c     H2O
            if (mol_name(mol)(1:strlen(mol_name(mol)))
     &           .eq.'H2O') then
               if (Plev(ilev).lt.1.0D+0) then
                  xlev(ilev,mol)=1.5D-6 ! 1.5 ppm
               endif
            endif ! H2O
c     SO2
            if (mol_name(mol)(1:strlen(mol_name(mol)))
     &           .eq.'SO2') then
               if ((Plev(ilev).le.4.0D+3).and.
     &              (Plev(ilev).gt.1.0D+2)) then
                  xlev(ilev,mol)=100.0D-9
     &                 -(100.0D-9-10.0D-9)
     &                 /(dlog(4.0D+3)-dlog(1.0D+2))
     &                 *(dlog(4.0D+3)-dlog(Plev(ilev))) ! linear with log(P) between 100 ppb @ P=4000 Pa and 10 ppb @ P=100 Pa
               else if (Plev(ilev).le.1.0D+2) then
                  xlev(ilev,mol)=0.3D-6 ! 0.3 ppm
               endif
            endif ! SO2
c     SO
c            if (mol_name(mol)(1:strlen(mol_name(mol)))
c     &           .eq.'SO') then
c               if ((Plev(ilev).le.1.0D+1).and.
c     &              (Plev(ilev).ge.1.0D-6)) then
c                  xlev(ilev,mol)=0.3D-6 ! 0.3 ppm between P=10 Pa and P=10^-6 Pa
c               else
c                  xlev(ilev,mol)=0.0D+0 ! absent
c               endif
c            endif ! SO
c
         enddo ! mol
      enddo                     ! ilev
c
c     =================================================================================================================
c     Composition file
c
      composition_file='./composition.in'
      call write_composition_file(composition_file,
     &     planet_descriptor,Nlev,Nmol,mol_name,
     &     zlev,Plev,Tlev,xlev)
c
c     =================================================================================================================
c     Isotopic composition file
c      
      do mol=1,Nmol
c     Set mol_index
         mol_index(mol)=mol
c     Special case of H2O
         if (mol_name(mol)(1:strlen(mol_name(mol))).eq.'H2O') then
c     Identify H2O within HITRAN molecules
            h2o_found=.false.
            do i=1,HNmol
               if (Hmol_name(i)(1:strlen(Hmol_name(i))).eq.'H2O') then
                  h2o_found=.true.
                  ih2o=i
                  goto 211
               endif
            enddo ! i
 211        continue
            if (.not.h2o_found) then
               call error(label)
               write(*,*) 'H2O could not be found in the '
     &              //'list of HITRAN molecules'
               stop
            endif
c     Leave only main isotope (161) and HDO (162)
            iso161_found=.false.
            do j=1,Hmol_niso(ih2o)
               if (Hiso_index(ih2o,j).eq.161) then
                  iso161_found=.true.
                  i161=j
                  goto 212
               endif
            enddo ! j
 212        continue
            if (.not.iso161_found) then
               call error(label)
               write(*,*) 'isotop 161 could not be found among:'
               do j=1,Hmol_niso(ih2o)
                  write(*,*) 'Isotope ',j,' has index:',
     &                 Hiso_index(ih2o,j)
               enddo ! j
               stop
            endif
c     Debug
c            write(*,*) 'i161=',i161
c     Debug
            iso162_found=.false.
            do j=1,Hmol_niso(ih2o)
               if (Hiso_index(ih2o,j).eq.162) then
                  iso162_found=.true.
                  i162=j
                  goto 213
               endif
            enddo ! j
 213        continue
            if (.not.iso162_found) then
               call error(label)
               write(*,*) 'isotop 162 could not be found among:'
               do j=1,Hmol_niso(ih2o)
                  write(*,*) 'Isotope ',j,' has index:',
     &                 Hiso_index(ih2o,j)
               enddo ! j
               stop
            endif
c     Debug
c            write(*,*) 'i162=',i162
c     Debug
c     Set mol_niso
            mol_niso(mol)=2
c     Set iso_index
            iso_index(mol,1)=161
            iso_index(mol,2)=162
c     Set iso_Qref
            iso_Qref(mol,1)=Hiso_Qref(ih2o,i161)
            iso_Qref(mol,2)=Hiso_Qref(ih2o,i162)
c     Set iso_g
            iso_g(mol,1)=Hiso_g(ih2o,i161)
            iso_g(mol,2)=Hiso_g(ih2o,i162)
c     Set iso_mass
            iso_mass(mol,1)=Hiso_mass(ih2o,i161)
            iso_mass(mol,2)=Hiso_mass(ih2o,i162)
c     Set iso_abundance
            do ilev=1,Nlev
               iso_abundance(ilev,mol,1)=1.0D+0-abundance_HDO(ilev)
               iso_abundance(ilev,mol,2)=abundance_HDO(ilev)
            enddo ! ilev
               
c     Other molecules
         else
c     Indentify molecule index 'mol' among HITRAN molecules
            imol_found=.false.
            do i=1,HNmol
               if (mol_name(mol)(1:strlen(mol_name(mol)))
     &              .eq.Hmol_name(i)(1:strlen(Hmol_name(i)))) then
                  imol_found=.true.
                  imol=i
                  goto 311
               endif
            enddo ! i
 311        continue
c     Set mol_niso, iso_index, iso_Qref, iso_g, iso_mass, iso_abundance
            mol_niso(mol)=Hmol_niso(imol)
            do j=1,mol_niso(mol)
               iso_index(mol,j)=Hiso_index(imol,j)
               iso_Qref(mol,j)=Hiso_Qref(imol,j)
               iso_g(mol,j)=Hiso_g(imol,j)
               iso_mass(mol,j)=Hiso_mass(imol,j)
               do ilev=1,Nlev
                  iso_abundance(ilev,mol,j)=Hiso_abundance(imol,j)
               enddo ! ilev
            enddo ! iso
         endif ! mol_name(mol)
      enddo ! mol

      iso_comp_file='./molparam.in'
      call write_isotopic_composition_file(iso_comp_file,
     &     Nlev,Nmol,mol_name,mol_index,mol_niso,iso_index,
     &     iso_abundance,iso_Qref,iso_g,iso_mass)
c
c     =================================================================================================================
c     Narrowbands file
c
      nu_inf=2.0D+1 ! inv. cm (lower limit of LW)
      nu_tr=11800.0D+0 ! inv. cm (transition between LW and SW)
      nu_sup=5.0D+4 ! inv. cm (higher limit of SW)
      dnu_lw=20.0D+0 ! inv. cm (LW band width)
      dnu_sw=200.0D+0 ! inv. cm (SW band width)

      band=0
      do while (nu_hi(band).lt.nu_tr)
         band=band+1
         nu_lo(band)=nu_inf+dnu_lw*(band-1)
         nu_hi(band)=nu_inf+dnu_lw*(band)
      enddo ! while
      Nband_lw=band
      do while (nu_hi(band).lt.nu_sup)
         band=band+1
         nu_lo(band)=nu_tr+dnu_sw*(band-Nband_lw-1)
         nu_hi(band)=nu_tr+dnu_sw*(band-Nband_lw)
      enddo ! while
      Nband=band
      Nband_sw=Nband-Nband_lw
      if (Nband.gt.Nbmx) then
         call error(label)
         write(*,*) 'Nband=',Nband
         write(*,*) 'while Nbmx=',Nbmx
         stop
      endif

      nu_file='./narrowbands.in'
      call write_narrowbands_file(nu_file,
     &     Nband,nu_lo,nu_hi)
      
      end