SUBROUTINE AbsorpCoeffiLBL(P, T, wavenumberMin, wavenumberMax, delta_wavenumber, HITRAN, kmax, SpLoop, path, absorb_number )
   use initialset
   implicit none

   INTRINSIC EXP
   integer,intent(in) :: kmax, SpLoop
   REAL(8), intent(in) :: wavenumberMin, wavenumberMax, delta_wavenumber  
   REAL(8), intent(in) :: P, T, HITRAN(linenumber(SpLoop),6)
   REAL(8), intent(in) :: path ! = number density * height * 5/3 
   REAL(8), intent(out) :: absorb_number(1:kmax)

   REAL(8) :: S, B, Q, SE ,voigt, x, y, hcok, Qiso(10) 
   REAL(8) :: gamma_D, gamma_L, wavenumber 
   integer :: i, k, l, IDmol, iso

   hcok = h* light_speed / bolztman

   !.... line parameter convert

   if ( SpLoop == 1 ) then
      IDmol = 2
   elseif ( SpLoop == 2 ) then
      IDmol = 1
   elseif ( SpLoop == 3 ) then
      IDmol = 4
   else
      write(*,'(a)') 'unknown molecude ID'  ; stop
   endif

   do iso = 1, int( maxval( HITRAN(:,6) ) )
      call TIPSglh (T, IDmol, iso, Qiso(iso) ) 
   enddo

   absorb_number(1:kmax) = 0.0

   l = 1

   loop_wavenumber : do k = 1, kmax

      wavenumber = wavenumberMin + delta_wavenumber * ( k - 0.5)
      if (wavenumber > wavenumberMax ) exit

      loop_line : do i = l, linenumber(SpLoop)

         if (  HITRAN(i,1) <  wavenumber - wing(SpLoop) ) then 
            l = l + 1

         elseif (  HITRAN(i,1) >  wavenumber + wing(SpLoop) ) then 
            exit

         elseif ( abs( HITRAN(i,1) - wavenumber ) <= wing(SpLoop) ) then 

            gamma_D = HITRAN(i,1) / light_speed * ( 2.0 * R_const * T / m(SpLoop) )**(0.5) 
            gamma_L = HITRAN(i,3) * P/P0 * (T0/T)**( HITRAN(i,5) )

            B = exp( - hcok * HITRAN(i,4)*1.0E+2 * ( 1.0E+0 / T - 1.0E+0/ T0 ) )
            Q = Qiso( int( hitran(i,6) ) )
            SE = (1.0 - exp ( - hcok * HITRAN(i,1) * 1.0E+2 / T) ) &
                 &   / (  1.0 - exp ( - hcok * HITRAN(i,1) * 1.0E+2 / T0 ) )

            S = HITRAN(i,2) * B * Q * SE
            x = ( wavenumber - HITRAN(i,1))/gamma_D
            y = gamma_L/gamma_D

            absorb_number(k) = absorb_number(k) + S * voigt(x,y) / gamma_D

         end if

      enddo loop_line
      
      absorb_number(k) = absorb_number(k) / pi**(0.5) * 1.0E-4 * path    ! 1.0E-4 [m^2] = 1 [cm^2]

      write(*,*) wavenumber, absorb_number(k) 
      
   enddo loop_wavenumber

END SUBROUTINE
