program main

  use vtype_module
  use netcdf
  use ni3_module

  use fi_module

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

  use io_profile, only : IOProfileOutput

  use constants0

  implicit none

  integer, parameter :: IDPtclNone  = 0
  integer, parameter :: IDPtclEarth = 1
  integer, parameter :: IDPtclMars  = 2

  integer            :: IDPtcl      = IDPtclNone

  integer                 :: NcID
  character(extstr)       :: Mode


  real(dp)                :: TempConst

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

  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), 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
  real(8)             , allocatable :: z_Press      (:)
  integer                           :: NPtcl
  integer             , allocatable :: a_PtclNum    (:)
  character(len=NChar), allocatable :: a_PtclName   (:)
  real(8)             , allocatable :: a_PtclDens   (:)
  real(8)             , allocatable :: za_PtclEffRad(:,:)
  real(8)             , allocatable :: za_PtclMMR   (:,:)

  ! variables for Mars
  real(dp) :: Grav
  real(dp) :: DustMMR0
  real(dp) :: DustVerDistCoef
  real(dp) :: DustMaxHeight
  real(dp) :: DustVerDistRefPress
  real(dp) :: DustEffRadius
  real(dp) :: DustDens
  real(dp) :: DustExtEffFact067
  real(dp) :: DODFactor
  real(dp) :: DustColMass
  real(dp) :: DustColNumDen

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

  integer :: FileUnit
  integer :: iostat

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


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

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

  TempConst      = -1.0d0

  InNcFn      = ''
  OutNcFn     = ''
  !
  NLev     = 0
  Sigma(:) = 0.0d0
  !
  NMol        = 0
  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 'NMolNum has to be positive.'
  end if
  allocate( m_MolNum(NMol) )
!!$  allocate( m_MolWt (NMolNum) )
  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


  !
  ! Input
  !

  if ( NLev <= 0 ) then

    Mode = "read"
    call ni3_open( InNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, PressName, NLev )

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

    call ni3_get_var( NcID, PressName , r_Press  )
    call ni3_get_var( NcID, TempName  , r_Temp   )
    call ni3_get_var( NcID, H2OVapName, r_MMRH2O )
    call ni3_get_var( NcID, O3Name    , r_MMRO3  )

    call ni3_close( NcID )


!!$    do k = 1, NLev
!!$!      r_Temp(k) = 220.0d0 - 20.0d0 * log10(r_Press(k)/r_Press(1))
!!$      r_Temp(k) = 280.0d0 + 20.0d0 * log10(r_Press(k)/r_Press(1))
!!$    end do

!!$    do k = 1+1, NLev
!!$      if ( r_Temp(k) > r_Temp(k-1) ) exit
!!$    end do
!!$    r_Temp(k:NLev) = r_Temp(k-1)



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

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

    deallocate( r_Sigma )

  end if


  allocate( rm_VMR(NLev, NMol) )
  !
  if ( NMol >= 1 ) then
    if ( VMR(1) >= 0.0d0 ) then ! H2O
      m = 1
      rm_VMR(:,m) = VMR(m)
      rm_VMR(:,m) = VMR(m) * ( 1.0d0 - rm_VMR(:,1) )
    else
      rm_VMR(:,1) = r_MMRH2O * MeanMolWt / m_MolWt(1)
    end if
  end if
  if ( NMol >= 2 ) then
    m = 2
    rm_VMR(:,m) = VMR(m)
!!$    rm_VMR(:,m) = VMR(m) * ( 1.0d0 - rm_VMR(:,1) )
  end if
  if ( NMol >= 3 ) then
    if ( VMR(3) >= 0.0d0 ) then ! O3
      m = 3
      rm_VMR(:,m) = VMR(m)
!!$      rm_VMR(:,m) = VMR(m) * ( 1.0d0 - rm_VMR(:,1) )
    else
      rm_VMR(:,3) = r_MMRO3  * MeanMolWt / m_MolWt(3)
    end if
  end if
  do m = 4, NMol
    rm_VMR(:,m) = VMR(m)
!!$    rm_VMR(:,m) = VMR(m) * ( 1.0d0 - rm_VMR(:,1) )
  end do


  ! Set particles
  !
  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
  select case ( IDPtcl )
  case ( IDPtclNone )
    NPtcl              = 0
    !
    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
  case ( IDPtclEarth )
    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

    ! 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
  case ( IDPtclMars )
    NPtcl              = 1
    !
    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

    Grav                = 3.72d0
    DustMMR0            = 1.0d0
    DustVerDistCoef     = 0.01d0
!!$    DustVerDistCoef     = 0.007d0
    DustMaxHeight       = 70.0d3
    DustVerDistRefPress = 7.0d2
    DustEffRadius       = 1.85d-6   ! Ockert-Bell et al. (1997)
    DustDens            = 2500.0_DP ! Pettengill and Ford (2000)
    DustExtEffFact067   = 3.04_DP   ! Ockert-Bell et al. (1997)

    ! dust, optical depth of 0.2
    a_PtclName(1)      = "mars_dust"
    a_PtclDens(1)      = DustDens
    za_PtclEffRad(:,1) = DustEffRadius

    m = 1
    do k = 1, NLev-1
      za_PtclMMR(k,m) = DustMMR0 &
        & * exp( DustVerDistCoef        &
        &        * ( 1.0_DP - ( DustVerDistRefPress / z_Press(k) )**(70.0d3/DustMaxHeight) ) &
        &      )
    end do
    DustColMass = 0.0d0
    do k = NLev-1, 1, -1
      DustColMass = DustColMass &
        & + za_PtclMMR(k,m) * ( r_Press(k) - r_Press(k+1) ) / Grav
    end do
    DustColNumDen = DustColMass / ( 4.0d0 / 3.0d0 * PI * DustEffRadius**3 * DustDens )
    DODFactor = DustColNumDen * PI * DustEffRadius**2 * DustExtEffFact067
    DODFactor = 0.2d0 / DODFactor * r_Press(1) / DustVerDistRefPress
    za_PtclMMR(:,m) = za_PtclMMR(:,m) * DODFactor
  case default
    write( 6, * ) 'Unexpected IDPtcl: ', IDPtcl
    stop
  end select


  call IOProfileOutput( &
    & OutNcFn, &
    & NLev, NMol, r_Press, m_MolNum, &
    & r_Temp, rm_VMR, Ps, Ts, &
    & z_Press, &
    & NChar, NPtcl, a_PtclNum, a_PtclName, &
    & a_PtclDens, za_PtclEffRad, za_PtclMMR &
    & )

end program main

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