program main

  use vtype_module
  use netcdf
  use ni3_module

  use fi_module

  use set_1d_profile, only : Set1DProfileAtm, Set1DProfileTemp
  use interpolate, only : Interpolate1D

  implicit none

  integer                 :: NcID
  character(extstr)       :: Mode
  character(extstr)       :: Name
  character(extstr)       :: StdName
  character(extstr)       :: LongName
  character(extstr)       :: Units
  integer                 :: NDims
  character(extstr), allocatable :: a_DimNames(:)


  real(dp)                :: TempConst

  real(dp)                :: Ps
  real(dp)                :: Ts

  integer , parameter     :: HITRANNMol = 49
  real(dp)                :: VMR   (HITRANNMol)
  integer                 :: MolNum(HITRANNMol)
!!$  character(8)            :: m_HITRANMolName(HITRANNMol)
!!$  real(sp)                :: m_HITRANMolWt  (HITRANNMol)

!!$  integer                 :: NLevTmp
!!$  real(dp), allocatable   :: r_PressTmp(:)

  integer                 :: NLev
  real(dp), allocatable   :: r_Sigma(:)

  real(dp), allocatable   :: r_Press(:)
  integer                 :: NMol
  integer , allocatable   :: m_MolNum(:)             ! Molecular number
!!$  real(dp), allocatable   :: m_MolWt(:)              ! Molecular weight
!!$  real(dp)                :: m_MolWt(1:3)              ! Molecular weight

  integer                 :: NCldNum

  real(dp), allocatable   :: r_Temp   (:)
  real(dp), allocatable   :: rm_VMR   (:,:)
  real(dp), allocatable   :: rc_CldMMR(:,:)

!!$  real(dp)                :: MolWt

  real(dp), parameter     :: MeanMolWt = 28.966d-3


  character(128) :: ctlfn = 'mkprofile.conf'

  character(128) :: InNcFn
  character(128) :: OutNcFn
  character(128) :: PressName
  character(128) :: TempName
  character(128) :: H2OVapName
  character(128) :: O3Name


  integer              , parameter   :: NChar = stdstr
  integer                            :: a_CharNum( NChar )
  real(8)              , allocatable :: z_Press      (:)
  integer                            :: NPtcl
  integer              , allocatable :: a_PtclNum    (:)
  character(len=stdstr), allocatable :: a_PtclName   (:)
  real(8)              , allocatable :: a_PtclDens   (:)
  real(8)              , allocatable :: za_PtclEffRad(:,:)
  real(8)              , allocatable :: za_PtclMMR   (:,:)


  integer, parameter :: MaxNmlArySize = 256
  real(8)            :: Sigma(MaxNmlArySize)

  integer :: FileUnit
  integer :: iostat

  integer :: k
  integer :: l
  integer :: m
  integer :: n


  namelist /file_nml/ &
    & InNcFn, OutNcFn, &
    & PressName, TempName, H2OVapName, O3Name

  namelist /atmosphere_nml/ &
    & NLev, Sigma, &
    & TempConst, &
    & NMol, MolNum, VMR

  TempConst      = -1.0d0

  InNcFn      = ''
  OutNcFn     = ''
  !
  NLev     = 0
  Sigma(:) = 0.0d0
  !
  NMol        = 0
  MolNum(:)   = 0
!!$  MolNum(1)   = 1
!!$  MolNum(2)   = 2
!!$  MolNum(3)   = 3
  VMR(:)      =  0.0d0
!!$  VMR(1)      = -1.0d0
!!$  VMR(2)      = 300.0d-6
!!$  VMR(3)      = -1.0d0


  call fi_open( CtlFN, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = file_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, file_nml, is something wrong.'
  end if
  write( 6, nml = file_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = atmosphere_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, atmosphere_nml, is something wrong.'
  end if
  write( 6, nml = atmosphere_nml )
  !
  close( FileUnit )


  if ( NMol <= 0 ) then
    stop 'NMol has to be positive.'
  end if
  !
  allocate( m_MolNum(NMol) )
!!$  allocate( m_MolWt (NMol) )
  !
  do m = 1, NMol
    m_MolNum(m) = MolNum(m)
!!$    m_MolWt (m) = m_HITRANMolWt(m_MolNum(m))
  end do


  !
  ! Input
  !

  if ( NLev <= 0 ) then

    Mode = "read"
    call ni3_open( InNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, "level", NLev )
    allocate( r_Press  ( NLev ) )
    allocate( r_Temp   ( NLev ) )
    call ni3_get_var( NcID, PressName , r_Press  )
    call ni3_get_var( NcID, TempName  , r_Temp   )
    call ni3_close( NcID )


    ! Set uniform temperature and mixing ratios
    !
    Ps = r_Press(1)
    if ( TempConst > 0.0d0 ) then
      r_Temp = TempConst
    end if
    Ts = r_Temp (1)

  else

    allocate( r_Sigma ( NLev ) )

!!$    r_Sigma = Sigma(1:NLev)

!!$    do k = 1, NLev-1
    do k = 1, NLev
      r_Sigma(k) = 1.0d0 * exp( - dble(k-1)/5.0d0 )
    end do
!!$    r_Sigma(NLev) = 0.0d0

    ! check
    if ( r_Sigma(1) /= 1.0d0 ) then
      write( 6, * ) 'Sigma(1) has to be 1, but is ', r_Sigma(1), '.'
      stop
    end if

    allocate( r_Press ( NLev ) )
    allocate( r_Temp  ( NLev ) )

!!$    call Set1DProfileAtm(      &
!!$      & InNcFN, PressName, TempName, H2OVapName, O3Name, &
!!$      & NLev, r_Sigma,       &
!!$      & Ps, r_Press, Ts, r_Temp, r_MMRH2O, r_MMRO3       &
!!$      & )
    call Set1DProfileTemp(      &
      & InNcFN, PressName, TempName, &
      & NLev, r_Sigma,       &
      & Ps, r_Press, Ts, r_Temp  &
      & )

    deallocate( r_Sigma )

  end if


  allocate( rm_VMR   ( NLev, NMol) )
  do m = 1, NMol
    m_MolNum(m) = MolNum(m)
    rm_VMR(:,m) = VMR(m)
  end do

  NCldNum = 0
  allocate( rc_CldMMR( NLev, NCldNum) )
  rc_CldMMR = 0.0d0


!!$  IndexH2O = -1
!!$  do m = 1, NMol
!!$    if ( m_MolNum(m) == 1 ) then
!!$      if ( IndexH2O > 0 ) then
!!$        write( 6, * ) 'ERROR: H2O is found twice.'
!!$        stop
!!$      end if
!!$      IndexH2O = m
!!$    end if
!!$  end do
!!$  do m = 1, NMol
!!$    if ( VMR(m) >= 0.0d0 ) then ! H2O
!!$      rm_VMR(:,m) = VMR(m)
!!$!      rm_VMR(:,m) = VMR(m) * ( 1.0d0 - VMR(:,IndexH2O) )
!!$    else
!!$      rm_VMR(:,m) = rm_VMR(:,m)
!!$    end if
!!$  end do


  ! Set particles
  !
  do l = 1, NChar
    a_CharNum(l) = l
  end do
  allocate( z_Press(1:NLev-1) )
  do k = 1, NLev-1
    z_Press(k) = ( r_Press(k) + r_Press(k+1) ) / 2.0d0
  end do

  ! temporarily set
!!$  NPtcl              = 0
!!$  NPtcl              = 2
  NPtcl              = NCldNum
  !
  allocate( a_PtclNum (NPtcl) )
  allocate( a_PtclName(NPtcl) )
  allocate( a_PtclDens(NPtcl) )
  allocate( za_PtclEffRad(NLev-1,NPtcl) )
  allocate( za_PtclMMR   (NLev-1,NPtcl) )
  !
  do n = 1, NPtcl
    a_PtclNum(n) = n
  end do
  a_PtclName         = ""
  a_PtclDens         = 0.0d0
  za_PtclEffRad      = 0.0d0
  za_PtclMMR         = 0.0d0

  do m = 1, NPtcl
    do k = 1, NLev-1
      za_PtclMMR(k,m) = ( rc_CldMMR(k,m) + rc_CldMMR(k+1,m) ) / 2.0d0
    end do
  end do



  ! Output

  call ni3_open( OutNcFN, "new", NcID )

  Name     = "Press"
  StdName  = Name
  LongName = "Pressure"
  Units    = "Pa"
  call ni3_set_dim( NcID, Name, NI3_DOUBLE, r_Press, &
    & StdName, LongName, Units )
  Name     = "MolNum"
  StdName  = Name
  LongName = "Molecular number"
  Units    = "1"
  call ni3_set_dim( NcID, Name, NI3_DOUBLE, m_MolNum, &
    & StdName, LongName, Units )
  Name     = "PtclNum"
  StdName  = Name
  LongName = "Particle number"
  Units    = "1"
  call ni3_set_dim( NcID, Name, NI3_INT, a_PtclNum, &
    & StdName, LongName, Units )
  Name     = "PressM"
  StdName  = Name
  LongName = "Pressure"
  Units    = "Pa"
  call ni3_set_dim( NcID, Name, NI3_DOUBLE, z_Press, &
    & StdName, LongName, Units )
  Name     = "CharNum"
  StdName  = Name
  LongName = "character number"
  Units    = "1"
  call ni3_set_dim( NcID, Name, NI3_INT, a_CharNum, &
    & StdName, LongName, Units )


  NDims = 1
  allocate( a_DimNames( NDims ) )
  !
  Name     = "Temp"
  StdName  = Name
  LongName = "Temperature"
  Units    = "K"
  a_DimNames( 1 ) = "Press"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, r_Temp )
!!$  !
!!$  Name     = "MolWt"
!!$  StdName  = Name
!!$  LongName = "Molecular weight"
!!$  Units    = "kg mol-1"
!!$  a_DimNames( 1 ) = "MolNum"
!!$  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
!!$    & StdName, LongName, Units )
!!$  call ni3_put_var( NcID, Name, m_MolWt )
  !
  deallocate( a_DimNames )


  NDims = 2
  allocate( a_DimNames( NDims ) )
  !
  Name     = "VMR"
  StdName  = Name
  LongName = "Volume mixing ratio"
  Units    = "1"
  a_DimNames( 1 ) = "Press"
  a_DimNames( 2 ) = "MolNum"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, rm_VMR )
  !
  deallocate( a_DimNames )

  !-------------------

  NDims = 0
  allocate( a_DimNames( NDims ) )
  !
  Name     = "Ts"
  StdName  = Name
  LongName = "Surface temperature"
  Units    = "K"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, Ts )
  !
  Name     = "Ps"
  StdName  = Name
  LongName = "Surface pressure"
  Units    = "Pa"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, Ps )
  !
  deallocate( a_DimNames )

  !-------------------

  NDims = 1
  allocate( a_DimNames( NDims ) )
  !
  Name     = "PtclDens"
  StdName  = Name
  LongName = "particle density"
  Units    = "kg m-3"
  a_DimNames( 1 ) = "PtclNum"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, a_PtclDens )
  !
  deallocate( a_DimNames )

  NDims = 2
  allocate( a_DimNames( NDims ) )
  !
  Name     = "PtclName"
  StdName  = Name
  LongName = "particle name"
  Units    = "m"
  a_DimNames( 1 ) = "CharNum"
  a_DimNames( 2 ) = "PtclNum"
  call ni3_def_var( NcID, Name, NI3_CHAR, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, a_PtclName )
  !
  Name     = "PtclEffRad"
  StdName  = Name
  LongName = "particle effective radius"
  Units    = "m"
  a_DimNames( 1 ) = "PressM"
  a_DimNames( 2 ) = "PtclNum"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, za_PtclEffRad )
  !
  Name     = "PtclMMR"
  StdName  = Name
  LongName = "particle mass mixing ratio"
  Units    = "1"
  a_DimNames( 1 ) = "PressM"
  a_DimNames( 2 ) = "PtclNum"
  call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
    & StdName, LongName, Units )
  call ni3_put_var( NcID, Name, za_PtclMMR )
  !
  deallocate( a_DimNames )

  !-------------------

  call ni3_close( NcID )


end program main
