program main

  use vtype_module

  use netcdf
  use ni3_module

  use fi_module

  use rad_planetconst, only : &
    & RadPlanetConstInit, &
    & MolWtNonRadAct , &
    & MolWtH2O

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

  use hitranconst, only : m_HITRANMolWt
  use convection, only : CalcMoistConvDTempDLogPress
  use saturation_simple, only : &
    & CalcPressH2OVapSat

  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(:)


  character(stdstr) :: CtlFN0 = 'planetconst.conf'
  real(DP)          :: GravSet
  real(DP)          :: CpSet
  real(DP)          :: MolWtNonRadActSet

  real(dp)                :: TempConst

  real(dp)                :: Ps
  real(dp)                :: Ts
  real(dp)                :: PTropo
  real(dp)                :: TTropo

  integer , parameter     :: HITRANNMol = 47
  real(dp)                :: VMR(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

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

  real(dp), allocatable   :: r_MMRH2O(:)
  real(dp), allocatable   :: r_MMRO3 (:)

  real(dp)                :: MolWt

  real(dp)                :: MeanMolWt
  real(DP)                :: LatentHeatPerMol
  real(DP), parameter :: LatentHeat     = 2.5d6     ! J kg-1

  real(DP) :: DTempDLogPress
  real(DP) :: VMRNonRadAct

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

  character(128) :: InNcFn
  character(128) :: OutNcFn

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

  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   (:,:)

  logical :: FlagMoistConvection

  integer :: FileUnit
  integer :: iostat

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


  namelist /planetconst_nml/ &
    & GravSet, CpSet, MolWtNonRadActSet

  namelist /file_nml/ &
    & OutNcFn

  namelist /atmosphere_nml/ &
    & Ps, Ts, PTropo, TTropo, &
    & NLev, Sigma, &
    & TempConst, &
    & NMol, VMR

  ! for Earth
  !
  GravSet           = 9.80665d0
  CpSet             = 1004.6d0
  MolWtNonRadActSet = 28.96d-3

  TempConst      = -1.0d0

  InNcFn      = ''
  OutNcFn     = ''
  !
  Ps       = 1000.0d2
  NLev     = 0
  Sigma(:) = 0.0d0
  !
  Ts       = 300.0d0
  !
  PTropo   = -1.0d0
  !
  TTropo   = 200.0d0
  !
  NMol        = 0
  VMR(:)      =  0.0d0
  VMR(1)      = -1.0d0
  VMR(2)      = 300.0d-6
  VMR(3)      = -1.0d0


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


  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 ( NLev <= 0 ) then
    stop 'NLev has to be positive.'
  end if

  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) = m
  end do
!!$  do m = 1, NMol
  do m = 1, min(NMol,3)
    m_MolNum(m) = m
    !
    select case ( m )
    case ( 1 )
      MolWt = 0.0180152636021376
    case ( 2 )
      MolWt = 0.044009767472744
    case ( 3 )
      MolWt = 0.0479970946907997
    case default
      write( 6, * ) 'Unexpected case'
      stop
    end select
    m_MolWt (m) = MolWt
  end do



  ! Initialize planetary constants
  !
  call RadPlanetConstInit( &
    & GravSet           = GravSet,           &
    & CpSet             = CpSet,             &
    & MolWtNonRadActSet = MolWtNonRadActSet  &
    & )

  if ( VMR(1) > 0.0d0 ) then
    FlagMoistConvection = .true.
  else
    FlagMoistConvection = .false.
  end if

  !
  ! Input
  !

  allocate( r_Sigma ( NLev ) )

  if ( all( Sigma == 0.0d0 ) ) then
    k = 1
    r_Sigma(k) = 1.0d0
    do k = 1+1, NLev-2
      r_Sigma(k) = r_Sigma(k-1) * exp( - 1.0d0 / 5.0d0 )
    end do
    k = NLev-1
    r_Sigma(k) = r_Sigma(k-1) / 2.0d0
    k = NLev
    r_Sigma(k) = 0.0d0
  else
    r_Sigma = Sigma(1:NLev)
    ! check
    if ( r_Sigma(1) /= 1.0d0 ) then
      write( 6, * ) 'Sigma(1) has to be 1, but is ', r_Sigma(1), '.'
      stop
    end if
  end if

!!$    r_Sigma(1) = 1.0d0
!!$    do k = 1+1, NLev
!!$      r_Sigma(k) = r_Sigma(k-1) * exp( - 1.0d0 / 20.0d0 )
!!$    end do
  !
!!$  r_Sigma = Sigma(1:NLev)
!!$  ! check
!!$  if ( r_Sigma(1) /= 1.0d0 ) then
!!$    write( 6, * ) 'Sigma(1) has to be 1, but is ', r_Sigma(1), '.'
!!$    stop
!!$  end if
  !
!!$    if ( r_Sigma(NLev) /= 0.0d0 ) then
!!$      write( 6, * ) 'Sigma(Nlev) has to be 0, but is ', r_Sigma(NLev), '.'
!!$      stop
!!$    end if

  allocate( r_Press ( NLev ) )
  allocate( r_Temp  ( NLev ) )
  allocate( r_MMRH2O( NLev ) )
  allocate( r_MMRO3 ( NLev ) )

  r_Press = Ps * r_Sigma

  deallocate( r_Sigma )


  allocate( rm_VMR(NLev, NMol) )
  !
  if ( NMol >= 1 ) then
    if ( VMR(1) >= 0.0d0 ) then ! H2O
      rm_VMR(:,1) = VMR(1)
    end if
  end if
  do m = 2, NMol
    rm_VMR(:,m) = VMR(m)
  end do



  LatentHeatPerMol = LatentHeat * MolWtH2O

  k = 1
  r_Temp(k) = Ts
  if ( VMR(1) >= 0.0d0 ) then ! H2O
    rm_VMR(k,1) = VMR(1)
  else
    rm_VMR(k,1) = CalcPressH2OVapSat( LatentHeatPerMol, r_Temp(k) ) / r_Press(k)
  end if

  if ( PTropo >= 0.0d0 ) then
    ! Tropopause is set by pressure.

    TTropo = -1.0d0

    do k = 1+1, NLev

      VMRNonRadAct = 1.0d0
      do m = 1, NMol
        VMRNonRadAct = VMRNonRadAct - rm_VMR(k-1,m)
      end do
      MeanMolWt = VMRNonRadAct * MolWtNonRadACt
      do m = 1, NMol
        MeanMolWt = MeanMolWt + rm_VMR(k-1,m) * m_HITRANMolWt(m_MolNum(m))
      end do

      call CalcMoistConvDTempDLogPress( &
        & FlagMoistConvection, &
        & r_Press(k-1), MeanMolWt, r_Temp(k-1), &
        & DTempDLogPress &
        & )
      if ( r_Press(k) > 0.0d0 ) then
        r_Temp(k) = r_Temp(k-1) + DTempDLogPress * log( r_Press(k-1) / r_Press(k) )
      else
        write( 6, * ) k, ' Temperature at zero pressure level is the same as that in a lower level.'
        r_Temp(k) = r_Temp(k-1)
      end if

      if ( VMR(1) >= 0.0d0 ) then ! H2O
        rm_VMR(k,1) = VMR(1)
      else
        rm_VMR(k,1) = CalcPressH2OVapSat( LatentHeatPerMol, r_Temp(k) ) / r_Press(k)
      end if
      if ( r_Press(k) < PTropo ) then
        if ( TTropo < 0.0d0 ) then
          TTropo = &
            &   ( r_Temp(k) - r_Temp(k-1) ) &
            & / log( r_Press(k) / r_Press(k-1) ) &
            & * log( PTropo     / r_Press(k-1) ) &
            & + r_Temp(k-1)
        end if
        r_Temp(k)   = TTropo
        if ( VMR(1) >= 0.0d0 ) then ! H2O
          rm_VMR(k,1) = VMR(1)
        else
          rm_VMR(k,1) = CalcPressH2OVapSat( LatentHeatPerMol, r_Temp(k) ) / PTropo
        end if
      end if

    end do

  else
    ! Tropopause is set by temperature

    PTropo = -1.0d0

    do k = 1+1, NLev
      VMRNonRadAct = 1.0d0
      do m = 1, NMol
        VMRNonRadAct = VMRNonRadAct - rm_VMR(k-1,m)
      end do
      MeanMolWt = VMRNonRadAct * MolWtNonRadACt
      do m = 1, NMol
        MeanMolWt = MeanMolWt + rm_VMR(k-1,m) * m_HITRANMolWt(m_MolNum(m))
      end do

      call CalcMoistConvDTempDLogPress( &
        & FlagMoistConvection, &
        & r_Press(k-1), MeanMolWt, r_Temp(k-1), &
        & DTempDLogPress &
        & )
      if ( r_Press(k) > 0.0d0 ) then
        r_Temp(k) = r_Temp(k-1) + DTempDLogPress * log( r_Press(k-1) / r_Press(k) )
      else
        write( 6, * ) k, ' Temperature at zero pressure level is the same as that in a lower level.'
        r_Temp(k) = r_Temp(k-1)
      end if

      if ( VMR(1) >= 0.0d0 ) then ! H2O
        rm_VMR(k,1) = VMR(1)
      else
        rm_VMR(k,1) = CalcPressH2OVapSat( LatentHeatPerMol, r_Temp(k) ) / r_Press(k)
      end if
      if ( r_Temp(k) < TTropo ) then
        if ( PTropo < 0.0d0 ) then
          ! check PTropo
          PTropo = log( r_Press(k) / r_Press(k-1) ) / ( r_Temp(k) - r_Temp(k-1) ) * ( TTropo - r_Temp(k-1) ) + log( r_Press(k-1) )
          PTropo = exp( PTropo )
        end if
        r_Temp(k)   = TTropo
        if ( VMR(1) >= 0.0d0 ) then ! H2O
          rm_VMR(k,1) = VMR(1)
        else
          rm_VMR(k,1) = CalcPressH2OVapSat( LatentHeatPerMol, r_Temp(k) ) / PTropo
        end if
      end if

    end do

  end if


  ! 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
  !
  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

  ! 10 g m-2 cloud water
!!$  a_PtclName(1)      = "water"
!!$  a_PtclName(2)      = "ice"
!!$  a_PtclDens(1)      = 1.0d3
!!$  a_PtclDens(2)      = 1.0d3
!!$  za_PtclEffRad(:,1) = 10.0d-6
!!$  za_PtclEffRad(:,2) = 50.0d-6
!!$  m = 1
!!$  k = 17
!!$  za_PtclMMR(k,m) = &
!!$    & 10.0d-3 &
!!$    &   / ( r_Press(k) - r_Press(k+1) ) / 9.8d0



  ! Output

  Mode = "new"
  call ni3_open( OutNcFN, Mode, NcID )

  Name     = "Press"
  StdName  = Name
  LongName = "Pressure"
  Units    = "Pa"
  call ni3_set_dim( NcID, Name, NF90_DOUBLE, r_Press, &
    & StdName, LongName, Units )
  Name     = "MolNum"
  StdName  = Name
  LongName = "Molecular number"
  Units    = "1"
  call ni3_set_dim( NcID, Name, NF90_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, NF90_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, NF90_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, NF90_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, NF90_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, NF90_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
