program main

  use vtype_module
  use netcdf
  use ni3_module

  use fi_module

  use set_1d_profile, only : Set1DProfileInterpolate
  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                            :: NLevIn
  real(8)              , allocatable :: r_PressIn  (:)
  real(8)              , allocatable :: r_TempIn   (:)
  real(8)              , allocatable :: rm_VMRIn   (:,:)
  real(8)              , allocatable :: rc_CldMMRIn(:,:)


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

  integer :: FileUnit
  integer :: iostat

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

  data m_HITRANMolName &
    / &
    & "H2O     ", "CO2     ", "O3      ",  &
    & "N2O     ", "CO      ", "CH4     ",  &
    & "O2      ", "NO      ", "SO2     ",  &
    & "NO2     ", "NH3     ", "HNO3    ",  &
    & "OH      ", "HF      ", "HCl     ",  &
    & "HBr     ", "HI      ", "ClO     ",  &
    & "OCS     ", "H2CO    ", "HOCl    ",  &
    & "N2      ", "HCN     ", "CH3Cl   ",  &
    & "H2O2    ", "C2H2    ", "C2H6    ",  &
    & "PH3     ", "COF2    ", "SF6     ",  &
    & "H2S     ", "HCOOH   ", "HO2     ",  &
    & "O       ", "ClONO2  ", "NO+     ",  &
    & "HOBr    ", "C2H4    ", "CH3OH   ",  &
    & "CH3Br   ", "CH3CN   ", "CF4     ",  &
    & "C4H2    ", "HC3N    ", "H2      ",  &
    & "CS      ", "SO3     ", "C2N2    ",  &
    & "COCl2   "   &
    /

  data m_HITRANMolWt &
    / &
    &    1.80152636E-02 ,    4.40097749E-02 ,    4.79970984E-02 ,  &
    &    4.40113507E-02 ,    2.80104484E-02 ,    1.60430688E-02 ,  &
    &    3.19984369E-02 ,    2.99942717E-02 ,    6.32543117E-02 ,  &
    &    4.56073023E-02 ,    1.70222521E-02 ,    6.25422969E-02 ,  &
    &    1.70005895E-02 ,    2.00063828E-02 ,    3.64607014E-02 ,  &
    &    8.09115693E-02 ,   0.127912432     ,    5.13256155E-02 ,  &
    &    5.99990338E-02 ,    3.00044846E-02 ,    5.23229167E-02 ,  &
    &    2.80182287E-02 ,    2.70203035E-02 ,    4.98920977E-02 ,  &
    &    3.38338241E-02 ,    2.60346141E-02 ,    3.00372038E-02 ,  &
    &    3.39813642E-02 ,    6.58462718E-02 ,   0.138690650     ,  &
    &    3.40629518E-02 ,    4.52647023E-02 ,    3.28361988E-02 ,  &
    &    1.59569755E-02 ,    9.63944271E-02 ,    2.98172217E-02 ,  &
    &    9.66614708E-02 ,    2.80324090E-02 ,    3.15756090E-02 ,  &
    &    9.38285887E-02 ,    3.99543643E-02 ,    8.70160088E-02 ,  &
    &    4.78148647E-02 ,    4.91411462E-02 ,    2.01596250E-03 ,  &
    &    4.40405197E-02 ,    7.54312649E-02 ,    5.04850745E-02 ,  &
    &    9.16662812E-02    &
    /

!!$  data a_CldH13Zb &
!!$    & / 49.0d3, 65.0d3, 49.0d3, 49.0d3 /
!!$  data a_CldH13Zc &
!!$    & / 16.0d3,  1.0d3, 11.0d3,  8.0d3 /
!!$  data a_CldH13Hup &
!!$    & /  5.0d3,  3.5d3,  1.0d3,  1.0d3 /
!!$  data a_CldH13Hlo &
!!$    & /  1.0d3,  3.0d3,  0.1d3,  0.5d3 /
!!$  data a_CldH13N0 &
!!$    & / 181.0d-6, 100.0d-6,  50.0d-6,  14.0d-6 /

  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, "plev", 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   )
    !
    allocate( rm_VMR   ( NLev, NMol) )
    do m = 1, NMol
      Name = 'VMR' // trim( m_HITRANMolName( m_MolNum(m) ) )
      call ni3_get_var( NcID, Name, rm_VMR(:,m) )
    end do

    call ni3_inq_dimlen( NcID, "cldnum", NCldNum )
    allocate( rc_CldMMR( NLev, NCldNum) )
    Name = 'CldMMR'
    call ni3_get_var( NcID, Name, rc_CldMMR )

    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

    Mode = "read"
    call ni3_open( InNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, "plev"  , NLevIn )

    allocate( r_PressIn  ( NLevIn ) )
    allocate( r_TempIn   ( NLevIn ) )
    allocate( rm_VMRIn   ( NLevIn, NMol ) )
    !
    call ni3_get_var( NcID, PressName , r_PressIn  )
    call ni3_get_var( NcID, TempName  , r_TempIn   )
    !
    do m = 1, NMol
      Name = 'VMR' // trim( m_HITRANMolName( m_MolNum(m) ) )
      call ni3_get_var( NcID, Name, rm_VMRIn(:,m) )
    end do

    call ni3_inq_dimlen( NcID, "cldnum", NCldNum )
    allocate( rc_CldMMRIn( NLevIn, NCldNum ) )
    Name = 'CldMMR'
    call ni3_get_var( NcID, Name, rc_CldMMRIn )

    call ni3_close( NcID )



    allocate( r_Sigma ( NLev ) )

    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

    allocate( r_Press  ( NLev ) )
    allocate( r_Temp   ( NLev ) )
    allocate( rm_VMR   ( NLev, NMol) )
    allocate( rc_CldMMR( NLev, NCldNum ) )


    Ps = r_PressIn(1)

    r_Press = Ps * r_Sigma

    if ( TempConst > 0.0d0 ) then
      r_Temp = TempConst
    else
      call Set1DProfileInterpolate(     &
        & NLevIn, r_PressIn, r_TempIn,  &
        & NLev  , r_Press  ,            &
        & r_Temp                        &
        & )
    end if
    Ts = r_Temp (1)

    do m = 1, NMol
      call Set1DProfileInterpolate(          &
        & NLevIn, r_PressIn, rm_VMRIn(:,m),  &
        & NLev  , r_Press  ,                 &
        & rm_VMR(:,m)                        &
        & )
    end do


    do m = 1, NCldNum
      call Set1DProfileInterpolate(             &
        & NLevIn, r_PressIn, rc_CldMMRIn(:,m),  &
        & NLev  , r_Press  ,                    &
        & rc_CldMMR(:,m)                        &
        & )
    end do


    deallocate( r_PressIn   )
    deallocate( r_TempIn    )
    deallocate( rm_VMRIn    )
    deallocate( rc_CldMMRIn )

    deallocate( r_Sigma )

  end if


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

  ! Venus cloud
  !
  ! Pollack et al. (1993)
!!$  if ( NPtcl /= 4 ) then
!!$    write( 6, * ) 'Unexpected NPtcl: ', NPtcl
!!$    stop
!!$  end if
!!$  a_PtclName(1)      = "vcldmode1"
!!$  a_PtclName(2)      = "vcldmode2"
!!$  a_PtclName(3)      = "vcldmode2'"
!!$  a_PtclName(4)      = "vcldmode3"
!!$  a_PtclDens         = 1.6692d3
!!$  za_PtclEffRad(:,1) = 0.49d-6
!!$  za_PtclEffRad(:,2) = 1.18d-6
!!$  za_PtclEffRad(:,3) = 1.40d-6
!!$  za_PtclEffRad(:,4) = 3.65d-6
  !
  ! Crisp (1986)
  if ( NPtcl /= 5 ) then
    write( 6, * ) 'Unexpected NPtcl: ', NPtcl
    stop
  end if
  a_PtclName(1)      = "venus_cloud_mode1"
  a_PtclName(2)      = "venus_cloud_mode2"
  a_PtclName(3)      = "venus_cloud_mode2'"
  a_PtclName(4)      = "venus_cloud_mode3"
  a_PtclName(5)      = "venus_uv_absorber"
  a_PtclDens         = 1.6692d3
  za_PtclEffRad(:,1) = 0.49d-6
  za_PtclEffRad(:,2) = 1.18d-6
  za_PtclEffRad(:,3) = 1.40d-6
  za_PtclEffRad(:,4) = 3.65d-6
  za_PtclEffRad(:,5) = 0.49d-6
  !
  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


!!$  do m = 1, NPtcl
!!$    if ( m == 1 ) cycle ! mode 1
!!$!     if ( m == 2 ) cycle ! mode 2
!!$!     if ( m == 3 ) cycle ! mode 2'
!!$!     if ( m == 4 ) cycle ! mode 3
!!$    if ( m == 5 ) cycle ! UV absorber
!!$    za_PtclMMR(:,m) = za_PtclMMR(:,m) * 0.5d0
!!$  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
