  module lbl_module

    use vtype_module

    implicit none

    private

    public :: LBLInit
    public :: LineParam
    public :: lbl_lineparam_init
    public :: lbl_lineparam_finalize
    public :: lbl_readfile
    public :: lbl_readfile_hitran04
    public :: lbl_writefile_netcdf
    public :: lbl_readfile_netcdf
    public :: lbl_select_lines
    public :: lbl_show_statistics
    public :: lbl_lineparam_correct
    public :: lbl_lineparam_sort
    public :: lbl_calc_ac


    ! IMol            : molecular number for specie
    ! nline           : number of lines
    ! iiso            : isotopomer index
    ! wnlc0(nline)    : line center wave number (m^-1)
    ! stren0(nline)   : line strength (m^-1/(molecule m^-2) @ 296K)
    ! lowene(nline)   : lower state energy (m^-1)
    ! tdep_abw(nline) : coefficient of temperature dependence of
    !                 : air-broadened halfwidth
    ! tdep_sbw(nline) : coefficient of temperature dependence of
    !                 : self-broadened halfwidth
    ! abwid0(nline)   : air-broadened half width (m^-1 @ 296K)
    ! sbwid0(nline)   : self-broadened half width (m^-1 @ 296K)
    ! apshift(nline)  : air pressure-induced line shift (m^-1 atm^-1 @ 296K)

    ! wnlc1(linen)    : wave number at line center after temperature/pressure
    !                 : correction (m^-1) 
    ! stren1(nline)     : line strength after temperature/pressure correction
    !                   : (m^-1 / molecule m^-2)
    ! pbwid1(nline)     : pressure broadened (lorentz) halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! dowid1(nline)     : doppler halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! abwid1(nline)     : air-broadened half width
    !                   : after temperature/pressure correction (m^-1)
    ! sbwid1(nline)     : self-broadened half width
    !                   : after temperature/pressure correction (m^-1)

    ! wnco(nline)       : cutoff wave number (m^-1)
    ! abwidfac(nline)   : air-broadened half width factor
    ! renfac(nline)     : renormalization factor

    type :: LineParam
      integer(i4b)          :: IMol
      integer(i4b)          :: nline
      integer(i4b), pointer :: iiso( : )
      real(dp)    , pointer :: wnlc0( : ), stren0( : ), lowene( : ), &
        & tdep_abw( : ), tdep_sbw( : ), apshift( : ), &
        & abwid0( : ), sbwid0( : )

      real(dp)              :: mmass

      integer(i4b)          :: km
      real(dp)    , pointer :: &
        & wnlc1( :, : ), stren1( :, : ), pbwid1( :, : ), dowid1( :, : ), &
        & abwid1( :, : ), sbwid1( :, : )
      real(dp)    , pointer :: wnco1    ( :, : )
      real(dp)    , pointer :: abwidfac1( :, : )
      real(dp)    , pointer :: renfac1  ( :, : )
    end type LineParam


    real(dp), parameter         :: LowEneMissingValue = -1.0d100

    integer              , save :: NLSFTblX
    integer              , save :: NLSFTblY
    real(dp), allocatable, save :: a_LSFTblX(:)
    real(dp), allocatable, save :: a_LSFTblY(:)
    real(dp), allocatable, save :: aa_LSFTblLSF(:,:)
    real(dp)             , save :: DelLog10LSFTblX
    real(dp)             , save :: Log10LSFTblX1
    real(dp)             , save :: DelLog10LSFTblY
    real(dp)             , save :: Log10LSFTblY1



    ! Parameters for sub-Lorentian
    ! see Perrin and Hartmann, 1989, JQSRT, 42, 311.
    !
    integer, parameter :: SLPH89NMol         = 2
    integer, parameter :: SLPH89NWaveNumBnds = 3
    integer, parameter :: SLPH89NB           = 3
    real(dp), save     :: wm_SLPH89WaveNumBnds(SLPH89NWaveNumBnds,SLPH89NMol)
    real(dp), save     :: bm_SLPH89Alpha      (SLPH89NB,SLPH89NMol)
    real(dp), save     :: bm_SLPH89Beta       (SLPH89NB,SLPH89NMol)
    real(dp), save     :: bm_SLPH89Eps        (SLPH89NB,SLPH89NMol)


    ! Unit of wm_SLPH89WaveNumBnds is cm-1
    data wm_SLPH89WaveNumBnds / 3.0d0, 30.0d0, 120.0d0, & ! CO2-CO2
      &                         3.0d0, 10.0d0,  70.0d0 /  ! CO2-N2
    data bm_SLPH89Alpha       /  0.0888d0 , 0.000d0  , 0.0232d0, & ! CO2-CO2
      &                          0.416d0  , 0.00167d0, 0.0200d0 /  ! CO2-N2
    data bm_SLPH89Beta        / -0.160d0  , 0.0526d0 , 0.0d0   , & ! CO2-CO2
      &                         -0.354d0  , 0.0421d0 , 0.0d0    /  ! CO2-N2
    data bm_SLPH89Eps         /  0.00410d0, 0.00152d0, 0.0d0   , & ! CO2-CO2
      &                          0.00386d0, 0.00248d0, 0.0d0    /  ! CO2-N2


    !**************************************************************************

  contains

    !**************************************************************************

    subroutine LBLInit

      call LBLVoigtH82TblInit

    end subroutine LBLInit

    !**************************************************************************

    subroutine lbl_lineparam_init( IMol, nline, km, lp )

      integer(i4b)     , intent(in   ) :: IMol, nline, km
      type( LineParam ), intent(inout) :: lp


      lp % IMol  = IMol
      lp % nline = nline

      allocate( &
           lp % iiso    ( nline ), &
           lp % wnlc0   ( nline ), &
           lp % stren0  ( nline ), &
           lp % lowene  ( nline ), &
           lp % tdep_abw( nline ), &
           lp % tdep_sbw( nline ), &
           lp % apshift ( nline ), &
           lp % abwid0  ( nline ), &
           lp % sbwid0  ( nline )  &
           )

      lp % km    = km

      allocate( &
           lp % wnlc1 ( nline, km ), &
           lp % stren1( nline, km ), &
           lp % pbwid1( nline, km ), &
           lp % dowid1( nline, km ), &
           lp % abwid1( nline, km ), &
           lp % sbwid1( nline, km ) &
           )

      allocate( &
           lp % wnco1    ( nline, km )    , &
           lp % abwidfac1( nline, km )    , &
           lp % renfac1  ( nline, km )  &
           )

    end subroutine lbl_lineparam_init

    !**************************************************************************

    subroutine lbl_lineparam_finalize( lp )

      type( LineParam ), intent(inout) :: lp


      lp % IMol  = -1
      lp % nline = -1

      deallocate( &
           lp % iiso, &
           lp % wnlc0, &
           lp % stren0, &
           lp % lowene, &
           lp % tdep_abw, &
           lp % tdep_sbw, &
           lp % apshift, &
           lp % abwid0, &
           lp % sbwid0 &
           )

      lp % km    = -1

      deallocate( &
           lp % wnlc1, &
           lp % stren1, &
           lp % pbwid1, &
           lp % dowid1, &
           lp % abwid1, &
           lp % sbwid1 &
           )

      deallocate( &
           lp % wnco1    , &
           lp % abwidfac1, &
           lp % renfac1    &
           )

    end subroutine lbl_lineparam_finalize

    !**************************************************************************

!!$    subroutine lbl_molmn( &
!!$      & MolNum, &
!!$      & MolMassNum &
!!$      & )
!!$
!!$      integer , intent(in ) :: MolNum
!!$      real(DP), intent(out) :: MolMassNum
!!$
!!$
!!$      ! Local variables
!!$      !
!!$      integer, parameter :: HITTblNMol    = 3
!!$      integer, parameter :: HITTblNIsoMax = 10
!!$      integer            :: im_HITTblIsoMolMN(HITTblNIsoMax, HITTblNMol)
!!$      real(DP)           :: im_HITTblIsoFrac (HITTblNIsoMax, HITTblNMol)
!!$      real(DP)           :: m_HITTblMolMN(HITTblNMol)
!!$
!!$      real(DP)           :: FracSum
!!$
!!$      integer :: i
!!$      integer :: m
!!$
!!$
!!$      ! see https://www.cfa.harvard.edu/hitran/facts.html
!!$      !     and https://www.cfa.harvard.edu/hitran/molecules.html
!!$
!!$      ! see also http://hitran.org/media/molparam.txt
!!$
!!$      data im_HITTblIsoFrac &
!!$        & / &
!!$        ! H2O
!!$        & 9.973d-1, 1.999d-3, 3.719d-4, 3.107d-4, 6.230d-7, &
!!$        & 1.158d-7, 0.0d0   , 0.0d0   , 0.0d0   , 0.0d0   , &
!!$        ! CO2
!!$        & 9.842d-1, 1.106d-2, 3.947d-3, 7.340d-4, 4.434d-5, &
!!$        & 8.246d-6, 3.957d-6, 1.472d-6, 1.368d-7, 4.446d-8, &
!!$        ! O3
!!$        & 9.929d-1, 3.982d-3, 1.991d-3, 7.405d-4, 3.702d-4, &
!!$        & 0.0d0   , 0.0d0   , 0.0d0   , 0.0d0   , 0.0d0     &
!!$        & /
!!$
!!$
!!$      do m = 1, HITTblNMol
!!$        do i = 1, HITTblNIsoMax
!!$          if ( im_HITTblIsoFrac(i,m) > 0.0d0 ) then
!!$            im_HITTblIsoMolMN(i,m) = lbl_iso_mw( m, i )
!!$          else
!!$            im_HITTblIsoMolMN(i,m) = 0
!!$          end if
!!$        end do
!!$      end do
!!$
!!$
!!$      do m = 1, HITTblNMol
!!$        m_HITTblMolMN(m) = 0.0d0
!!$        FracSum = 0.0d0
!!$        do i = 1, HITTblNIsoMax
!!$          m_HITTblMolMN(m) = m_HITTblMolMN(m) &
!!$            & + im_HITTblIsoMolMN(i,m) * im_HITTblIsoFrac(i,m)
!!$          FracSum = FracSum + im_HITTblIsoFrac(i,m)
!!$        end do
!!$        m_HITTblMolMN(m) = m_HITTblMolMN(m) / FracSum
!!$      end do
!!$
!!$      MolMassNum = m_HITTblMolMN(MolNum)
!!$
!!$
!!$    end subroutine lbl_molmn

    !**************************************************************************
    ! subroutine readfile
    ! read molecular parameters (line strength, halfwidth, etc) from file
    !**************************************************************************
    ! linen           : number of line
    ! wnlc0(linen)    : wave number at line center (m^-1)
    ! stren0(linen)   : line strength (m^-1 / (molecule m^-2))
    ! abwid0(linen)   : air-broadened half width (m^-1 @ 296K)
    ! sbwid0(linen)   : self-broadened half width (m^-1 @ 296K)
    ! lowene(linen)   : lower level energy of the transition (m^-1)
    ! tdep_abw(linen) : coefficient of temperature dependence of
    !                 : air-broadened halfwidth
    ! tdep_sbw(linen) : coefficient of temperature dependence of
    !                 : self-broadened halfwidth
    ! fn              : input filename
    !**************************************************************************
    ! CAUTION: these variables (wnlc0,lstren,ahawid,etc) are returned
    ! in MKS units (these are in CGS unit in the input file)
    !**************************************************************************

    subroutine lbl_readfile( fn, lp )

      use fi_module

      character(len=*) , intent(in   ) :: fn
      type( LineParam ), intent(inout) :: lp


      !
      ! local variables
      !
      integer(i4b)     :: mol
      real(dp)         :: trprob
      integer(i4b)     :: gqiup, gqilo
      character(len=9) :: lqup, lqlo
      integer(i4b)     :: ier
      integer(i4b)     :: iref

      integer(i4b)     :: fu
      integer(i4b)     :: l


      !
      ! input from file
      !
      call fi_open( fn, "read", fu )

      do l = 1, lp % nline
         read( fu, 1000 ) mol, &
              lp % iiso    ( l ), &
              lp % wnlc0   ( l ), &
              lp % stren0  ( l ), &
              trprob            , &
              lp % abwid0  ( l ), &
              lp % sbwid0  ( l ), &
              lp % lowene  ( l ), &
              lp % tdep_abw( l ), &
              lp % apshift ( l ), &
              gqiup, gqilo, lqup, lqlo, ier, iref

         ! check
         if ( l == 1 ) then
           if ( lp % imol /= mol ) then
             write( 6, * ) 'In lbl_readfile: '
             write( 6, * ) '  Molecular number in ', trim( fn ), &
               & ' does not match with that specifined in a variable, lp.'
             write( 6, * ) '  ', trim( fn ), ':', mol
             write( 6, * ) '  lp:', lp % imol
             stop
           end if
         end if

         if( lp % sbwid0( l ) .eq. 0.0 ) then
            write( 6, * ) l,"th line's self-broadened half width is zero"
            if( lp % abwid0( l ) .eq. 0.0 ) then
               write( 6, * ) &
                    l,"th line's air-broadened  half width is also zero"
            end if
            lp % sbwid0( l ) = lp % abwid0( l )
!!$            ! from reference in Hourdin, 1992
!!$            lowid( l ) = ahawid * 1.8d0
!!$         else
!!$            lowid( l ) = shawid
         end if

         lp % tdep_sbw( l ) = lp % tdep_abw( l )

         ! transform units of return values from CGS to MKS
         lp % wnlc0  ( l ) = lp % wnlc0  ( l ) * 1.0d2
         lp % stren0 ( l ) = lp % stren0 ( l ) * 1.0d-2
         lp % apshift( l ) = lp % apshift( l ) * 1.0d2
         lp % abwid0 ( l ) = lp % abwid0 ( l ) * 1.0d2
         lp % sbwid0 ( l ) = lp % sbwid0 ( l ) * 1.0d2
         if ( lp % lowene ( l ) == -1.0d0 ) then
           write( 6, * ) 'Lower state energy of ', l, 'th line is missing.'
           lp % lowene ( l ) = LowEneMissingValue
         else
           lp % lowene ( l ) = lp % lowene ( l ) * 1.0d2
         end if
      end do
      close( fu )

1000  format( i2, i1, f12.6, 2e10.3, 2f5.3, f10.4, f4.2, f8.5, 2i3, 2a9, i3, i6 )

    end subroutine lbl_readfile

    !**************************************************************************

    !**************************************************************************
    ! subroutine readfile
    ! read molecular parameters (line strength, halfwidth, etc) from file
    !**************************************************************************
    ! linen           : number of line
    ! wnlc0(linen)    : wave number at line center (m^-1)
    ! stren0(linen)   : line strength (m^-1 / (molecule m^-2))
    ! abwid0(linen)   : air-broadened half width (m^-1 @ 296K)
    ! sbwid0(linen)   : self-broadened half width (m^-1 @ 296K)
    ! lowene(linen)   : lower level energy of the transition (m^-1)
    ! tdep_abw(linen) : coefficient of temperature dependence of
    !                 : air-broadened halfwidth
    ! tdep_sbw(linen) : coefficient of temperature dependence of
    !                 : self-broadened halfwidth
    ! fn              : input filename
    !**************************************************************************
    ! CAUTION: these variables (wnlc0,lstren,ahawid,etc) are returned
    ! in MKS units (these are in CGS unit in the input file)
    !**************************************************************************

    subroutine lbl_readfile_hitran04( fn, lp )

      use fi_module

      character(len=*) , intent(in   ) :: fn
      type( LineParam ), intent(inout) :: lp


      !
      ! local variables
      !
      integer(i4b)      :: mol
      character         :: CIISO
      integer(i4b)      :: ICIISO
      real(dp)          :: eacoef
      character(len=15) :: gqup, gqlo, lqup, lqlo
      integer(i4b)      :: ier( 6 ), iref( 6 )
      character(len=1)  :: flag
      real(dp)          :: swup, swlo

      integer(i4b)      :: fu
      integer(i4b)      :: l


      !
      ! declaration of function
      !
!!$      real(dp)          :: uncertainty


      !
      ! input from file
      !
      call fi_open( fn, "read", fu )
      do l = 1, lp % nline
         read( fu, 1000 ) mol, &
!!$              lp % iiso    ( l ), &
              CIISO, &
              lp % wnlc0   ( l ), &
              lp % stren0  ( l ), &
              eacoef            , &
              lp % abwid0  ( l ), &
              lp % sbwid0  ( l ), &
              lp % lowene  ( l ), &
              lp % tdep_abw( l ), &
              lp % apshift ( l ), &
              gqup, gqlo, lqup, lqlo, &
              ier (1), ier (2), ier (3), ier (4), ier (5), ier (6), &
              iref(1), iref(2), iref(3), iref(4), iref(5), iref(6), &
              flag, swup, swlo

         ! check
         if ( l == 1 ) then
           if ( lp % imol /= mol ) then
             write( 6, * ) 'In lbl_readfile_hitran04: '
             write( 6, * ) '  Molecular number in ', trim( fn ), &
               & ' does not match with that specifined in a variable, lp.'
             write( 6, * ) '  ', trim( fn ), ':', mol
             write( 6, * ) '  lp:', lp % imol
             stop
           end if
         end if

         ICIISO = iachar( CIISO )
         if ( ICIISO < 48 ) then
           write( 6, * ) 'Unexpected isotope number: ', CIISO
           stop
         else if ( ( 48 <= ICIISO ) .and. ( ICIISO <=  57 ) ) then
           lp % iiso( l ) = ICIISO - 48
         else if ( ( 58 <= ICIISO ) .and. ( ICIISO <=  64 ) ) then
           write( 6, * ) 'Unexpected isotope number: ', CIISO
           stop
         else if ( ( 65 <= ICIISO ) .and. ( ICIISO <=  90 ) ) then
           lp % iiso( l ) = ICIISO - 65 + 11  ! 'A' = 11
         else
           write( 6, * ) 'Unexpected isotope number: ', CIISO
           stop
         end if

         ! NOTE:
         !   Isopomer index of 0 is changed to 10 for clarity.
         if( lp % iiso( l ) == 0 ) then
           lp % iiso( l ) = 10
         end if

         if( lp % sbwid0( l ) .eq. 0.0 ) then
            write( 6, * ) l,"th line's self-broadened half width is zero"
            if( lp % abwid0( l ) .eq. 0.0 ) then
               write( 6, * ) &
                    l,"th line's air-broadened  half width is also zero"
            end if
            lp % sbwid0( l ) = lp % abwid0( l )
!!$            ! from reference in Hourdin, 1992
!!$            lowid( l ) = ahawid * 1.8d0
!!$         else
!!$            lowid( l ) = shawid
         end if

         lp % tdep_sbw( l ) = lp % tdep_abw( l )

         ! transform units of return values from CGS to MKS
         lp % wnlc0   ( l ) = lp % wnlc0   ( l ) * 1.0d2
         lp % stren0  ( l ) = lp % stren0  ( l ) * 1.0d-2
!!$         lp % stren0 ( l ) = lp % stren0 ( l ) * 1.0d-2 * uncertainty( ier( 2 ) )
         lp % apshift ( l ) = lp % apshift ( l ) * 1.0d2
         lp % abwid0  ( l ) = lp % abwid0  ( l ) * 1.0d2
         lp % sbwid0  ( l ) = lp % sbwid0  ( l ) * 1.0d2
         if ( lp % lowene ( l ) == -1.0d0 ) then
           write( 6, * ) 'Lower state energy of ', l, 'th line is missing.'
           lp % lowene ( l ) = LowEneMissingValue
         else
           lp % lowene  ( l ) = lp % lowene  ( l ) * 1.0d2
         end if
      end do
      close( fu )

      ! see https://www.cfa.harvard.edu/hitran/formats.html
!1000  format( I2, I1, F12.6, E10.3, E10.3, F5.4, F5.4, F10.4, F4.2, F8.6, A15, A15, A15, A15, 6I1, 6I2, A1, F7.1, F7.1 )
1000  format( I2, A1, F12.6, E10.3, E10.3, F5.4, F5.4, F10.4, F4.2, F8.6, A15, A15, A15, A15, 6I1, 6I2, A1, F7.1, F7.1 )

    end subroutine lbl_readfile_hitran04

    !**************************************************************************

    subroutine lbl_select_lines( LineStrenMin, lp1, lp2 )


      real(DP)         , intent(in ) :: LineStrenMin
      type( LineParam ), intent(in ) :: lp1
      type( LineParam ), intent(out) :: lp2


      !
      ! local variables
      !
      integer(i4b)              :: NLine2
      logical     , allocatable :: a_FlagOutput(:)
      integer(i4b)              :: iLine
      integer(i4b)              :: l


      allocate( a_FlagOutput( lp1%NLine ) )

      NLine2 = 0
      do l = 1, lp1%NLine
!!$        if ( .true. ) then
        if ( lp1%stren0(l) >= LineStrenMin ) then
          a_FlagOutput( l ) = .true.
          NLine2 = NLine2 + 1
        else
          a_FlagOutput( l ) = .false.
        end if
      end do
      write( 6, '(a,i9)' ) " Number of selected lines: ", NLine2

      call lbl_lineparam_init( lp1%IMol, NLine2, lp1%km, lp2 )

      iLine = 0
      do l = 1, lp1%NLine
        if ( a_FlagOutput(l) ) then
          iLine = iLine + 1
          if ( iLine > NLine2 ) then
            write( 6, * ) "Unexpected number of selected lines: ", NLine2
            stop
          end if
          lp2%iiso    (iLine) = lp1%iiso    (l)
          lp2%wnlc0   (iLine) = lp1%wnlc0   (l)
          lp2%stren0  (iLine) = lp1%stren0  (l)
          lp2%lowene  (iLine) = lp1%lowene  (l)
          lp2%tdep_abw(iLine) = lp1%tdep_abw(l)
          lp2%tdep_sbw(iLine) = lp1%tdep_sbw(l)
          lp2%apshift (iLine) = lp1%apshift (l)
          lp2%abwid0  (iLine) = lp1%abwid0  (l)
          lp2%sbwid0  (iLine) = lp1%sbwid0  (l)
        end if
      end do

      if ( iLine /= NLine2 ) then
        write( 6, * ) "Number of selected lines is not correct: ", iLine, NLine2
        stop
      end if

      deallocate( a_FlagOutput )


    end subroutine lbl_select_lines

    !**************************************************************************

    subroutine lbl_writefile_netcdf( fn, Src, Ref, lp )

      use ni3_module

      character(len=*) , intent(in   ) :: fn
      character(len=*) , intent(in   ) :: Src
      character(len=*) , intent(in   ) :: Ref
      type( LineParam ), intent(inout) :: lp


      !
      ! local variables
      !
      integer(i4b)     , allocatable :: a_iLine(:)

      integer(i4b)                   :: NDims
      character(stdstr), allocatable :: a_DimNames(:)
      integer(i4b)                   :: NcID
      character(stdstr)              :: Mode
      character(stdstr)              :: Name
      character(stdstr)              :: StdName
      character(extstr)              :: LongName
      character(extstr)              :: Unit

      integer(i4b)      :: iLine
      integer(i4b)      :: l


      allocate( a_iLine(lp%NLine) )
      do l = 1, lp%NLine
        a_iLine(l) = l
      end do


      Mode = "new"
      call ni3_open( fn, Mode, NcID )

      call ni3_set_ga( NcID, src = Src, ref = Ref )

      Name     = "Line"
      StdName  = Name
      LongName = "line number"
      Unit     = "1"
      call ni3_set_dim( NcID, Name, NI3_INT, a_iLine, &
        & StdName, LongName, Unit )
!!$      call ni3_def_dim( NcID, Name, NI3_INT, NI3_UNLIMITED, &
!!$        & StdName, LongName, Unit )

      NDims = 0
      allocate( a_DimNames( NDims ) )

      ! IMol            : molecular number for specie
      Name     = "imol"
      StdName  = Name
      LongName = "line number"
      Unit     = "1"
      call ni3_def_var( NcID, Name, NI3_INT, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%IMol )

      deallocate( a_DimNames )


      NDims = 1
      allocate( a_DimNames( NDims ) )
      a_DimNames(1) = "Line"

      ! iiso            : isotopomer index
      Name     = "iiso"
      StdName  = Name
      LongName = "isotopomer index"
      Unit     = "1"
      call ni3_def_var( NcID, Name, NI3_INT, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%iiso )

      ! wnlc0(nline)    : line center wave number (m^-1)
      Name     = "wnlc0"
      StdName  = Name
      LongName = "line center wave number"
      Unit     = "m-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%wnlc0 )

      ! stren0(nline)   : line strength (m^-1/(molecule m^-2) @ 296K)
      Name     = "stren0"
      StdName  = Name
      LongName = "line strength @ 296 K"
      Unit     = "m-1 (molecule m-2)-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%stren0 )

      ! lowene(nline)   : lower state energy (m^-1)
      Name     = "lowene0"
      StdName  = Name
      LongName = "lower state energy"
      Unit     = "m-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_att( NcID, Name, 'missing_value', LowEneMissingValue )
      call ni3_put_var( NcID, Name, lp%lowene )

      ! tdep_abw(nline) : coefficient of temperature dependence of
      !                 : air-broadened halfwidth
      Name     = "tdep_abw"
      StdName  = Name
      LongName = "coefficient of temperature dependence of air-broadened halfwidth"
      Unit     = "1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%tdep_abw )

      ! tdep_sbw(nline) : coefficient of temperature dependence of
      !                 : self-broadened halfwidth
      Name     = "tdep_sbw"
      StdName  = Name
      LongName = "coefficient of temperature dependence of self-broadened halfwidth"
      Unit     = "1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%tdep_sbw )

      ! abwid0(nline)   : air-broadened half width (m^-1 @ 296K)
      Name     = "abwid0"
      StdName  = Name
      LongName = "air-broadened halfwidth @ 296 K"
      Unit     = "m-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%abwid0 )

      ! sbwid0(nline)   : self-broadened half width (m^-1 @ 296K)
      Name     = "sbwid0"
      StdName  = Name
      LongName = "self-broadened halfwidth @ 296 K"
      Unit     = "m-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%sbwid0 )

      ! apshift(nline)  : air pressure-induced line shift (m^-1 atm^-1 @ 296K)
      Name     = "apshift"
      StdName  = Name
      LongName = "air pressure-induced line shift @ 296 K"
      Unit     = "m-1 atm-1"
      call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
        & StdName, LongName, Unit )
      call ni3_put_var( NcID, Name, lp%apshift )


!!$      iLine = 1
!!$      do l = 1, lp%NLine
!!$
!!$        if ( .true. ) then
!!$          Name     = "Line"
!!$          call ni3_put_varss( NcID, Name, iLine, a_iLine(iLine) )
!!$          Name     = "iiso"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%iiso(iLine) )
!!$          Name     = "wnlc0"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%wnlc0(iLine) )
!!$          Name     = "stren0"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%stren0(iLine) )
!!$          Name     = "lowene0"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%lowene(iLine) )
!!$          Name     = "tdep_abw"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%tdep_abw(iLine) )
!!$          Name     = "tdep_sbw"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%tdep_sbw(iLine) )
!!$          Name     = "abwid0"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%abwid0(iLine) )
!!$          Name     = "sbwid0"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%sbwid0(iLine) )
!!$          Name     = "apshift"
!!$          call ni3_put_varss( NcID, Name, iLine, lp%apshift(iLine) )
!!$
!!$          iLine = iLine + 1
!!$        end if
!!$
!!$      end do


      call ni3_close( NcID )

      deallocate( a_iLine    )
      deallocate( a_DimNames )


    end subroutine lbl_writefile_netcdf

    !**************************************************************************

    subroutine lbl_readfile_netcdf( fn, km, lp )

      use ni3_module

      character(len=*) , intent(in ) :: fn
      integer          , intent(in ) :: km
      type( LineParam ), intent(out) :: lp


      !
      ! local variables
      !
      real(dp)          :: LowEneMissingValueInput
      integer(i4b)      :: NcID
      character(stdstr) :: Mode
      character(stdstr) :: Name
      integer(i4b)      :: IMol
      integer(i4b)      :: NLine
      integer(i4b)      :: l


      Mode = "read"
      call ni3_open( fn, Mode, NcID )

      ! IMol            : molecular number for specie
      Name     = "imol"
      call ni3_get_var( NcID, Name, IMol )

      Name     = "Line"
      call ni3_inq_dimlen( NcID, Name, NLine )

      call lbl_lineparam_init( IMol, NLine, km, lp )

      ! iiso            : isotopomer index
      Name     = "iiso"
      call ni3_get_var( NcID, Name, lp%iiso )

      ! wnlc0(nline)    : line center wave number (m^-1)
      Name     = "wnlc0"
      call ni3_get_var( NcID, Name, lp%wnlc0 )

      ! stren0(nline)   : line strength (m^-1/(molecule m^-2) @ 296K)
      Name     = "stren0"
      call ni3_get_var( NcID, Name, lp%stren0 )

      ! lowene(nline)   : lower state energy (m^-1)
      Name     = "lowene0"
      call ni3_get_att( NcID, Name, 'missing_value', LowEneMissingValueInput )
      if ( LowEneMissingValue /= LowEneMissingValueInput ) then
        write( 6, * ) &
          & 'A missing value for lower state energy in file does not match with that defined in a module.'
        write( 6, * ) LowEneMissingValueInput, LowEneMissingValue
        stop
      end if
      call ni3_get_var( NcID, Name, lp%lowene )

      ! tdep_abw(nline) : coefficient of temperature dependence of
      !                 : air-broadened halfwidth
      Name     = "tdep_abw"
      call ni3_get_var( NcID, Name, lp%tdep_abw )

      ! tdep_sbw(nline) : coefficient of temperature dependence of
      !                 : self-broadened halfwidth
      Name     = "tdep_sbw"
      call ni3_get_var( NcID, Name, lp%tdep_sbw )

      ! abwid0(nline)   : air-broadened half width (m^-1 @ 296K)
      Name     = "abwid0"
      call ni3_get_var( NcID, Name, lp%abwid0 )

      ! sbwid0(nline)   : self-broadened half width (m^-1 @ 296K)
      Name     = "sbwid0"
      call ni3_get_var( NcID, Name, lp%sbwid0 )

      ! apshift(nline)  : air pressure-induced line shift (m^-1 atm^-1 @ 296K)
      Name     = "apshift"
      call ni3_get_var( NcID, Name, lp%apshift )

      call ni3_close( NcID )


    end subroutine lbl_readfile_netcdf

    !**************************************************************************

    subroutine lbl_show_statistics( lp )


      type( LineParam ), intent(in ) :: lp


      !
      ! local variables
      !
      integer(i4b)              :: NLine
      integer(i4b), allocatable :: a_strenOrder(:)
      integer(i4b), allocatable :: a_strenHist (:)
      real(dp)                  :: strenMin
      real(dp)                  :: strenMax
      integer(i4b)              :: OrderMin
      integer(i4b)              :: OrderMax
      integer(i4b)              :: i
      integer(i4b)              :: l


      l = 1
      strenMax = lp%stren0(l)
      strenMin = strenMax
      do l = 1, lp%NLine
        strenMax = max( strenMax, lp%stren0(l) )
        strenMin = min( strenMin, lp%stren0(l) )
      end do

      OrderMin = nint( log10( strenMin ) )
      OrderMax = nint( log10( strenMax ) )

      allocate( a_strenOrder(OrderMin:OrderMax) )
      allocate( a_strenHist (OrderMin:OrderMax) )

      do i = OrderMin, OrderMax
        a_strenOrder(i) = i
      end do

      a_strenHist = 0
      do l = 1, lp%NLine
        i = nint( log10( lp%stren0(l) ) )
        a_strenHist( i ) = a_strenHist( i ) + 1
      end do

      NLine = 0
      do i = OrderMin, OrderMax
        NLine = NLine + a_strenHist(i)
      end do

      if ( NLine /= lp%NLine ) then
        write( 6, * ) "Histgram does not include all lines: ", NLine, lp%NLine
        stop
      end if

      write( 6, '(a,i9)' ) " Statistics : Number of lines: ", NLine
      do i = OrderMin, OrderMax
        write( 6, '(a,i4,a,i9)' ) "    ", i, " : ", a_strenHist(i)
      end do

    end subroutine lbl_show_statistics

    !**************************************************************************

    function uncertainty( ierr ) result( derr )


      integer(i4b), intent(in) :: ierr

      real(dp)                 :: derr


      derr = 0.0d0

      select case( ierr )
      case( 0 )
         derr = 0.0d0
      case( 1 )
         derr = 0.0d0
      case( 2 )
         derr = 0.0d0
      case( 3 )
         derr = 0.20d0
      case( 4 )
         derr = 0.10d0
      case( 5 )
         derr = 0.05d0
      case( 6 )
         derr = 0.02d0
      case( 7 )
         derr = 0.01d0
      case( 8 )
         derr = 0.0d0
      case default
         write( 6, * ) 'Unexpected uncertainty code : ', ierr
         stop
      end select

      derr = 1.0d0 + derr


    end function uncertainty




    !**************************************************************************
    ! subroutine lbl_lineparam_correct
    ! calculate line strength, lorentz halfwidth, doppler halfwidth
    !**************************************************************************
    ! tpress            : total pressure (Pa)
    ! ppress            : partial pressure (Pa)
    ! temp              : temperature (K)
    ! mmass             : molecular mass (kg)
    ! IMol              : molecular number
    ! linen             : number of lines
    ! iiso              : isotopomer index
    ! wnlc0(linen)      : wave number at line center (m^-1)
    ! stren0(linen)     : line strength (m^-1 / molecule m^-2)
    ! lowene(linen)     : lower level energy of the transition (m^-1)
    ! abwid0(linen)     : air-broadened half width (m^-1 @ 296K)
    ! sbwid0(linen)     : self-broadened half width (m^-1 @ 296K)
    ! tdep_abw(linen)   : coefficient of temperature dependence of
    !                   : air-broadened halfwidth
    ! tdep_sbw(linen)   : coefficient of temperature dependence of
    !                   : self-broadened halfwidth
    ! apshift(nline)    : air pressure-induced line shift (m^-1 atm^-1 @ 296K)
    ! wnlc1(linen)      : wave number at line center after temperature/pressure
    !                   : correction (m^-1) 
    ! stren1(linen)     : line strength after temperature/pressure correction 
    !                   : (m^-1 / molecule m^-2)
    ! lowid1(linen)     : lorentz halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! dowid1(linen)     : doppler halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! kmingn
    ! kmaxgn
    !**************************************************************************

    subroutine lbl_lineparam_correct( &
         tpress, ppress, temp, dwn, &
         IMol, linen, iiso, wnlc0, stren0, lowene, &
         abwid0, sbwid0, tdep_abw, tdep_sbw, apshift, &
         & wnco, abwidfac, &
         wnlc1, stren1, pbwid1, dowid1, abwid1, sbwid1, renfac1, &
         ls, le )

      use constants0, only : &
        & c=>SpeedOfLight, &
        & Boltz, &
        & Planc=>PlanckConst, &
!!$        & amu=>AtomMassConst, &
        & AvogConst

      use hitranconst, only : &
        & NHITRANMolNum, &
        & HITRANConstIsoMolWt

      real(dp)    , intent(in ) :: tpress, ppress, temp, dwn
      integer(i4b), intent(in ) :: IMol
      integer(i4b), intent(in ) :: linen
      integer(i4b), intent(in ) :: iiso( linen )
      real(dp)    , intent(in ) :: wnlc0( linen ), stren0( linen ), &
           lowene( linen ), abwid0( linen ), sbwid0( linen ), &
           tdep_abw( linen ), tdep_sbw( linen ), &
           apshift( linen )
      real(dp)    , intent(in ) :: wnco    ( linen )
      real(dp)    , intent(in ) :: abwidfac( linen )
      real(dp)    , intent(out) :: &
           wnlc1( linen ), stren1( linen ), pbwid1( linen ), &
           dowid1( linen ), abwid1( linen ), sbwid1( linen ), renfac1( linen )
      integer(i4b), intent(in ) :: ls, le


      !
      ! local variables
      !
      real(dp)                  :: press0 = 101325.0d0, temp0 = 296.0d0
      real(dp)                  :: c2 = planc * c / boltz
      real(dp)    , allocatable :: tips0( : ), tips( : )
      real(dp)    , allocatable :: mmass( : )
      integer(i4b)              :: niso
      real(dp)                  :: gi
      integer(i4b)              :: i, l


      niso = 1
      do l = ls, le
         niso = max( iiso( l ), niso )
      end do
      allocate( tips0( niso ), tips( niso ) )
      do i = 1, niso
!!$         call BD_TIPS_2003( IMol, temp0, i, gi, tips0( i ) )
!!$         call BD_TIPS_2003( IMol, temp , i, gi, tips ( i ) )
        !
!!$        call BD_TIPS_2011( &
!!$          & IMol, temp0, i, &
!!$          & gi, tips0( i ) &
!!$          & )
!!$        call BD_TIPS_2011( &
!!$          & IMol, temp , i, &
!!$          & gi, tips ( i ) &
!!$          & )
        !
        call BD_TIPS_2017( &
          & IMol, temp0, i, &
          & gi, tips0( i ) &
          & )
        call BD_TIPS_2017( &
          & IMol, temp , i, &
          & gi, tips ( i ) &
          & )
      end do
!!$      tips0( : ) = 1.0d0
!!$      tips ( : ) = 1.0d0


      allocate( mmass( niso ) )
      do i = 1, niso
!!$         mmass( i ) = lbl_iso_mw( IMol, i ) * amu
!!$        mmass( i ) = lbl_iso_mw( IMol, i ) / AvogConst
        mmass( i ) = HITRANConstIsoMolWt( IMol, i ) / AvogConst
      end do

!!$!      stren1( : ) = stren0( : ) * ( temp0 / temp ) &
!!$!           * exp( planc * c / boltz * lowene( : ) &
!!$!           * ( 1.0d0 / temp0 - 1.0d0 / temp ) )
!!$      stren1( : ) = stren0( : ) &
!!$           * tips0( iiso( : ) ) / tips( iiso( : ) ) &
!!$           * exp( c2 * lowene( : ) * ( 1.0d0 / temp0 - 1.0d0 / temp ) ) &
!!$           * ( 1.0d0 - exp( -c2 * wnlc0( : ) / temp  ) ) &
!!$           / ( 1.0d0 - exp( -c2 * wnlc0( : ) / temp0 ) )
!!$
!!$      pbwid1( : ) &
!!$           = abwid0(:)*((tpress-ppress)/press0) * (temp0/temp)**tdep_abw(:) &
!!$           + sbwid0(:)*(ppress         /press0) * (temp0/temp)**tdep_sbw(:)
!!$
!!$      dowid1( : ) = wnlc0( : ) / c * sqrt( 2.0d0 * boltz * temp / mmass )



      do l = 1, ls-1
         wnlc1 ( l ) = 1.0d100
         stren1( l ) = 1.0d100
         pbwid1( l ) = 1.0d100      
         dowid1( l ) = 1.0d100
         abwid1( l ) = 1.0d100      
         sbwid1( l ) = 1.0d100      
      end do
      do l = ls, le
!!$         wnlc1 ( l ) = wnlc0( l )
         wnlc1 ( l ) = wnlc0( l ) + apshift( l ) * ( tpress / press0 )

         if ( lowene(l) == LowEneMissingValue ) then
           stren1( l ) = stren0( l )
         else
           stren1( l ) = stren0( l ) &
             * tips0( iiso( l ) ) / tips( iiso( l ) ) &
             * exp( c2 * lowene( l ) * ( 1.0d0 / temp0 - 1.0d0 / temp ) ) &
             * ( 1.0d0 - exp( -c2 * wnlc1( l ) / temp  ) ) &
             / ( 1.0d0 - exp( -c2 * wnlc1( l ) / temp0 ) )
         end if

!!$         pbwid1( l ) &
!!$              = abwid0(l)*((tpress-ppress)/press0)*(temp0/temp)**tdep_abw(l) &
!!$              + sbwid0(l)*(ppress         /press0)*(temp0/temp)**tdep_sbw(l)
         !
         abwid1( l ) &
              = abwid0(l)*((tpress-ppress)/press0)*(temp0/temp)**tdep_abw(l)
         abwid1( l ) = abwid1( l ) * abwidfac( l )
         sbwid1( l ) &
              = sbwid0(l)*(ppress         /press0)*(temp0/temp)**tdep_sbw(l)
         pbwid1( l ) = abwid1(l) + sbwid1(l)

         dowid1( l ) = wnlc1( l ) / c &
              * sqrt( 2.0d0 * boltz * temp / mmass( iiso( l ) ) )
      end do
      do l = le+1, linen
         wnlc1 ( l ) = 1.0d100
         stren1( l ) = 1.0d100
         pbwid1( l ) = 1.0d100
         dowid1( l ) = 1.0d100
         abwid1( l ) = 1.0d100
         sbwid1( l ) = 1.0d100
      end do


      deallocate( tips0, tips )

      deallocate( mmass )


!!$      do l = 1, linen
!!$         renfac( l ) &
!!$              = lbl_lineshape_renorm( pbwid1( l ), dowid1( l ), wnco( l ), dowid1( l ) / 5.0d0 )
!!$!              = lbl_lineshape_renorm( pbwid1( l ), dowid1( l ), wnco( l ), dwn )
!!$
!!$!         renfac( l ) = lbl_lineshape_voigt_for_normalize &
!!$!              ( pbwid1( l ), dowid1( l ), wnco( l ), dwn )
!!$      end do

      renfac1( : ) = 1.0d0


    end subroutine lbl_lineparam_correct

    !**************************************************************************

    subroutine lbl_lineparam_sort( lp )

      use sort_module, only : sort_heap

      type( LineParam ), intent(inout) :: lp


      ! Local variables
      !
      type( LineParam ) :: lpBK
      integer           :: aa_INdex(lp%nline,lp%km)
      integer           :: k
      integer           :: l


      do k = 1, lp%km
        call sort_heap( lp%nline, lp%wnlc1(:,k), aa_Index(:,k) )
      end do

      call lbl_lineparam_init( lp%IMol, lp%nline, lp%km, lpBK )
      lpBK%stren1    = lp%stren1
      lpBK%pbwid1    = lp%pbwid1
      lpBK%dowid1    = lp%dowid1
      lpBK%abwid1    = lp%abwid1
      lpBK%sbwid1    = lp%sbwid1
      lpBK%wnco1     = lp%wnco1
      lpBK%abwidfac1 = lp%abwidfac1
      lpBK%renfac1   = lp%renfac1
      do k = 1, lp%km
        do l = 1, lp%nline
          lp%stren1 (l,k) = lpBK%stren1 (aa_Index(l,k),k)
          lp%pbwid1 (l,k) = lpBK%pbwid1 (aa_Index(l,k),k)
          lp%dowid1 (l,k) = lpBK%dowid1 (aa_Index(l,k),k)
          lp%abwid1 (l,k) = lpBK%abwid1 (aa_Index(l,k),k)
          lp%sbwid1 (l,k) = lpBK%sbwid1 (aa_Index(l,k),k)
        end do
      end do
      do k = 1, lp%km
        do l = 1, lp%nline
          lp%wnco1    (l,k) = lpBK%wnco1    (aa_Index(l,k),k)
          lp%abwidfac1(l,k) = lpBK%abwidfac1(aa_Index(l,k),k)
          lp%renfac1  (l,k) = lpBK%renfac1  (aa_Index(l,k),k)
        end do
      end do
      call lbl_lineparam_finalize( lpBK )


    end subroutine lbl_lineparam_sort

    !**************************************************************************

!!$    function OLD_lbl_iso_mw( IMol, iiso ) result( mw )
!!$
!!$      integer(i4b), intent(in ) :: IMol
!!$      integer(i4b), intent(in ) :: iiso
!!$
!!$      integer(i4b)              :: mw
!!$
!!$
!!$      ! see http://hitran.org/docs/iso-meta/
!!$
!!$      select case( IMol )
!!$      case( 1 ) ! H2O
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 1 + 16 + 1
!!$         case( 2 ) 
!!$            mw = 1 + 18 + 1
!!$         case( 3 ) 
!!$            mw = 1 + 17 + 1
!!$         case( 4 ) 
!!$            mw = 1 + 16 + 2
!!$         case( 5 ) 
!!$            mw = 1 + 18 + 2
!!$         case( 6 ) 
!!$            mw = 1 + 17 + 2
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 2 ) ! CO2
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 16 + 12 + 16
!!$         case( 2 ) 
!!$            mw = 16 + 13 + 16
!!$         case( 3 ) 
!!$            mw = 16 + 12 + 18
!!$         case( 4 ) 
!!$            mw = 16 + 12 + 17
!!$         case( 5 ) 
!!$            mw = 16 + 13 + 18
!!$         case( 6 ) 
!!$            mw = 16 + 13 + 17
!!$         case( 7 ) 
!!$            mw = 18 + 12 + 18
!!$         case( 8 ) 
!!$            mw = 18 + 12 + 17
!!$         case( 9 ) 
!!$            mw = 17 + 12 + 17
!!$          case( 10 )
!!$            ! NOTE:
!!$            !   iiso: Isotopomer index
!!$            !   iiso for this isotope is zero in HITRAN database
!!$            !   iiso is modified for clarity in a subroutine reading a database
!!$            mw = 18 + 13 + 18
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 3 ) ! O3
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 16 + 16 + 16
!!$         case( 2 ) 
!!$            mw = 16 + 16 + 18
!!$         case( 3 ) 
!!$            mw = 16 + 18 + 16
!!$         case( 4 ) 
!!$            mw = 16 + 16 + 17
!!$         case( 5 ) 
!!$            mw = 16 + 17 + 16
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 4 ) ! N2O
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 14 + 14 + 16
!!$         case( 2 ) 
!!$            mw = 14 + 15 + 16
!!$         case( 3 ) 
!!$            mw = 15 + 14 + 16
!!$         case( 4 ) 
!!$            mw = 14 + 14 + 18
!!$         case( 5 ) 
!!$            mw = 14 + 14 + 17
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 5 ) ! CO
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 12 + 16
!!$         case( 2 ) 
!!$            mw = 13 + 16
!!$         case( 3 ) 
!!$            mw = 12 + 18
!!$         case( 4 ) 
!!$            mw = 12 + 17
!!$         case( 5 ) 
!!$            mw = 13 + 18
!!$         case( 6 ) 
!!$            mw = 13 + 17
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 6 ) ! CH4
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 12 + 1 + 1 + 1 + 1
!!$         case( 2 ) 
!!$            mw = 13 + 1 + 1 + 1 + 1
!!$         case( 3 ) 
!!$            mw = 12 + 1 + 1 + 1 + 2
!!$         case( 4 ) 
!!$            mw = 13 + 1 + 1 + 1 + 2
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case( 7 ) ! O2
!!$         select case( iiso )
!!$         case( 1 ) 
!!$            mw = 16 + 16
!!$         case( 2 ) 
!!$            mw = 16 + 18
!!$         case( 3 ) 
!!$            mw = 16 + 17
!!$         case default
!!$            write( 6, * ) 'Unexpected isotope number, ', iiso, &
!!$                 ', of molecular number, ', IMol, '.'
!!$            stop
!!$         end select
!!$      case default
!!$         write( 6, * ) 'Unexpected molecular number, ', IMol, '.'
!!$         stop
!!$      end select
!!$
!!$
!!$    end function OLD_lbl_iso_mw

    !**************************************************************************

    subroutine lbl_search_indices( nline, wnlc, wn, wnco, ls, le )

      integer(i4b), intent(in ) :: nline
      real(dp)    , intent(in ) :: wnlc( nline )
      real(dp)    , intent(in ) :: wn, wnco( nline )
      integer(i4b), intent(out) :: ls, le


      !
      ! local variables
      !
      integer(i4b) :: l


      search_ls : do l = 1, nline
         if( wnlc( l ) .ge. ( wn - wnco( l ) ) ) exit search_ls
      end do search_ls
      ls = l
      !
      ! NOTE:
      ! Start index for the next loop "search_le : do l = ls, nline"
      ! must be l=ls, instead of l=ls+1. 
      ! Start index of l=ls can treat successfully the case that there
      ! are no absorption lines between wn-wnco and wn+wnco. 
      !
      search_le : do l = ls, nline
         if( wnlc( l ) .gt. ( wn + wnco( l ) ) ) exit search_le
      end do search_le
      le = l - 1


    end subroutine lbl_search_indices

    !**************************************************************************

    subroutine lbl_calc_ac( &
      & Temp, lparam, &
      & nwnsl, wn, ac &
      & )

      use constants0, only : &
        & PI, &
        & Boltz, &
        & Planc=>PlanckConst


      real(dp)         , intent(in ) :: Temp
      type( LineParam ), intent(in ) :: lparam
      integer(i4b)     , intent(in ) :: nwnsl
      real(dp)         , intent(in ) :: wn( nwnsl )
      real(dp)         , intent(out) :: ac( nwnsl )


      !
      ! local variables
      !
!!$      real(dp)     :: MolMassNum

      integer(i4b) :: iwnsls
      integer(i4b) :: iwnsle

      real(dp)     :: fredist( nwnsl ), lsf, lorwid, dopwid, dv
      integer(i4b) :: k, l
      integer(i4b) :: l_start
      integer(i4b) :: iwnsl, ls( nwnsl ), le( nwnsl )

      integer  :: l_edge1
      integer  :: l_edge2
      integer  :: l_mp

      real(dp) :: CutOffFactor
      real(dp) :: VVHNormFactor

!!$      logical            :: FlagCO2SubLorentzianFactor = .false.
      logical            :: FlagCO2SubLorentzianFactor = .true.

      real(dp)           :: AbsDelWaveNum
      real(dp)           :: bm_SLPH89B    (SLPH89NB,SLPH89NMol)
      real(dp)           :: m_SLPH89Factor(SLPH89NMol)
      integer            :: m
      !
      real(dp) :: FactorCO2CO2
      real(dp) :: FactorCO2N2

      !
      ! varialbes for inline expansion
      !
      real(dp)    :: x, y, s
      complex(dp) :: w4, t, u


      iwnsls = 1
      iwnsle = nwnsl


!!$      ! Set molecular mass number
!!$      !
!!$      call lbl_molmn( &
!!$        & lparam % IMol, &
!!$        & MolMassNum &
!!$        & )


      ! parameters for sub-Lorentian
      if ( lparam % IMol == 2 ) then
        bm_SLPH89B = bm_SLPH89Alpha + bm_SLPH89Beta * exp( - bm_SLPH89Eps * Temp )
      end if


      k = 1

      do iwnsl = iwnsls, iwnsle
        ac( iwnsl ) = 0.0d0
      end do

      do iwnsl = iwnsls, iwnsle
!!$              call lbl_search_indices( &
!!$                   lparam( n ) % nline, &
!!$                   lparam( n ) % wnlc0, &
!!$                   wn( iwnsl ), lparam( n ) % wnco, ls( iwnsl ), le( iwnsl ) )

        ! The lbl_search_indices subroutine sometimes cause very 
        ! bad computational performance. Sometimes, it does not a 
        ! matter.
        ! The reason is unknown. 

        l_start = 1
        !
!        l_edge1 = 1
        l = 1
        do
          if ( lparam%wnlc1(l,k) > 0.0d0 ) exit
          if ( l > lparam%nline ) exit
          l = l + 1
        end do
        l_edge1 = l
        !
        l_edge2 = lparam % nline
        l_mp    = ( l_edge1 + l_edge2 ) / 2
        do
          l = l_mp
          if( lparam%wnlc1(l,k) >= wn(iwnsl) - lparam%wnco1(l,k) ) then
            l_edge2 = l_mp
          else
            l_edge1 = l_mp
          end if
          if ( l_edge2 - l_edge1 <= 1 ) exit
          l_mp = ( l_edge1 + l_edge2 ) / 2
        end do
        l_start = l_edge2


        search_ls : do l = l_start, lparam % nline
          if( lparam % wnlc1( l, k ) .ge. &
            ( wn( iwnsl ) - lparam % wnco1( l, k ) ) ) &
            exit search_ls
        end do search_ls
        ls( iwnsl ) = l
        !
        ! NOTE:
        ! Start index for the next loop "search_le : do l = ls, nline"
        ! must be l=ls, instead of l=ls+1. 
        ! Start index of l=ls can treat successfully the case that there
        ! are no absorption lines between wn-wnco and wn+wnco. 
        !
        search_le : do l = ls( iwnsl ), lparam % nline
          if( lparam % wnlc1( l, k ) .gt. &
            ( wn( iwnsl ) + lparam % wnco1( l, k ) ) ) &
            exit search_le
        end do search_le
        le( iwnsl ) = l - 1

      end do

      do iwnsl = iwnsls, iwnsle
        do l = ls( iwnsl ), le( iwnsl )
          fredist( iwnsl ) = lparam % wnlc1( l, k ) - wn( iwnsl )


!!$                 lsf = &
!!$                      lbl_voigt_h82 &
!!$                      lbl_voigt_fs85 &
!!$                      ( lparam( n ) % pbwid1( l, k ), &
!!$                        lparam( n ) % dowid1( l, k ), &
!!$                        abs( fredist( iwnsl ) ) )

          lorwid = lparam % pbwid1( l, k )
          dopwid = lparam % dowid1( l, k )
          dv     = abs( fredist( iwnsl ) )


          !-----
          include 'lbl_voigt_h82.f90'
          lsf = lsf / ( dopwid * sqrt( pi ) )
          !-----
!!$          lsf = lbl_voigt_h82 ( lorwid, dopwid, dv )
!!$          lsf = lsf / ( dopwid * sqrt( pi ) )
          !-----
!!$          lsf = LBLVoigtH82CalcWithTbl( LorWid, DopWid, DV )
!!$          lsf = lsf / ( dopwid * sqrt( pi ) )
          !-----
!!$                 lsf    = lbl_voigt_fs85( lorwid, dopwid, dv )
          !-----
!!$          lsf = lbl_lorentz( lorwid, dv )
          !-----


          ! Van Vleck-Huber normalization factor
          VVHNormFactor = &
            &   ( wn(iwnsl)         * tanh( Planc*wn(iwnsl)/(2.0_DP*Boltz*Temp) ) )         &
            & / ( lparam%wnlc1(l,k) * tanh( Planc*lparam%wnlc1(l,k)/(2.0_DP*Boltz*Temp) ) )
          VVHNormFactor = 1.0_DP

          ! This cutoff factor is deduced from difference in line shapes
          ! between LBLRTM and my calculation.
          !
          CutOffFactor = &
            & - ( fredist( iwnsl ) / ( lparam % wnco1( l, k ) ) )**2 + 1.0d0
!!$          if ( wn( iwnsl ) < 1300.0d2 ) then
!!$            ! This scale of 4 cm-1 was found by trial and error. 2017/07/31.
!!$            CutOffFactor = &
!!$              & - ( fredist( iwnsl ) / ( 4.0d2 ) )**2 + 1.0d0
!!$          end if
          CutOffFactor = max( CutOffFactor, 0.0d0 )



          if ( FlagCO2SubLorentzianFactor ) then
            ! If this is applied, CutOffFactor is overwritten.

            if ( lparam % IMol == 2 ) then

              if ( lparam % pbwid1( l, k ) > 0.0d0 ) then

                !=========================================================
!!$                call lbl_CO2SubLorentianFactor( &
!!$                  & Temp, DV, &
!!$                  & FactorCO2CO2, FactorCO2N2 &
!!$                  & )
                !---------------------------------------------------------
                ! Unit of AbsDelWaveNum is cm-1
                AbsDelWaveNum = abs( DV ) * 1.0d-2
                do m = 1, SLPH89NMol
                  if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(1,m) ) then
                    m_SLPH89Factor(m) = 0.0d0
                  else if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(2,m) ) then
                    m_SLPH89Factor(m) = &
                      & - bm_SLPH89B(1,m) * ( AbsDelWaveNum - wm_SLPH89WaveNumBnds(1,m) )
                  else if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(3,m) ) then
                    m_SLPH89Factor(m) = &
                      & - bm_SLPH89B(1,m) * ( wm_SLPH89WaveNumBnds(2,m) - wm_SLPH89WaveNumBnds(1,m) ) &
                      & - bm_SLPH89B(2,m) * ( AbsDelWaveNum             - wm_SLPH89WaveNumBnds(2,m) )
                  else
                    m_SLPH89Factor(m) = &
                      & - bm_SLPH89B(1,m) * ( wm_SLPH89WaveNumBnds(2,m) - wm_SLPH89WaveNumBnds(1,m) ) &
                      & - bm_SLPH89B(2,m) * ( wm_SLPH89WaveNumBnds(3,m) - wm_SLPH89WaveNumBnds(2,m) ) &
                      & - bm_SLPH89B(3,m) * ( AbsDelWaveNum             - wm_SLPH89WaveNumBnds(3,m) )
                  end if
                end do
                m_SLPH89Factor = exp( m_SLPH89Factor )
                !
                FactorCO2CO2 = m_SLPH89Factor(1)
                FactorCO2N2  = m_SLPH89Factor(2)
                !=========================================================

                CutOffFactor = &
                  &   (   lparam % abwid1( l, k ) * FactorCO2N2    &
                  &     + lparam % sbwid1( l, k ) * FactorCO2CO2 ) &
                  &   / lparam % pbwid1( l, k )

              end if

            end if

          end if

          lsf = lsf * VVHNormFactor

          lsf = lsf * CutOffFactor


          lsf = lsf * lparam % renfac1( l, k )

          ! absorption coefficient per unit mass
          ! WARNINGWARNINGWARNING: This is wrong: need Avogadro constant.
!!$          ac( iwnsl ) = ac( iwnsl ) &
!!$            & + lparam % stren1( l, k ) * lsf / ( MolMassNum * amu )
          ! absorption coefficient per molecule
          ac( iwnsl ) = ac( iwnsl ) &
            & + lparam % stren1( l, k ) * lsf

        end do ! do l = ls, le
      end do


    end subroutine lbl_calc_ac

    !**************************************************************************


    subroutine lbl_CO2SubLorentianFactor( &
      & Temp, DV, &
      & FactorCO2CO2, FactorCO2N2 &
      & )

      real(dp), intent(in ) :: Temp
      real(dp), intent(in ) :: DV
      real(dp), intent(out) :: FactorCO2CO2
      real(dp), intent(out) :: FactorCO2N2


      ! Local variables
      !

      real(dp)           :: AbsDelWaveNum
      real(dp)           :: bm_SLPH89B    (SLPH89NB,SLPH89NMol)
      real(dp)           :: m_SLPH89Factor(SLPH89NMol)
      integer            :: m


      bm_SLPH89B = bm_SLPH89Alpha + bm_SLPH89Beta * exp( - bm_SLPH89Eps * Temp )


      ! Unit of AbsDelWaveNum is cm-1
      AbsDelWaveNum = abs( DV ) * 1.0d-2

      do m = 1, SLPH89NMol
        if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(1,m) ) then
          m_SLPH89Factor(m) = 0.0d0
        else if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(2,m) ) then
          m_SLPH89Factor(m) = &
            & - bm_SLPH89B(1,m) * ( AbsDelWaveNum - wm_SLPH89WaveNumBnds(1,m) )
        else if ( AbsDelWaveNum <= wm_SLPH89WaveNumBnds(3,m) ) then
          m_SLPH89Factor(m) = &
            & - bm_SLPH89B(1,m) * ( wm_SLPH89WaveNumBnds(2,m) - wm_SLPH89WaveNumBnds(1,m) ) &
            & - bm_SLPH89B(2,m) * ( AbsDelWaveNum             - wm_SLPH89WaveNumBnds(2,m) )
        else
          m_SLPH89Factor(m) = &
            & - bm_SLPH89B(1,m) * ( wm_SLPH89WaveNumBnds(2,m) - wm_SLPH89WaveNumBnds(1,m) ) &
            & - bm_SLPH89B(2,m) * ( wm_SLPH89WaveNumBnds(3,m) - wm_SLPH89WaveNumBnds(2,m) ) &
            & - bm_SLPH89B(3,m) * ( AbsDelWaveNum             - wm_SLPH89WaveNumBnds(3,m) )
        end if
      end do
      m_SLPH89Factor = exp( m_SLPH89Factor )



      FactorCO2CO2 = m_SLPH89Factor(1)
      FactorCO2N2  = m_SLPH89Factor(2)


    end subroutine lbl_CO2SubLorentianFactor

    !**************************************************************************
    ! function lbl_lineshape_voigt
    ! calculate line shape factor of voigt type
    !**************************************************************************
    ! lorwid        : lorentz halfwidth (m^-1)
    ! dopwid        : doppler halfwidth (m^-1)
    ! dv            : dv=(v-v0) : v  : wave number
    !               :           : v0 : wave number at line center
    !               : this value will be casted absolute value
    !**************************************************************************
    function lbl_voigt_h82( lorwid, dopwid, dv ) result( lsf )
      !************************************************************************
      !****
      !****  routine computes the voigt function:                     ***
      !****     y/pi*integral from - to + infinity of                 ***
      !                              exp(-t*t)/(y*y+(x-t)*(x-t)) dt   ***
      !****
      !****  Further, return value is divided by dopwid * sqrt(pi).
      !
      ! This method is described by Humlicek, JQSRT, 27, 437, 1982
      !

      real(dp), intent(in) :: lorwid, dopwid, dv

      real(dp)             :: lsf


      !
      ! local variables
      !
      real(dp)    :: x, y, s
      complex(dp) :: w4, t, u


      include 'lbl_voigt_h82.f90'


    end function lbl_voigt_h82

    !**************************************************************************

    subroutine LBLVoigtH82TblInit

      real(dp)             :: lsf


      !
      ! local variables
      !

      integer  :: NXPer10
      integer  :: NYPer10

      real(dp) :: MinDopWid
      real(dp) :: MaxLorWid
      real(dp) :: MaxDelWN
      real(dp) :: MaxX
      real(dp) :: MaxY

      real(dp) :: DopWid
      real(dp) :: LorWid
      real(dp) :: DelWN

      real(dp) :: LSFTrue
      real(dp) :: LSFTbl
      real(dp) :: LSFError

      integer :: i
      integer :: j


      !         dowid1( l ) = wnlc1( l ) / c &
      !              * sqrt( 2.0d0 * boltz * temp / mmass( iiso( l ) ) )
      MinDopWid = &
        & 1.0d0 / 3.0d8 &
        &   * sqrt( 2.0d0 * 1.0d-23 * 10.0d0 / ( 100.0d0 * 2.0d-27 ) )
      MaxLorWid = 100.0d0
      MaxDelWN = 25.0d2

      MaxX = MaxDelWN  / MinDopWid
      MaxY = MaxLorWid / MinDopWid

      NXPer10 = 200
      NLSFTblX = ( log10( MaxX ) - log10( MinDopWid / 100.0d0 ) ) * NXPer10 + 1
      NYPer10 = 200
      NLSFTblY = ( log10( MaxY ) - log10( MinDopWid / 100.0d0 ) ) * NYPer10 + 1

      write( 6, * ) 'NLSFTblX: ', NLSFTblX
      write( 6, * ) 'NLSFTblY: ', NLSFTblY

      allocate( a_LSFTblX( 0:NLSFTblX ) )
      allocate( a_LSFTblY( 0:NLSFTblY ) )
      allocate( aa_LSFTblLSF( 0:NLSFTblX, 0:NLSFTblY ) )

      DelLog10LSFTblX = &
        & ( log10( MaxX ) - log10( MinDopWid / 100.0d0 ) ) / NLSFTblX
      Log10LSFTblX1   = log10( MinDopWid / 100.0d0 )
      DelLog10LSFTblY = &
        & ( log10( MaxY ) - log10( MinDopWid / 100.0d0 ) ) / NLSFTblY
      Log10LSFTblY1   = log10( MinDopWid / 100.0d0 )

      i = 0
      a_LSFTblX(i) = 0.0d0
      do i = 1, NLSFTblX
        a_LSFTblX(i) = DelLog10LSFTblX * ( i - 1 ) + Log10LSFTblX1
        a_LSFTblX(i) = 10.0d0**a_LSFTblX(i)
      end do
      j = 0
      a_LSFTblY(j) = 0.0d0
      do j = 1, NLSFTblY
        a_LSFTblY(j) = DelLog10LSFTblY * ( j - 1 ) + Log10LSFTblY1
        a_LSFTblY(j) = 10.0d0**a_LSFTblY(j)
      end do


      do j = 1, NLSFTblY
        do i = 1, NLSFTblX
          DopWid = 1.0d0
          LorWid = a_LSFTblY(j) * DopWid
          DelWN  = a_LSFTblX(i) * DopWid
          aa_LSFTblLSF(i,j) = lbl_voigt_h82( LorWid, DopWid, DelWN )
        end do
      end do

      do j = 1+1, NLSFTblY
        do i = 1+1, NLSFTblX
          DopWid = 1.0d0
          LorWid = ( a_LSFTblY(j-1) + a_LSFTblY(j) ) / 2.0d0 * DopWid
          DelWN  = ( a_LSFTblX(i-1) + a_LSFTblX(i) ) / 2.0d0 * DopWid
          LSFTrue = lbl_voigt_h82( LorWid, DopWid, DelWN )
          LSFTbl  = LBLVoigtH82CalcWithTbl( LorWid, DopWid, DelWN )
          LSFError = ( LSFTrue - LSFTbl ) / LSFTrue
          if ( ( abs( LSFError ) > 1.0d-3 ) .and. ( LSFTrue > 1.0d-6 ) ) then
            write( 6, * ) "Error is too large: ", LSFTrue, LSFTbl, LSFError
!!$            write( 6, * ) "Error is too large: ", LorWid, DopWid, DelWN
          end if
        end do
      end do


    end subroutine LBLVoigtH82TblInit

    !**************************************************************************

    subroutine LBLVoigtH82TblFinalize


      !
      ! local variables
      !


      deallocate( a_LSFTblX )
      deallocate( a_LSFTblY )
      deallocate( aa_LSFTblLSF )


    end subroutine LBLVoigtH82TblFinalize

    !**************************************************************************

    function LBLVoigtH82CalcWithTbl( LorWid, DopWid, DV ) result( LSF )

      real(dp), intent(in) :: LorWid
      real(dp), intent(in) :: DopWid
      real(dp), intent(in) :: DV

      real(dp)             :: LSF


      !
      ! local variables
      !
      real(dp)    :: XVal
      real(dp)    :: YVal
      real(dp)    :: Log10X
      real(dp)    :: Log10Y

      integer, parameter :: NIntpol = 3
      real(dp)           :: Val
      real(dp)           :: a_Grid   (NIntpol)
      real(dp)           :: aa_Val   (NIntpol,NIntpol)
      real(dp)           :: a_Weight (NIntpol)
      real(dp)           :: a_WeightX(NIntpol)
      real(dp)           :: a_WeightY(NIntpol)

      integer     :: IndexX
      integer     :: IndexY
      integer     :: m
      integer     :: n
      integer     :: n1
      integer     :: n2
      integer     :: l1
      integer     :: l2


      XVal = abs( DV ) / DopWid
      YVal = LorWid    / DopWid

      ! check table size
      if ( XVal > a_LSFTblX(NLSFTblX) ) then
        write( 6, * ) 'DV/DopWid is out of range of table: ', &
          & XVal, a_LSFTblX(NLSFTblX)
        stop
      end if
      if ( YVal > a_LSFTblY(NLSFTblY) ) then
        write( 6, * ) 'LorWid/DopWid is out of range of table: ', &
          & YVal, a_LSFTblY(NLSFTblY)
        stop
      end if

      if ( XVal > 0.0d0 ) then
        Log10X = log10( XVal )
        if ( ( log10X - Log10LSFTblX1 ) > 0.0d0 ) then
          IndexX = int( ( log10X - Log10LSFTblX1 ) / DelLog10LSFTblX ) + 1
          IndexX = min( IndexX, NLSFTblX-NIntpol )
        else
          IndexX = 0
        end if
      else
        IndexX = 0
      end if

      if ( YVal > 0.0d0 ) then
        Log10Y = log10( YVal )
        if ( ( log10Y - Log10LSFTblY1 ) > 0.0d0 ) then
          IndexY = int( ( log10Y - Log10LSFTblY1 ) / DelLog10LSFTblY ) + 1
          IndexY = min( IndexY, NLSFTblY-NIntpol )
        else
          IndexY = 0
        end if
      else
        IndexY = 0
      end if


      Val = XVal
      do n = 1, NIntpol
        a_Grid(n) = a_LSFTblX(n+IndexX-1)
      end do
      !
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do m = 1, NIntpol
          if ( m == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(m) ) / ( a_Grid(n) - a_Grid(m) )
        end do
      end do
      !
      a_WeightX = a_Weight

      Val = YVal
      do n = 1, NIntpol
        a_Grid(n) = a_LSFTblY(n+IndexY-1)
      end do
      !
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do m = 1, NIntpol
          if ( m == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(m) ) / ( a_Grid(n) - a_Grid(m) )
        end do
      end do
      !
      a_WeightY = a_Weight


      do n2 = 1, NIntpol
        do n1 = 1, NIntpol
          l1 = n1+IndexX-1
          l2 = n2+IndexY-1
          aa_Val(n1,n2) = aa_LSFTblLSF(l1,l2)
        end do
      end do

      LSF = 0.0d0
      do n2 = 1, NIntpol
        do n1 = 1, NIntpol
          LSF = LSF + aa_Val(n1,n2) * a_WeightX(n1) * a_WeightY(n2)
        end do
      end do
      LSF = max( LSF, 0.0d0 )


    end function LBLVoigtH82CalcWithTbl

    !**************************************************************************

    function zominmax( val, minval, maxval ) result( zo )

      real(dp), intent(in) :: val, minval, maxval

      real(dp)             :: zo

      zo   = ( 1.0d0 + sign( 1.0d0, val - minval ) ) / 2.0d0 &
           * ( 1.0d0 - sign( 1.0d0, val - maxval ) ) / 2.0d0

    end function zominmax

    !**************************************************************************

    function zomin( val, minval ) result( zo )

      real(dp), intent(in) :: val, minval

      real(dp)             :: zo

      zo   = ( 1.0d0 + sign( 1.0d0, val - minval ) ) / 2.0d0

    end function zomin

    !**************************************************************************

    function zomax( val, maxval ) result( zo )

      real(dp), intent(in) :: val, maxval

      real(dp)             :: zo

      zo   = ( 1.0d0 - sign( 1.0d0, val - maxval ) ) / 2.0d0

    end function zomax

    !**************************************************************************
    ! function lbl_lineshape_voigt
    ! calculate line shape factor of voigt type
    !**************************************************************************
    ! lorwid        : lorentz halfwidth (m^-1)
    ! dopwid        : doppler halfwidth (m^-1)
    ! dv            : dv=(v-v0) : v  : wave number
    !               :           : v0 : wave number at line center
    !               : this value will be casted absolute value
    !**************************************************************************

    !**************************************************************************
    function lbl_voigt_fs85( lorwid, dopwid, dv ) result( voigt )
      !************************************************************************
      !****
      !****  routine computes the voigt function:                     ***
      !****     y/pi*integral from - to + infinity of                 ***
      !                              exp(-t*t)/(y*y+(x-t)*(x-t)) dt   ***
      !****
      !****  Return value seems to be divided by dopwid * sqrt(pi).
      !
      ! This is written in Liou, Radiation and Cloud Processes in the 
      ! Atmosphere, p31. 

      use constants0, only : pi

      real(dp), intent(in) :: lorwid, dopwid, dv

      real(dp)             :: voigt


      !
      ! local variables
      !
      real(dp)            :: ln2, c1, voiwid, xi, eta, etasq, pivoiwid


      ln2      = log( 2.0d0 )
      c1       = sqrt( lorwid * lorwid + 4.0d0 * ln2 * dopwid * dopwid )
      voiwid   = 0.5d0 * ( lorwid + c1 ) &
           + 0.05d0 * lorwid * ( 1.0d0 - ( lorwid + lorwid ) / ( lorwid + c1 ) )
      xi       = lorwid / voiwid
      eta      = dv / voiwid
      etasq    = eta * eta
      pivoiwid = pi * voiwid

      voigt = sqrt( ln2 / pi ) / voiwid &
           * ( 1.0d0 - xi ) * exp( -ln2 * etasq ) &
           + xi / pivoiwid &
           * ( 1.0d0 / ( 1.0d0 + etasq ) &
           - ( 1.0d0 - xi ) * ( 1.5d0 / ln2 + 1.0d0 + xi ) &
           * ( 0.066d0 * exp( -0.4d0 * etasq ) &
           - 1.0d0 / ( 40.0d0 - 5.5d0 * etasq + etasq * etasq ) ) )


    end function lbl_voigt_fs85

    !**************************************************************************
    ! function lbl_lineshape_voigt
    ! calculate line shape factor of voigt type
    !**************************************************************************
    ! lorwid        : lorentz halfwidth (m^-1)
    ! dopwid        : doppler halfwidth (m^-1)
    ! dv            : dv=(v-v0) : v  : wave number
    !               :           : v0 : wave number at line center
    !               : this value will be casted absolute value
    !**************************************************************************

!!$    !**************************************************************************
!!$    function lbl_lineshape_voigt_for_normalize( lorwid, dopwid, wnco, dwn ) &
!!$         result( renorm )
!!$      !************************************************************************
!!$      !****
!!$      !****  routine computes the voigt function:                     ***
!!$      !****     y/pi*integral from - to + infinity of                 ***
!!$      !                              exp(-t*t)/(y*y+(x-t)*(x-t)) dt   ***
!!$      !****
!!$      !****  Return value seems to be divided by dopwid * sqrt(pi).
!!$      !
!!$      ! This is written in Liou, Radiation and Cloud Processes in the 
!!$      ! Atmosphere, p31. 
!!$
!!$      use const_module, only : pi
!!$
!!$      real(dp), intent(in) :: lorwid, dopwid, wnco, dwn
!!$
!!$      real(dp)             :: renorm
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp)            :: ln2, c1, voiwid, xi, eta, etasq, pivoiwid, dv
!!$      real(dp)            :: fac1, fac2, fac3
!!$      real(dp)            :: voigt
!!$      real(dp)            :: l_dwn
!!$      integer(i4b)        :: l, lm
!!$
!!$
!!$      ln2      = log( 2.0d0 )
!!$      c1       = sqrt( lorwid * lorwid + 4.0d0 * ln2 * dopwid * dopwid )
!!$      voiwid   = 0.5d0 * ( lorwid + c1 ) &
!!$           + 0.05d0 * lorwid * ( 1.0d0 - ( lorwid + lorwid ) / ( lorwid + c1 ) )
!!$      xi       = lorwid / voiwid
!!$      pivoiwid = pi * voiwid
!!$
!!$
!!$      fac1 = sqrt( ln2 / pi ) / voiwid * ( 1.0d0 - xi )
!!$      fac2 = xi / pivoiwid
!!$      fac3 = ( 1.0d0 - xi ) * ( 1.5d0 / ln2 + 1.0d0 + xi )
!!$
!!$
!!$      lm     = int( wnco / dwn + 0.5 )
!!$      l_dwn  = wnco / dble( lm )
!!$
!!$
!!$      renorm = 0.0d0
!!$      do l = 1, lm
!!$         dv       = l_dwn * dble( l - 1 ) + l_dwn * 0.5d0
!!$
!!$         eta      = dv / voiwid
!!$         etasq    = eta * eta
!!$
!!$         voigt = fac1 * exp( -ln2 * etasq ) &
!!$              + fac2 &
!!$              * ( 1.0d0 / ( 1.0d0 + etasq ) &
!!$              - fac3 &
!!$              * ( 0.066d0 * exp( -0.4d0 * etasq ) &
!!$              - 1.0d0 / ( 40.0d0 - 5.5d0 * etasq + etasq * etasq ) ) )
!!$
!!$         renorm = renorm + voigt * l_dwn
!!$      end do
!!$
!!$
!!$    end function lbl_lineshape_voigt_for_normalize

!!$    !**************************************************************************
!!$    ! real*8 function voigt
!!$    ! calculate line shape factor of voigt type
!!$    !**************************************************************************
!!$    ! lorwid        : lorentz halfwidth (m^-1)
!!$    ! dopwid        : doppler halfwidth (m^-1)
!!$    ! dv            : dv=(v-v0) : v  : wave number
!!$    !               :           : v0 : wave number at line center
!!$    !               : this value will be casted absolute value
!!$    !**************************************************************************
!!$
!!$    !**************************************************************************
!!$    FUNCTION lbl_VOIGT_d76( lorwid, dopwid, dv ) result( voigt )
!!$      !************************************************************************
!!$      !****
!!$      !****  ROUTINE COMPUTES THE VOIGT FUNCTION: Y/PI*INTEGRAL FROM  ***
!!$      !****   - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT       ***
!!$      !****
!!$
!!$      !****  Return value seems to be divided by dopwid * sqrt(pi).
!!$
!!$      use const_module, only : pi
!!$
!!$      real(dp), intent(in) :: lorwid, dopwid, dv
!!$
!!$      real(dp)             :: voigt
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      REAL(dp) :: B(22),RI(15),XN(15),YN(15)
!!$      REAL(dp) :: D0(25),D1(25),D2(25),D3(25),D4(25)
!!$      REAL(dp) :: HN(25),H,XX(3),HH(3),NBY2(19),C(21)
!!$      LOGICAL  :: TRU
!!$      DATA B(1)/0.0/,B(2)/0.7093682E-7/
!!$      DATA XN/10.0,9.0,2*8.0,7.0,6.0,5.0,4.0,7*3.0/
!!$      DATA YN/3*0.6,0.5,2*0.4,4*0.3,1.0,0.9,0.8,2*0.7/
!!$      DATA H/0.201/
!!$      DATA XX/0.5246476,1.65068,0.7071068/
!!$      DATA HH/0.2562121,0.2588268E-1,0.2820948/
!!$      DATA NBY2/9.5,9.0,8.5,8.0,7.5,7.0,6.5,6.0,5.5,5.0,4.5,4.0,3.5,3.0, &
!!$           2.5,2.0,1.5,1.0,0.5/
!!$      DATA C / 0.7093682E-7,-0.2518434E-6, 0.8566874E-6,-0.2787638E-5, &
!!$               0.8668740E-5,-0.2565551E-4, 0.7228775E-4,-0.1933631E-3, &
!!$               0.4899520E-3,-0.1173267E-2, 0.2648762E-2,-0.5623190E-2, &
!!$               0.1119681E-1,-0.2084976E-1, 0.3621573E-1,-0.5851412E-1, &
!!$               0.8770816E-1,-0.121664, 0.15584, -0.184, 0.2           /
!!$      DATA TRU/.FALSE./
!!$
!!$      real(dp)     :: x, y
!!$      real(dp)     :: co, uu, vv, u, y2, d, dx, v
!!$      integer(i4b) :: i, j, max, min, n
!!$
!!$      x = abs( dv ) / dopwid
!!$      y = lorwid / dopwid
!!$
!!$      IF (TRU) GO TO 104
!!$      !****
!!$      !****  REGION I: COMPUTE DAWSON'S FUNCTION AT MESH POINTS. ***
!!$      !****
!!$      TRU=.TRUE.
!!$      DO I=1,15
!!$        RI(I)=-I/2.0
!!$     END DO
!!$      DO I=1,25
!!$        HN(I)=H*(I-0.5)
!!$        CO=4*HN(I)*HN(I)/25-2
!!$        DO J=2,21
!!$          B(J+1)=CO*B(J)-B(J-1)+C(J)
!!$       END DO
!!$        D0(I)=HN(I)*(B(22)-B(21))/5
!!$        D1(I)=1-2*HN(I)*D0(I)
!!$        D2(I)=(HN(I)*D1(I)+D0(I))/RI(2)
!!$        D3(I)=(HN(I)*D2(I)+D1(I))/RI(3)
!!$        D4(I)=(HN(I)*D3(I)+D2(I))/RI(4)
!!$     END DO
!!$  104 CONTINUE
!!$      IF (X.GE.5.0) GO TO 112
!!$      IF (Y.LE.1.0) GO TO 110
!!$      IF (X.GT.1.85*(3.6-Y)) GO TO 112
!!$      !****
!!$      !****  REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***
!!$      !****
!!$      IF (Y.LT.1.45) GO TO 107
!!$      I=Y+Y
!!$      GO TO 108
!!$  107 CONTINUE
!!$      I=11*Y
!!$  108 CONTINUE
!!$      J=X+X+1.85
!!$      MAX=XN(J)*YN(I)+0.46
!!$      MIN=MIN0(16,21-2*MAX)
!!$      !****
!!$      !****  EVALUATE CONTINUED FRACTION. ***
!!$      !****
!!$      UU=Y
!!$      VV=X
!!$      DO J=MIN,19
!!$        U=NBY2(J)/(UU*UU+VV*VV)
!!$        UU=Y+U*UU
!!$        VV=X-U*VV
!!$     END DO
!!$      VOIGT=UU/(UU*UU+VV*VV)/1.772454
!!$      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
!!$      RETURN
!!$      !****
!!$  110 CONTINUE
!!$      Y2=Y*Y
!!$      IF (X+Y.GE.5.0) GO TO 113
!!$      !****
!!$      !****  REGION I: COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***
!!$      !****
!!$      N=X/H
!!$      D=X-HN(N+1)
!!$      U=(((D4(N+1)*DX+D3(N+1))*DX+D2(N+1))*DX+D1(N+1))*DX+D0(N+1)
!!$      V=1-2*X*U
!!$      !****
!!$      !****  TAYLOR SERIES EXPANSION ABOUT Y=0.0. ***
!!$      !****
!!$      VV=EXP(Y2-X*X)*COS(2*X*Y)/1.128379-Y*V
!!$      UU=-Y
!!$      MAX=5+(12.5-X)*0.8*Y
!!$      DO I=2,MAX,2
!!$        U=(X*V+U)/RI(I)
!!$        V=(X*U+V)/RI(I+1)
!!$        UU=-UU*Y2
!!$        VV=VV+V*UU
!!$     END DO
!!$      VOIGT=1.128379*VV
!!$      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
!!$      RETURN
!!$      !****
!!$  112 CONTINUE
!!$      Y2=Y*Y
!!$      IF (Y.LT.11-0.6875*X) GO TO 113
!!$      !****
!!$      !****  REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***
!!$      !****
!!$      U=X-XX(3)
!!$      V=X+XX(3)
!!$      VOIGT=Y*(HH(3)/(Y2+U*U)+HH(3)/(Y2+V*V))
!!$      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
!!$      RETURN
!!$      !****
!!$      !****  REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***
!!$      !****
!!$  113 CONTINUE
!!$      U=X-XX(1)
!!$      V=X+XX(1)
!!$      UU=X-XX(2)
!!$      VV=X+XX(2)
!!$      VOIGT=Y*(HH(1)/(Y2+U*U)+HH(1)/(Y2+V*V) &
!!$           +HH(2)/(Y2+UU*UU)+HH(2)/(Y2+VV*VV))
!!$      voigt=1.0e0/(dopwid*sqrt(pi))*voigt
!!$      RETURN
!!$    END FUNCTION LBL_VOIGT_D76

    !**************************************************************************
    function lbl_lorentz( lorwid, dv ) result( lorentz )
      !************************************************************************
      !****
      !****  routine computes the Lorentz function: 
      !****           1/pi*lorwid/(lorwid**2+dv**2)
      !****

      use constants0, only : pi

      real(dp), intent(in) :: lorwid, dv

      real(dp)             :: lorentz


      !
      ! local variables
      !

      lorentz = lorwid / ( pi * ( lorwid * lorwid + dv * dv ) )


    end function lbl_lorentz

    !**************************************************************************
    ! real*8 function renorm
    ! Renormalize Line Shape
    !**************************************************************************

    function lbl_lineshape_renorm( lowid, dowid, wnmax, dwn ) result( renorm )


      real(dp), intent(in ) :: lowid, dowid, wnmax, dwn

      real(dp)              :: renorm


      !
      ! local variables
      !
      integer(i4b) :: l, lm
      real(dp)     :: l_dwn
      real(dp)     :: wn
      real(dp)     :: val1, val2


!!$      dwn    = dowid / 5.0d0
      lm     = int( wnmax / dwn + 0.5 ) + 1
      l_dwn  = wnmax / dble( lm - 1 )

      renorm = 0.0d0

      l = 1
      wn   = l_dwn * dble( l - 1 )
      val1 = lbl_voigt_fs85( lowid, dowid, wn )
      val2 = val1
      do l = 1+1, lm
         wn     = l_dwn * dble( l - 1 )
         val1   = lbl_voigt_fs85( lowid, dowid, wn )
         renorm = renorm + ( val1 + val2 ) * l_dwn * 0.5d0
         val2   = val1


!!$         val1   = lbl_lineshape_voigt( lowid, dowid, wn )
!!$         val2   = lbl_lineshape_voigt2( lowid, dowid, wn )
!!$         renorm = lbl_lineshape_lorentz( lowid, wn )
!!$         write( 6, * ) wn, val1, val2, renorm
!!$         if( l .eq. lm ) stop


      end do
      renorm = renorm + renorm

      renorm = 1.0d0 / renorm


    end function lbl_lineshape_renorm

    !**************************************************************************
    ! function lbl_cf
    !**************************************************************************
    ! chi-factor for carbon dioxide
    ! This routine is obtained from the subroutine chi_fn in src/oprop.f
    ! that is a part of AER LBLRTM v10.1. 
    !**************************************************************************

    function lbl_cf( wn ) result( chi )

      real(dp), intent(in) :: wn
      real(dp)             :: chi


      !
      ! local variables
      !
      real(dp), parameter :: asubl = 0.8d0, bsubl = 10.0d0
      real(sp)            :: x0, y0, f, y1, y2, z0, z1, z2, c6, c4, c2, fi


      ! SET UP CHI SUB-LORENTZIAN FORM FACTOR FOR CARBON DIOXIDE
      ! POLYNOMIAL MATCHED TO AN EXPONENTIAL AT X0 = 10 CM-1

      ! 0 - 25 cm-1 on 0.1 cm-1 grid


      X0 = 10.
      Y0 = asubl*EXP(-x0/bsubl)
      F  = 1./bsubl
      Y1 = -F*Y0
      Y2 = Y1*((BSUBL-1)/X0-F)
      Z0 = (Y0-1)/X0**2
      Z1 = Y1/(2*X0)
      Z2 = Y2/2.
      C6 = (Z0-Z1+(Z2-Z1)/4.)/X0**4
      C4 = (Z1-Z0)/X0**2-2.*X0**2*C6
      C2 = Z0-X0**2*C4-X0**4*C6

      fi = wn * 1.0d-2
      IF (FI.LT.X0) THEN
         CHI = 1.+C2*FI**2+C4*FI**4+C6*FI**6
      ELSE
         CHI = asubl*EXP(-fi/bsubl)
      ENDIF


    end function lbl_cf

    !**************************************************************************
!!$
!!$    subroutine chi_fn (chi)
!!$
!!$      dimension chi(0:250)
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp), parameter :: asubl = 0.8d0, bsubl = 10.0d0
!!$
!!$
!!$      ! SET UP CHI SUB-LORENTZIAN FORM FACTOR FOR CARBON DIOXIDE
!!$      ! POLYNOMIAL MATCHED TO AN EXPONENTIAL AT X0 = 10 CM-1
!!$
!!$      ! 0 - 25 cm-1 on 0.1 cm-1 grid
!!$
!!$
!!$      dvchi = 0.1d0
!!$
!!$
!!$      X0 = 10.
!!$      Y0 = asubl*EXP(-x0/bsubl)
!!$      F  = 1./bsubl
!!$      Y1 = -F*Y0
!!$      Y2 = Y1*((BSUBL-1)/X0-F)
!!$      Z0 = (Y0-1)/X0**2
!!$      Z1 = Y1/(2*X0)
!!$      Z2 = Y2/2.
!!$      C6 = (Z0-Z1+(Z2-Z1)/4.)/X0**4
!!$      C4 = (Z1-Z0)/X0**2-2.*X0**2*C6
!!$      C2 = Z0-X0**2*C4-X0**4*C6
!!$
!!$      do ISUBL = 0, 250
!!$         FI = DVCHI* REAL(ISUBL)
!!$         IF (FI.LT.X0) THEN
!!$            CHI(ISUBL) = 1.+C2*FI**2+C4*FI**4+C6*FI**6
!!$         ELSE
!!$            CHI(ISUBL) = asubl*EXP(-fi/bsubl)
!!$         ENDIF
!!$      end do
!!$
!!$
!!$    end subroutine chi_fn
!!$
    !**************************************************************************

 end module lbl_module
