module particle_ec_io_module

  use vtype_module
  use ni3_module
  use netcdf

  implicit none

  private

  public :: particle_ec_io_output_init
  public :: particle_ec_io_output_end
  public :: particle_ec_io_output
  public :: particle_ec_io_input_init
  public :: particle_ec_io_input
  public :: particle_ec_io_input_end


  integer          , save :: OutNcID
  integer          , save :: InNcID

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

contains

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

  subroutine particle_ec_io_output_init( &
    & FileName, NWN, &
    & NMolNum, m_MolNum, &
    & NParticle, a_ParticleName, a_ParticleNRadius, &
    & NRadiusMax, ra_ParticleRadius, &
    & MissingValue, DelWaveNum &
    & )


    character(*), intent(in ) :: FileName
    integer     , intent(in ) :: nwn
    integer     , intent(in ) :: NMolNum
    integer     , intent(in ) :: m_MolNum(NMolNum)
    integer     , intent(in ) :: NParticle
    character(*), intent(in ) :: a_ParticleName   (NParticle)
    integer     , intent(in ) :: a_ParticleNRadius(NParticle)
    integer     , intent(in ) :: NRadiusMax
    real(dp)    , intent(in ) :: ra_ParticleRadius(NRadiusMax,NParticle)
    real(DP)    , intent(in ) :: MissingValue
    real(DP)    , intent(in ) :: DelWaveNum


    !
    ! local variables
    !
    integer(i4b)         , parameter :: ndims1 = 1
    character(len=extstr)            :: dname1( ndims1 )
    integer(i4b)         , parameter :: ndims2 = 2
    character(len=extstr)            :: dname2( ndims2 )
    integer(i4b)         , parameter :: ndims3 = 3
    character(len=extstr)            :: dname3( ndims3 )
    character(len=extstr)            :: name
    character(len=extstr)            :: stdname
    character(len=extstr)            :: units
    integer(i4b)         , parameter :: NChar = stdstr
    integer(i4b)                     :: a_CharNum(NChar)
    integer(i4b)                     :: a_ParticleNum(NParticle)
    integer(i4b)                     :: r_ParticleRadiusNum(NRadiusMax)
    integer(i4b)                     :: l
    integer(i4b)                     :: m
    integer(i4b)                     :: n


    do l = 1, NChar
      a_CharNum(l) = l
    end do
    do n = 1, NParticle
      a_ParticleNum(n) = n
    end do
    do l = 1, NRadiusMax
      r_ParticleRadiusNum(l) = l
    end do


    call ni3_open( FileName, "new", OutNcID )

    call ni3_put_att( OutNcID, "_global_", "comment_dwn_m-1", DelWaveNum )

    name    = "WaveNum"
    stdname = "wavenumber"
    units   = "m-1"
    call ni3_def_dim( OutNcID, name, NI3_DOUBLE, nwn , &
!!$    call ni3_def_dim( OutNcID, name, NI3_DOUBLE, NI3_UNLIMITED , &
      & stdname = stdname, units  = units )

    name    = "MolNum"
    stdname = "molecular number"
    units   = "1"
    call ni3_set_dim( OutNcID, name, NI3_INT, m_MolNum, &
      & stdname = stdname, units  = units )

    if ( NParticle > 0 ) then
      name    = "CharNum"
      stdname = "character number"
      units   = "1"
      call ni3_set_dim( OutNcID, name, NI3_INT, a_CharNum, &
        & stdname = stdname, units  = units )
      !
      name    = "PtclNum"
      stdname = "particle number"
      units   = "1"
      call ni3_set_dim( OutNcID, name, NI3_INT, a_ParticleNum, &
        & stdname = stdname, units  = units )
      !
      name    = "PtclRadiusNum"
      stdname = "particle radius number"
      units   = "1"
      call ni3_set_dim( OutNcID, name, NI3_INT, r_ParticleRadiusNum, &
        & stdname = stdname, units  = units )


      dname1( 1 ) = "PtclNum"
      !
      name    = "PtclNRadius"
      stdname = "name of radius in each particle"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_INT, ndims1, dname1, &
        & stdname = stdname, units = units )
      call ni3_put_var( OutNcID, name, a_ParticleNRadius )
      !
      dname2( 1 ) = "CharNum"
      dname2( 2 ) = "PtclNum"
      !
      name    = "PtclName"
      stdname = "name of particles"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_CHAR, ndims2, dname2, &
        & stdname = stdname, units = units )
      call ni3_put_var( OutNcID, name, a_ParticleName )
      !
      dname2( 1 ) = "PtclRadiusNum"
      dname2( 2 ) = "PtclNum"
      !
      name    = "PtclRadius"
      stdname = "particle radius"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims2, dname2, &
        & stdname = stdname, units = units )
      call ni3_put_var( OutNcID, name, ra_ParticleRadius )
      call ni3_put_att( OutNcID, name, "missing_value", MissingValue )
    end if


    dname1( 1 ) = "WaveNum"
    !
    name    = "RayScatCoefNonRadAct"
    stdname = "Rayleigh_scattering_coefficient"
!!$    units   = "(kg m-2)-1"
      units   = "(molecule m-2)-1"
    call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims1, dname1, &
      & stdname = stdname, units = units )

    dname2( 1 ) = "WaveNum"
    dname2( 2 ) = "MolNum"
    !
    name    = "RayScatCoef"
    stdname = "Rayleigh_scattering_coefficient"
!!$    units   = "(kg m-2)-1"
      units   = "(molecule m-2)-1"
    call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims2, dname2, &
      & stdname = stdname, units = units )


    if ( NParticle > 0 ) then
      dname3( 1 ) = "WaveNum"
      dname3( 2 ) = "PtclRadiusNum"
      dname3( 3 ) = "PtclNum"
      !
      name    = "QExt"
      stdname = "Extinction efficiency factor of particles"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims3, dname3, &
        & stdname = stdname, units = units )
      call ni3_put_att( OutNcID, name, "missing_value", MissingValue )
      !
      name    = "SSA"
      stdname = "Single scattering albedo of particles"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims3, dname3, &
        & stdname = stdname, units = units )
      call ni3_put_att( OutNcID, name, "missing_value", MissingValue )
      !
      name    = "AF"
      stdname = "Asymmetry factor of particles"
      units   = "1"
      call ni3_def_var( OutNcID, name, NI3_DOUBLE, ndims3, dname3, &
        & stdname = stdname, units = units )
      call ni3_put_att( OutNcID, name, "missing_value", MissingValue )
    end if


  end subroutine particle_ec_io_output_init

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

  subroutine particle_ec_io_output_end


    call ni3_close( OutNcID )


  end subroutine particle_ec_io_output_end

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

  subroutine particle_ec_io_output( &
    & nwn, wn, w_RayScatCoefNonRadAct, NMolNum, wm_RayScatCoef, &
    & NRadiusMax, NParticle, a_ParticleName, wra_QExt, wra_SSA, wra_AF &
    & )


    integer(i4b), intent(in ) :: nwn
    real(dp)    , intent(in ) :: wn( nwn )
    real(dp)    , intent(in ) :: w_RayScatCoefNonRadAct( nwn )
    integer(i4b), intent(in ) :: NMolNum
    real(dp)    , intent(in ) :: wm_RayScatCoef( nwn, NMolNum )
    integer(i4b), intent(in ) :: NRadiusMax
    integer(i4b), intent(in ) :: NParticle
    character(*), intent(in ) :: a_ParticleName( NParticle )
    real(dp)    , intent(in ) :: wra_QExt( nwn, NRadiusMax, NParticle )
    real(dp)    , intent(in ) :: wra_SSA ( nwn, NRadiusMax, NParticle )
    real(dp)    , intent(in ) :: wra_AF  ( nwn, NRadiusMax, NParticle )


    !
    ! local variables
    !
    character(len=extstr)       :: err_mes = "in ac_io_output"


    call ni3_put_var( OutNcID, "WaveNum", wn )

    call ni3_put_var( OutNcID, 'RayScatCoefNonRadAct', w_RayScatCoefNonRadAct )
    call ni3_put_var( OutNcID, 'RayScatCoef'         , wm_RayScatCoef )

    if ( NParticle > 0 ) then
      call ni3_put_var( OutNcID, "QExt", wra_QExt )
      call ni3_put_var( OutNcID, "SSA" , wra_SSA  )
      call ni3_put_var( OutNcID, "AF"  , wra_AF   )
    end if

    call ni3_redef( OutNcID )
    call ni3_enddef( OutNcID )


  end subroutine particle_ec_io_output

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

  subroutine particle_ec_io_input_init( &
    & FileName, NWN &
    & )


    character(*), intent(in ) :: FileName
    integer     , intent(out) :: nwn


    !
    ! local variables
    !
    character(len=extstr)            :: name


    call ni3_open( FileName, "read", InNcID )

    name    = "WaveNum"
    call ni3_inq_dimlen( InNcID, name, NWN )


  end subroutine particle_ec_io_input_init

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

  subroutine particle_ec_io_input_end


    call ni3_close( InNcID )


  end subroutine particle_ec_io_input_end

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

  subroutine particle_ec_io_input( &
    & nwnsl, VarName, wn, sc, iwnsls, iwnsle, iwnl &
    & )


    integer(i4b), intent(in ) :: nwnsl
    character(*), intent(in ) :: VarName
    real(dp)    , intent(out) :: wn( nwnsl )
    real(dp)    , intent(out) :: sc( nwnsl )
    integer(i4b), intent(in ) :: iwnsls
    integer(i4b), intent(in ) :: iwnsle
    integer(i4b), intent(in ) :: iwnl


    !
    ! local variables
    !
    character(len=extstr)       :: name
    integer(i4b)                :: varid_wn
    integer(i4b)                :: varid_sc
    integer(i4b)                :: istatus
    character(len=extstr)       :: err_mes = "in ac_io_output"
    integer(i4b)                :: st1( 1 ), co1( 1 ), st2( 2 ), co2( 2 )


    name = "WaveNum"
    call ni3_inq_var( InNcID, name, varid_wn )
!!$    name = "AbsCoef"
    name = VarName
    call ni3_inq_var( InNcID, name, varid_sc )
    call ni3_enddef( InNcID )


    st1( 1 ) = iwnl + iwnsls - 1
    co1( 1 ) = iwnsle - iwnsls + 1
    istatus = nf90_get_var( InNcID, varid_wn, wn, start = st1, count = co1 )
    call ni3_handle_err( istatus, err_mes )


    st1( 1 ) = iwnl + iwnsls - 1
    co1( 1 ) = iwnsle - iwnsls + 1
    istatus = nf90_get_var( InNcID, varid_sc, sc, start = st1, count = co1 )
    call ni3_handle_err( istatus, err_mes )


  end subroutine particle_ec_io_input

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

end module particle_ec_io_module
