program main

  use vtype_module

#ifdef MPI
  use mpi
#endif

  use fi_module

  use io_profile, only : &
    & IOProfileInqDimLen, &
    & IOProfileRead

  use ac_calc, only : ACCalc

  use ac_calc_utils, only : &
    & ACCalcUtilsModifyWavenumberRange1, &
    & ACCalcUtilsModifyWavenumberRange2


  implicit none


  character(extstr) :: CtlFN  = 'calc_ac.conf'
  character(extstr) :: CtlFN2 = 'calc_ac_file.conf'

  integer           :: iostat

  logical , allocatable :: m_FlagCalcLine(:)
  logical , allocatable :: m_FlagCalcCont(:)

  integer, parameter :: MaxNmlArySize = 64


  logical                            :: FlagNetCDF

  character(len=extstr)              :: HITFileName(MaxNmlArySize)
  integer                            :: NHIT
  character(len=extstr), allocatable :: a_HITFileName(:)
  character(len=extstr)              :: OutputFileName

  integer                       :: NLev

  !
  ! NMolNum : number of radiatively active species
  !
  integer                        :: NMolNum


  real(dp)                           :: DelWaveNum
  real(dp)                           :: WaveNumS
  real(dp)                           :: WaveNumE


  integer(i4b)              :: FileUnit



  character(len=extstr)     :: InFileName
  character(len=extstr)     :: PressName
  character(len=extstr)     :: MolNumName
  character(len=extstr)     :: TempName
  character(len=extstr)     :: SurfTempName
  character(len=extstr)     :: VMRName

  integer     , allocatable :: m_MolNum( : )
  real(dp)    , allocatable :: r_Press(:)
  real(dp)    , allocatable :: r_Temp(:)
  real(dp)                  :: SurfTemp
  real(dp)    , allocatable :: rm_VMR(:,:)


  real(dp)    , allocatable :: z_Press(:)
  real(dp)    , allocatable :: z_Temp(:)
  real(dp)    , allocatable :: zm_VMR(:,:)
  real(dp)    , allocatable :: zm_VMRShrink(:,:)


  logical                            :: FlagCalcLine(MaxNmlArySize)
  logical                            :: FlagCalcCont(MaxNmlArySize)

  real(dp)                           :: LineCutOffWaveNum(MaxNmlArySize)
  real(dp), allocatable              :: m_LineCutOffWaveNum(:)

  real(dp)                           :: AirBroadHalfWidFac(MaxNmlArySize)
  real(dp), allocatable              :: m_AirBroadHalfWidFac(:)

  character(len=extstr)              :: CIAFileName(MaxNmlArySize)
  integer                            :: NCIA
  character(len=extstr), allocatable :: a_CIAFileName(:)

  character(len=extstr)              :: XSCFileName(MaxNmlArySize)
  integer                            :: NXSC
  character(len=extstr), allocatable :: a_XSCFileName(:)



  integer :: ks
  integer :: ke

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


  integer :: NProcs
  integer :: MyRank
#ifdef MPI
  integer :: iRank
  integer :: IErr
#endif

  integer, parameter :: FileNumS = 0


  namelist /wn_nml/ &
    & DelWaveNum, WaveNumS, WaveNumE
  namelist /gas_nml/ &
    & FlagCalcLine, FlagCalcCont, &
    & FlagNetCDF, &
    & HITFileName, CIAFileName, XSCFileName, &
    & LineCutOffWaveNum, &
    & AirBroadHalfWidFac
  namelist /file_nml/ &
    & InFileName, &
    & OutputFileName, &
    & PressName, MolNumName, TempName, SurfTempName, VMRName


#ifdef MPI
  call mpi_init( IErr )
  call mpi_comm_size( mpi_comm_world, NProcs, IErr )
  call mpi_comm_rank( mpi_comm_world, MyRank, IErr )
#else
  NProcs = 1
  MyRank = 0
#endif



  FlagCalcLine(:) = .true.
  FlagCalcCont(:) = .true.

  FlagNetCDF = .true.
  HITFileName(:) = ""
  !
  CIAFileName(:)  = ""
!  CIAFileName(1)  = "../data/CIA/NC_data/02-02_cia.nc"
  !
  XSCFileName(:)  = ""
  !
  OutputFileName = "./out/ac.nc"

!!$      ! for 15 mincron band
!!$      WaveNumS =  400.0d0 * 1.0d2
!!$      WaveNumE =  900.0d0 * 1.0d2
!!$      ! for CO2 500-900 cm^-1 for validation
!!$      WaveNumS =  500.0d0 * 1.0d2
!!$      WaveNumE =  900.0d0 * 1.0d2
  ! for 15 mincron band for validation
  DelWaveNum = 1.0d-2 * 1.0d2
  WaveNumS =  550.0d0 * 1.0d2 + DelWaveNum / 2.0d0
  WaveNumE =  750.0d0 * 1.0d2 - DelWaveNum / 2.0d0

  WaveNumS =  600.0d0 * 1.0d2 + DelWaveNum / 2.0d0
  WaveNumE =  601.0d0 * 1.0d2 - DelWaveNum / 2.0d0


!!$      WaveNumS =    0.0d0 * 1.0d2
!!$      WaveNumE = 3000.0d0 * 1.0d2

  LineCutOffWaveNum(:) = 25.0d2

  AirBroadHalfWidFac = 1.0d0

  !======================
  InFileName   = "./data.nc.nc"
  PressName    = "Press"
  MolNumName   = "MolNum"
  TempName     = "Temp"
  SurfTempName = "Ts"
  VMRName      = "VMR"
  !======================



  ! Read a namelist file

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

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



  call IOProfileInqDimLen(      &
    & InFileName, PressName, &
    & NLev &
    & )
  call IOProfileInqDimLen(      &
    & InFileName, MolNumName, &
    & NMolNum &
    & )
  allocate( m_MolNum( NMolNum ) )
  allocate( r_Press ( NLev    ) )
  allocate( r_Temp  ( NLev    ) )
  allocate( rm_VMR  ( NLev, NMolNum ) )
  call IOProfileRead(      &
    & InFileName, PressName, MolNumName, TempName, SurfTempName, VMRName, &
    & NLev, NMolNum, &
    & r_Press, m_MolNum, r_Temp, SurfTemp, rm_VMR     &
    & )

  allocate( z_Press ( NLev-1    ) )
  allocate( z_Temp  ( NLev-1    ) )
  allocate( zm_VMR  ( NLev-1, NMolNum ) )
  do k = 1, NLev-1
    z_Press ( k ) = ( r_Press(k) + r_Press(k+1) ) / 2.0d0
    z_Temp  ( k ) = ( r_Temp (k) + r_Temp (k+1) ) / 2.0d0
  end do
  do m = 1, NMolNum
    do k = 1, NLev-1
      zm_VMR  ( k, m ) = ( rm_VMR(k,m) + rm_VMR(k+1,m) ) / 2.0d0
    end do
  end do


  allocate( m_FlagCalcLine(NMolNum) )
  allocate( m_FlagCalcCont(NMolNum) )
  m_FlagCalcLine = FlagCalcLine(1:NMolNum)
  m_FlagCalcCont = FlagCalcCont(1:NMolNum)
  write( 6, * ) 'Flags for calculation of lines and continuum:'
  do m = 1, NMolNum
    write( 6, '(a,l,a,l,a,i2)' ) &
      & '  Line = ', m_FlagCalcLine(m), &
      & ', Continuum = ', m_FlagCalcCont(m), &
      & ' for molecule ', m_MolNum(m)
  end do

  allocate( m_LineCutOffWaveNum(NMolNum) )
  m_LineCutOffWaveNum = LineCutOffWaveNum(1:NMolNum)
  write( 6, * ) 'Line cut off wavenumbers:'
  do m = 1, NMolNum
    write( 6, '(a,i4,a,i2)' ) &
      & '  Cut-off at ', int( m_LineCutOffWaveNum(m)*1d-2 ), &
      & ' cm-1 for molecule ', m_MolNum(m)
  end do

  allocate( m_AirBroadHalfWidFac(NMolNum) )
  m_AirBroadHalfWidFac = AirBroadHalfWidFac(1:NMolNum)
  write( 6, * ) 'Air-broadening half width factors:'
  do m = 1, NMolNum
    write( 6, '(a,f4.2,a,i2)' ) &
      & '  Air-broadening half width is multiplied by ', m_AirBroadHalfWidFac(m), &
      & ' for molecule ', m_MolNum(m)
  end do


  NHIT = 0
  do n = 1, MaxNmlArySize
    if ( len_trim( HITFileName(n) ) == 0 ) exit
  end do
  NHIT = n - 1
  allocate( a_HITFileName(NHIT) )
  do n = 1, NHIT
    a_HITFileName(n) = HITFileName(n)
  end do
  ! check
  if ( NHIT /= NMolNum ) then
    write( 6, * ) 'Number of HITRAN data files specified in gas_nml does not match to number of species in file, ', &
      trim( InFileName ), '.'
    write( 6, * ) 'Number of HITRAN data files: ', NHIT
    write( 6, * ) 'Number of species          : ', NMolNum
    stop
  end if

  NCIA = 0
  do n = 1, MaxNmlArySize
    if ( len_trim( CIAFileName(n) ) == 0 ) exit
  end do
  NCIA = n - 1
  allocate( a_CIAFileName(NCIA) )
  do n = 1, NCIA
    a_CIAFileName(n) = CIAFileName(n)
  end do

  NXSC = 0
  do n = 1, MaxNmlArySize
    if ( len_trim( XSCFileName(n) ) == 0 ) exit
  end do
  NXSC = n - 1
  allocate( a_XSCFileName(NXSC) )
  do n = 1, NXSC
    a_XSCFileName(n) = XSCFileName(n)
  end do


#ifdef MPI

!!$  if ( .true. ) then
  if ( .false. ) then
    write( 6, * ) 'Merge wavenumbers'
    !
    ! MPI for wavenumbers
    !
!  call ACCalcUtilsModifyWavenumberRange1( &
!    & NProcs, MyRank, &
!    & DelWaveNum, &
!    & WaveNumS, WaveNumE &
!    & )
    call ACCalcUtilsModifyWavenumberRange2( &
      & FlagNetCDF, NHIT, a_HITFileName, &
      & NProcs, MyRank, &
      & NMolNum, m_MolNum, &
      & DelWaveNum, &
      & WaveNumS, WaveNumE &
      & )

    write( OutputFileName, '(a,a,i6.6,a)' ) &
      & OutputFileName(1:len_trim(OutputFileName)-3), &
      & "_", &
      & FileNumS + MyRank, &
      & ".nc"
    if ( MyRank == 0 ) then
      write( 6, * ) "# Wavenumber range is divided by MPI."
    end if
    write( 6, * ) "#", MyRank, ":", WaveNumS, WaveNumE, trim( OutputFileName )

    ks = 1
    ke = NLev-1

  else
    write( 6, * ) 'Merge vertical domain'
    !
    ! MPI for vertical levels
    !
    if ( (NLev-1) < NProcs ) then
      write( 6, * ) 'Number of processes is too large: ', NLev-1, NProcs
      stop
    end if
    !
    ke = 0
    do iRank = 0, MyRank
      ks = ke + 1
      ke = ks + (NLev-1) / NProcs - 1
      if ( iRank < mod( (NLev-1), NProcs ) ) then
        ke = ke + 1
      end if
    end do
    ! check
    if ( ke > (NLev-1) ) then
      write( 6, * ) 'Unexpected number of maximum of loop: ', ke, NLev-1
    end if

    write( OutputFileName, '(a,a,i6.6,a)' ) &
      & OutputFileName(1:len_trim(OutputFileName)-3), &
      & "_", &
      & FileNumS + MyRank, &
      & ".nc"
    if ( MyRank == 0 ) then
      write( 6, * ) "# Vertical domain is divided by MPI."
    end if
    write( 6, * ) "#", MyRank, ":", ks, ke, trim( OutputFileName )

  end if

#else
  ks = 1
  ke = NLev-1
#endif

  allocate( zm_VMRShrink( ks:ke, NMolNum ) )

  do m = 1, NMolNum
    do k = ks, ke
      zm_VMRShrink(k,m) = zm_VMR(k,m)
    end do
  end do


!!$  call LBLInit

  call ACCalc( &
    & DelWaveNum, WaveNumS, WaveNumE, &
    & FlagNetCDF, NHIT, a_HITFileName, &
    & NCIA, a_CIAFileName, &
    & NXSC, a_XSCFileName, &
    & OutputFileName, &
!!$    & NLev, NMolNum, r_Press, m_MolNum, r_Temp, rm_VMR, &
!!$    & NLev-1, NMolNum, z_Press, m_MolNum, z_Temp, zm_VMR, &
    & ke-ks+1, NMolNum, z_Press(ks:ke), m_MolNum, z_Temp(ks:ke), zm_VMRShrink(ks:ke,:), &
    & m_FlagCalcLine, m_FlagCalcCont, &
    & m_LineCutOffWaveNum, m_AirBroadHalfWidFac &
    & )

  deallocate( z_Press  )
  deallocate( z_Temp   )
  deallocate( zm_VMR   )

  deallocate( zm_VMRShrink )

  deallocate( a_HITFileName )
  deallocate( a_CIAFileName )
  deallocate( a_XSCFileName )
  deallocate( m_MolNum )
  deallocate( r_Press  )
  deallocate( r_Temp   )
  deallocate( rm_VMR   )
  deallocate( m_LineCutOffWaveNum )
  deallocate( m_FlagCalcLine )
  deallocate( m_FlagCalcCont )

  write( 6, * ) 'Finish'


#ifdef MPI
  call mpi_finalize( IErr )
#endif


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

!!$contains
!!$
!!$  !----------------------------------------------------------------------------
!!$
!!$  subroutine ModifyWavenumberRange1( &
!!$    & NProcs, MyRank, &
!!$    & DelWaveNum, &
!!$    & WaveNumS, WaveNumE &
!!$    & )
!!$
!!$    integer , intent(in   ) :: NProcs
!!$    integer , intent(in   ) :: MyRank
!!$    real(dp), intent(in   ) :: DelWaveNum
!!$    real(dp), intent(inout) :: WaveNumS
!!$    real(dp), intent(inout) :: WaveNumE
!!$
!!$
!!$    ! Local variables
!!$    !
!!$    integer :: NWaveNum
!!$    integer :: iRank
!!$
!!$
!!$    NWaveNum = int( ( WaveNumE - WaveNumS ) / DelWaveNum ) + 1
!!$    WaveNumE = WaveNumS - DelWaveNum
!!$    do iRank = 0, MyRank
!!$      WaveNumS = WaveNumE + DelWaveNum
!!$      WaveNumE = WaveNumS + ( int( NWaveNum / NProcs ) - 1 ) * DelWaveNum
!!$      if ( iRank+1 <= mod( NWaveNum, NProcs ) ) then
!!$        WaveNumE = WaveNumE + DelWaveNum
!!$      end if
!!$    end do
!!$
!!$
!!$  end subroutine ModifyWavenumberRange1
!!$
!!$  !----------------------------------------------------------------------------
!!$
!!$  subroutine ModifyWavenumberRange2( &
!!$    & FlagNetCDF, HITRANDir, HITRANFileNameExt, &
!!$    & NProcs, MyRank, &
!!$    & NMolNum, m_MolNum, &
!!$    & DelWaveNum, &
!!$    & WaveNumS, WaveNumE &
!!$    & )
!!$
!!$    use lbl_module
!!$
!!$    logical         , intent(in   ) :: FlagNetCDF
!!$    character(len=*), intent(in   ) :: HITRANDir
!!$    character(len=*), intent(in   ) :: HITRANFileNameExt
!!$    integer         , intent(in   ) :: NProcs
!!$    integer         , intent(in   ) :: MyRank
!!$    integer         , intent(in   ) :: NMolNum
!!$    integer         , intent(in   ) :: m_MolNum(NMolNum)
!!$    real(dp)        , intent(in   ) :: DelWaveNum
!!$    real(dp)        , intent(inout) :: WaveNumS
!!$    real(dp)        , intent(inout) :: WaveNumE
!!$
!!$
!!$    ! Local variables
!!$    !
!!$    integer, parameter :: MolNumChk = 2
!!$    logical            :: FlagClear
!!$
!!$    type( LineParam )         :: LParam
!!$    character(len=extstr)     :: FileName
!!$
!!$    integer  :: iLineS
!!$    integer  :: iLineE
!!$    integer  :: NLine
!!$    real(dp) :: a_WaveNumBnds(2,0:NProcs-1)
!!$
!!$    integer :: l
!!$    integer :: m
!!$
!!$
!!$    FlagClear = .false.
!!$    do m = 1, NMolNum
!!$      if ( m_MolNum(m) == MolNumChk ) then
!!$        FlagClear = .true.
!!$      end if
!!$    end do
!!$    if ( .not. FlagClear ) then
!!$      write( 6, * ) 'Unable to find ', MolNumChk, 'th molecule.'
!!$      write( 6, * ) ' ', m_MolNum
!!$      stop
!!$    end if
!!$
!!$
!!$    m = MolNumChk
!!$
!!$    write( FileName, '(a,i2.2,a)' ) &
!!$      & trim( HITRANDir ) // "/", m_MolNum( m ), trim( HITRANFileNameExt )
!!$    write( 6, '(a)' ) '  Line File: ' // trim( FileName )
!!$
!!$    if ( FlagNetCDF ) then
!!$      call lbl_readfile_netcdf( FileName, 1, LParam )
!!$      write( 6, '(a,i9)' ) '             Number of lines : ', LParam%NLine
!!$    else
!!$      write( 6, * ) 'Ascii format HITRAN data file is not supported, now.'
!!$      stop
!!$    end if
!!$
!!$
!!$    do l = 1, LParam%NLine
!!$      if ( WaveNumS <= LParam%wnlc0(l) ) exit
!!$    end do
!!$    iLineS = l
!!$    iLineS = min( max( iLineS, 1 ), LParam%NLine )
!!$    do l = iLineS, LParam%NLine
!!$      if ( WaveNumE <= LParam%wnlc0(l) ) exit
!!$    end do
!!$    iLineE = l
!!$    iLineE = min( max( iLineE, iLineS ), LParam%NLine )
!!$
!!$    NLine = iLineE - iLineS + 1
!!$
!!$    if ( NLine < NProcs ) then
!!$      write( 6, * ) 'Number of processes are too large.'
!!$      write( 6, * ) '  ', NLine, NProcs
!!$      stop
!!$    end if
!!$
!!$    do n = 0, NProcs-1
!!$      if ( n == 0 ) then
!!$        a_WaveNumBnds(1,n) = WaveNumS
!!$      else
!!$        a_WaveNumBnds(1,n) = a_WaveNumBnds(2,n-1) + DelWaveNum
!!$      end if
!!$      if ( n == (NProcs-1) ) then
!!$        a_WaveNumBnds(2,n) = WaveNumE
!!$      else
!!$        l = iLineS + int( NLine / NProcs ) * ( n+1 )
!!$        a_WaveNumBnds(2,n) = WaveNumS &
!!$          & + int( max( LParam%wnlc0(l) - WaveNumS, 0.0d0 ) / DelWaveNum ) * DelWaveNum
!!$      end if
!!$    end do
!!$
!!$    WaveNumS = a_WaveNumBnds(1,MyRank)
!!$    WaveNumE = a_WaveNumBnds(2,MyRank)
!!$
!!$    if ( WaveNumS == WaveNumE ) then
!!$      write( 6, * ) 'Unexpected zero wavenumber range for process ', MyRank
!!$      stop
!!$    end if
!!$
!!$    call lbl_lineparam_finalize( LParam )
!!$
!!$
!!$  end subroutine ModifyWavenumberRange2

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

end program main
