module ac_calc

  use vtype_module

  implicit none

  private

  public :: ACCalc


  !
  ! NWaveNumSubLoop : number of wavenumber sub loop 
  !                 : (inner most loop for optimization)
  !
!!$  integer, parameter :: NWaveNumSubLoop =   10
  integer, parameter :: NWaveNumSubLoop = 4096 ! for pi-computer


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

contains

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

  subroutine ACCalc( &
    & DelWaveNum, WaveNumS, WaveNumE, &
    & FlagNetCDF, NHIT, a_HITFileName, &
    & NCIA, a_CIAFileName, &
    & NXSC, a_XSCFileName, &
    & OutputFileName, &
    & NLev, NMolNum, z_Press, m_MolNum, z_Temp, zm_VMR, &
    & m_FlagCalcLine, m_FlagCalcCont, &
    & m_LineCutOffWaveNum, m_AirBroadHalfWidFac &
    & )

    use constants0

    use fi_module
    use ac_io_module

    use lbl_module
    use mtckd_wrapper3
    use cia_module
    use xsc_module


    real(dp)    , intent(in ) :: DelWaveNum
    real(dp)    , intent(in ) :: WaveNumS
    real(dp)    , intent(in ) :: WaveNumE
    logical     , intent(in ) :: FlagNetCDF
    integer     , intent(in ) :: NHIT
    character(*), intent(in ) :: a_HITFileName( NHIT )
    integer     , intent(in ) :: NCIA
    character(*), intent(in ) :: a_CIAFileName( NCIA )
    integer     , intent(in ) :: NXSC
    character(*), intent(in ) :: a_XSCFileName( NXSC )
    character(*), intent(in ) :: OutputFileName
    integer     , intent(in ) :: NLev
    integer     , intent(in ) :: NMolNum
    real(dp)    , intent(in ) :: z_Press( NLev )
    integer     , intent(in ) :: m_MolNum( NMolNum )
    real(dp)    , intent(in ) :: z_Temp( NLev )
    real(dp)    , intent(in ) :: zm_VMR( NLev, NMolNum )
    logical     , intent(in ) :: m_FlagCalcLine( NMolNum )
    logical     , intent(in ) :: m_FlagCalcCont( NMolNum )
    real(dp)    , intent(in ) :: m_LineCutOffWaveNum( NMolNum )
    real(dp)    , intent(in ) :: m_AirBroadHalfWidFac( NMolNum )

!!$    allocate( z_Press ( NLev    ) )
!!$    allocate( m_MolNum( NMolNum ) )
!!$    allocate( z_Temp  ( NLev    ) )
!!$    allocate( zm_VMR  ( NLev, NMolNum ) )


    ! Local variables
    !
    logical                        :: FlagCIA = .true.
    logical                        :: FlagXSC = .true.

    type( LineParam )              :: LParam
    integer                        :: NLine

    type( CIAParam )               :: a_CParam( NCIA )
    integer                        :: IndexMolOpp
    type( XSCParam )               :: a_XSCP  ( NXSC )


    integer                   :: iwnsl
    integer                   :: iwnsls
    integer                   :: iwnsle


    real(dp)             , allocatable :: w_WaveNum( : )
    integer(i4b)                       :: NWaveNum


    character(len=extstr)     :: FileName

    real(dp)                  :: DelWaveNumMod




    !
    ! ach( NWaveNumSubLoop ) :: absorption coefficient
    !
    real(dp)    , allocatable :: w_AbsCoef(:)

    real(dp)    , allocatable :: w_AbsCoefLine   (:)
    real(dp)    , allocatable :: w_AbsCoefCont   (:)
    real(dp)    , allocatable :: w_AbsCoefCIA    (:)
    real(dp)    , allocatable :: w_AbsCoefCIAEach(:)
    real(dp)    , allocatable :: w_AbsCoefXSC    (:)
    real(dp)    , allocatable :: w_AbsCoefXSCEach(:)

    real(dp) :: VMRH2O
    real(dp) :: VMRCO2
    real(dp) :: VMRO3
    real(dp) :: VMRO2
    real(dp) :: VMRN2
    real(dp) :: VMROTHERS

    real(dp) :: NumDen
    real(dp) :: PathLength
    real(dp) :: ColDen

    logical :: FlagOutputWaveNum

    integer            :: iWaveNum

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



    allocate( w_WaveNum ( NWaveNumSubLoop ) )

    allocate( w_AbsCoef( NWaveNumSubLoop ) )


    NWaveNum   = int( ( WaveNumE - WaveNumS ) / DelWaveNum ) + 1
!!$    DelWaveNumMod = ( WaveNumE - WaveNumS ) / dble( NWaveNum )
    DelWaveNumMod = DelWaveNum


    allocate( w_AbsCoefLine   ( NWaveNumSubLoop ) )
    allocate( w_AbsCoefCont   ( NWaveNumSubLoop ) )
    allocate( w_AbsCoefCIA    ( NWaveNumSubLoop ) )
    allocate( w_AbsCoefCIAEach( NWaveNumSubLoop ) )
    allocate( w_AbsCoefXSC    ( NWaveNumSubLoop ) )
    allocate( w_AbsCoefXSCEach( NWaveNumSubLoop ) )

    call ac_io_output_init( &
      & OutputFileName, NMolNum, m_MolNum, NLev, z_Press, zm_VMR, z_Temp, DelWaveNum &
      & )


    if ( FlagCIA ) then
      do n = 1, NCIA
        call CIAReadFileNetcdf( a_CIAFileName(n), a_CParam(n) )
      end do
    end if

    if ( FlagXSC ) then
      do n = 1, NXSC
        call XSCReadFileNetcdf( a_XSCFileName(n), a_XSCP(n) )
      end do
    end if


    write( 6, * )

    ! check
    if ( NHIT /= NMolNum ) then
      write( 6, * ) 'Number of HITRAN data files does not match to number of species.'
      stop
    end if

    do m = 1, NMolNum
      write( 6, '(a)' ) '================================================='
      write( 6, '(a,i2,a,i2)' ) 'm = ', m, ', Molecular number: ', m_MolNum(m)

      if ( m_FlagCalcLine(m) ) then

        FileName = a_HITFileName(m)
        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
          NLine = fi_lc( FileName )
          write( 6, '(a,i9)' ) '             Number of lines : ', NLine

!!$          call lbl_lineparam_init( m_MolNum( m ), NLine, NLev, LParam )
          call lbl_lineparam_init( m_MolNum( m ), NLine, 1, LParam )

!!$          call lbl_readfile         ( FileName, LParam )
          call lbl_readfile_hitran04( FileName, LParam )

        end if

        ! check
        !   This is not necessary because this is checked in lbl_readfile_*.
        if ( m_MolNum(m) /= LParam%IMol ) then
          write( 6, * ) 'File, ', trim( FileName ), ' does not include data for moleculer #', m_MolNum(m), '.'
          stop
        end if

        !
        ! wavenumber width over integration (m^-1)
        !
        do k = 1, LParam%km
!!$    LParam % wnco( : ) = 3.0d2
!!$    LParam % wnco( : ) = 25.0d2
          LParam % wnco1    ( :, k ) = m_LineCutOffWaveNum (m)
          LParam % abwidfac1( :, k ) = m_AirBroadHalfWidFac(m)
        end do

      end if

      if ( FlagCIA ) then
        do n = 1, NCIA
          if ( a_CParam(n)%IMol == m_MolNum(m) ) then
            write( 6, '(a)' ) &
              & '  CIA File : ' // trim( a_CIAFileName(n) )
            write( 6, '(a,i8,a,i8,a)' ) &
              & '             ', &
              & int( a_CParam(n)%w_WaveNum(1) ), '-', &
              & int( a_CParam(n)%w_WaveNum(a_CParam(n)%NWaveNum) ), ' m-1'
          end if
        end do
      end if

      if ( FlagXSC ) then
        do n = 1, NXSC
          if ( a_XSCP(n)%IMol == m_MolNum(m) ) then
            write( 6, '(a)' ) &
              & '  XSC File : ' // trim( a_XSCFileName(n) )
            write( 6, '(a,i8,a,i8,a)' ) &
              & '             ', &
              & int( a_XSCP(n)%w_WaveNum(1) ), '-', &
              & int( a_XSCP(n)%w_WaveNum(a_XSCP(n)%NWaveNum) ), ' m-1'
          end if
        end do
      end if


      write( 6, '(a)' ) '    ---------------------------------------------'
      do k = 1, NLev
        if( mod( k, max(NLev/10,1) ) == 0 ) then
          write( 6, '(2i6,a,i5,a)' ) k, NLev, &
            & "  ---  ", int( real( k ) / NLev * 100 ), "%"
        end if

        if ( m_FlagCalcLine(m) ) then
          !
          ! pressure and temperature correction of line parameters
          !
          call lbl_lineparam_correct( &
            z_Press(k), z_Press(k) * zm_VMR(k,m), z_Temp(k), &
            DelWaveNumMod, &
            LParam % IMol         , &
            LParam % NLine        , &
            LParam % iiso         , &
            LParam % wnlc0        , &
            LParam % stren0       , &
            LParam % lowene       , &
            LParam % abwid0       , &
            LParam % sbwid0       , &
            LParam % tdep_abw     , &
            LParam % tdep_sbw     , &
            LParam % apshift      , &
            LParam % wnco1    (:,1), &
            LParam % abwidfac1(:,1), &
!!$        LParam % wnlc1   (:,k), &
!!$        LParam % stren1  (:,k), &
!!$        LParam % pbwid1  (:,k), &
!!$        LParam % dowid1  (:,k), &
!!$        LParam % renfac  (:,k), &
            LParam % wnlc1   (:,1), &
            LParam % stren1  (:,1), &
            LParam % pbwid1  (:,1), &
            LParam % dowid1  (:,1), &
            LParam % abwid1  (:,1), &
            LParam % sbwid1  (:,1), &
            LParam % renfac1 (:,1), &
            1, LParam % NLine &
            & )

          call lbl_lineparam_sort( LParam )

        end if

        !
        ! loop for wavenumber
        !

        do iWaveNum = 1, NWaveNum, NWaveNumSubLoop
!!$        if( mod( iWaveNum, NWaveNum/NWaveNumSubLoop/10 ) == 0 ) then
!!$          write( 6, '(2i6,a,i5,a)' ) iWaveNum-1, NWaveNum, &
!!$            & "  ---  ", int( real( iWaveNum-1 ) / NWaveNum * 100 ), "%"
!!$        end if

          iwnsls = 1
          iwnsle = min( NWaveNum - ( iWaveNum - 1 ), NWaveNumSubLoop )

          do iwnsl = iwnsls, iwnsle
            w_WaveNum( iwnsl ) = &
              & WaveNumS + DelWaveNumMod * dble( iWaveNum - 1 + iwnsl - 1 )
          end do


          if ( m_FlagCalcLine(m) ) then
!!$            call lbl_calc_ac( &
!!$              & z_Temp(k), LParam, &
!!$              & NWaveNumSubLoop, iwnsls, iwnsle, w_WaveNum, w_AbsCoefLine &
!!$              & )
            call lbl_calc_ac( &
              & z_Temp(k), LParam, &
              & (iwnsle-iwnsls+1), w_WaveNum(iwnsls:iwnsle), w_AbsCoefLine(iwnsls:iwnsle) &
              & )
          else
            w_AbsCoefLine = 0.0d0
          end if
!!$        call ac_io_output( m, k, NWaveNumSubLoop, w_WaveNum, w_AbsCoefLine, iwnsls, iwnsle, iWaveNum )



          if ( m_FlagCalcCont(m) ) then

            PathLength = 1.0d0
            NumDen = z_Press(k) / ( Boltz * z_Temp(k) )
!!$            NumDen = NumDen * ( 1.0d0 - zm_VMR(k,1) ) * zm_VMR(k,m)
            NumDen = NumDen * zm_VMR(k,m)
            ColDen = NumDen * PathLength

            if ( ( NumDen > 0.0d0 ) .and. &
              &  ( ( m_MolNum(m) == 1 ) .or. ( m_MolNum(m) ==  2 ) .or. &
              &    ( m_MolNum(m) == 3 ) .or. &
              &    ( m_MolNum(m) == 7 ) .or. ( m_MolNum(m) == 22 ) ) ) then

              VMRH2O = 0.0d0
              VMRCO2 = 0.0d0
              VMRO3  = 0.0d0
              VMRO2  = 0.0d0
              VMRN2  = 0.0d0
              do m2 = 1, NMolNum
                select case ( m_MolNum(m2) )
                case ( 1 )
                  VMRH2O = zm_VMR(k,m2)
                case ( 2 )
                  VMRCO2 = zm_VMR(k,m2)
                case ( 3 )
                  VMRO3  = zm_VMR(k,m2)
                case ( 7 )
                  VMRO2  = zm_VMR(k,m2)
                case ( 22 )
                  VMRN2  = zm_VMR(k,m2)
                end select
              end do
              VMROTHERS = 1.0d0 - VMRH2O - VMRCO2 - VMRO3 - VMRO2 - VMRN2
              call call_mtckd_wrapper( &
                & z_Press(k), z_Temp(k), VMRH2O, PathLength, &
                & m_MolNum(m), &
                & VMRCO2, VMRO3, VMRO2, VMRN2, VMROTHERS, &
!!$          & WaveNumS, DelWaveNumMod, &
                & WaveNumS + DelWaveNumMod * dble( iWaveNum - 1 + iwnsls - 1 ), DelWaveNumMod, &
                & NWaveNumSubLoop, &
                & w_AbsCoefCont &
                & )
              ! IMPORTANT:
              ! This conversion does not seem to be needed based on comparison
              ! between results by this program and LBLRTM. (2017/08/06).
              !
              ! multiplying a factor to convert definition of absorption
              ! coefficient (see p.6 of Ptashnik et al. (2011).
              ! Reference
              !   Ptashnik et al., 2011, JGR, 116, D16305, doi:10.1029/2011JD015603
              !
!!$              w_AbsCoefCont = w_AbsCoefCont * ( 296.0d0 / z_Temp(k) )

              ! Abs. coef. is divided by column density of absorbing molecule
              w_AbsCoefCont = w_AbsCoefCont / ColDen

            else

              w_AbsCoefCont = 0.0d0

            end if

          else

            w_AbsCoefCont = 0.0d0

          end if


          if ( FlagCIA ) then
            w_AbsCoefCIA(iwnsls:iwnsle) = 0.0d0
            do n = 1, NCIA
              if ( a_CParam(n)%IMol == m_MolNum(m) ) then
                IndexMolOpp = CIAGetMolNumIndex( NMolNum, m_MolNum, a_CParam(n)%IMolOpp )
                call CIACalcAC( &
                  & a_CParam(n), zm_VMR(k,IndexMolOpp), z_Press(k) / ( Boltz * z_Temp(k) ), z_Temp(k), &
                  & (iwnsle-iwnsls+1), w_WaveNum(iwnsls:iwnsle),&
                  & w_AbsCoefCIAEach(iwnsls:iwnsle) &
                  & )
                w_AbsCoefCIA(iwnsls:iwnsle) = &
                  &   w_AbsCoefCIA    (iwnsls:iwnsle) &
                  & + w_AbsCoefCIAEach(iwnsls:iwnsle)
              end if
            end do
            w_AbsCoefCont(iwnsls:iwnsle) = &
              & w_AbsCoefCont(iwnsls:iwnsle) + w_AbsCoefCIA(iwnsls:iwnsle)
          end if

          if ( FlagXSC ) then
            w_AbsCoefXSC(iwnsls:iwnsle) = 0.0d0
            do n = 1, NXSC
              if ( a_XSCP(n)%IMol == m_MolNum(m) ) then
                call XSCCalcAC( &
                  & a_XSCP(n), z_Temp(k), &
                  & (iwnsle-iwnsls+1), w_WaveNum(iwnsls:iwnsle),&
                  & w_AbsCoefXSCEach(iwnsls:iwnsle) &
                  & )
                w_AbsCoefXSC(iwnsls:iwnsle) = &
                  &   w_AbsCoefXSC    (iwnsls:iwnsle) &
                  & + w_AbsCoefXSCEach(iwnsls:iwnsle)
              end if
            end do
            w_AbsCoefCont(iwnsls:iwnsle) = &
              & w_AbsCoefCont(iwnsls:iwnsle) + w_AbsCoefXSC(iwnsls:iwnsle)
          end if

          if ( ( m == 1 ) .and. ( k == 1 ) ) then
            FlagOutputWaveNum = .true.
          else
            FlagOutputWaveNum = .false.
          end if

!!$          call ac_io_output( m, k, NWaveNumSubLoop, &
!!$            & "AbsCoefLine", w_WaveNum, w_AbsCoefLine, &
!!$            & iwnsls, iwnsle, iWaveNum )
!!$
!!$          call ac_io_output( m, k, NWaveNumSubLoop, &
!!$            & "AbsCoefCont", w_WaveNum, w_AbsCoefCont, &
!!$            & iwnsls, iwnsle, iWaveNum )

          call ac_io_output( m, k, NWaveNumSubLoop, &
            & "AbsCoef", w_WaveNum, w_AbsCoefLine+w_AbsCoefCont, &
            & iwnsls, iwnsle, iWaveNum, FlagOutputWaveNum )

        end do

      end do

      if ( m_FlagCalcLine(m) ) then
        call lbl_lineparam_finalize( LParam )
      end if
      if ( FlagCIA ) then
        do n = 1, NCIA
          if ( a_CParam(n)%IMol == m_MolNum(m) ) then
            call CIAParamFinalize( a_CParam(n) )
          end if
        end do
      end if
      if ( FlagXSC ) then
        do n = 1, NXSC
          if ( a_XSCP(n)%IMol == m_MolNum(m) ) then
            call XSCParamFinalize( a_XSCP(n) )
          end if
        end do
      end if

    end do

    call ac_io_output_end

    deallocate( w_AbsCoefLine    )
    deallocate( w_AbsCoefCont    )
    deallocate( w_AbsCoefCIA     )
    deallocate( w_AbsCoefCIAEach )
    deallocate( w_AbsCoefXSC     )
    deallocate( w_AbsCoefXSCEach )

  end subroutine ACCalc

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

end module ac_calc
