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_lineparam_correct 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) ! wnco(nline) : cutoff wave number (m^-1) ! 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( :, : ) real(dp) , pointer :: wnco( : ), renfac( :, : ) end type LineParam 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 !************************************************************************** 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 ) & ) allocate( & lp % wnco ( nline ) , & lp % renfac( 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 & ) deallocate( & lp % wnco, & lp % renfac & ) 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 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 lp % lowene ( l ) = lp % lowene ( l ) * 1.0d2 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 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 ), & 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 ! 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, * ) & "For ", l, "th line's self-broadened half width is zero" !!$ write( *, * ) l,"th line's self-broadened half width is zero" !!$ ! 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 lp % lowene ( l ) = lp % lowene ( l ) * 1.0d2 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 ) end subroutine lbl_readfile_hitran04 !************************************************************************** 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, & wnlc1, stren1, pbwid1, dowid1, renfac, & ls, le ) use lbl_const_module 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 ), wnco( linen ) real(dp) , intent(out) :: & wnlc1( linen ), stren1( linen ), pbwid1( linen ), & dowid1( linen ), renfac( 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 ) & & ) end do !!$ tips0( : ) = 1.0d0 !!$ tips ( : ) = 1.0d0 allocate( mmass( niso ) ) do i = 1, niso mmass( i ) = lbl_iso_mw( IMol, i ) * amu 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 end do do l = ls, le !!$ wnlc1 ( l ) = wnlc0( l ) wnlc1 ( l ) = wnlc0( l ) + apshift( l ) * ( tpress / press0 ) 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 ) ) pbwid1( l ) & = abwid0(l)*((tpress-ppress)/press0)*(temp0/temp)**tdep_abw(l) & + sbwid0(l)*(ppress /press0)*(temp0/temp)**tdep_sbw(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 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 renfac( : ) = 1.0d0 end subroutine lbl_lineparam_correct !************************************************************************** function 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 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( & & lparam, & & nwnsl, iwnsls, iwnsle, wn, ac & & ) !!$ use const_module, only : & !!$ & pi, & !!$ & amu use lbl_const_module, only : & & PI type( LineParam ), intent(in ) :: lparam integer(i4b) , intent(in ) :: nwnsl, iwnsls, iwnsle real(dp) , intent(in ) :: wn( nwnsl ) real(dp) , intent(out) :: ac( nwnsl ) ! ! local variables ! !!$ real(dp) :: MolMassNum 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 ! ! varialbes for inline expansion ! real(dp) :: x, y, s complex(dp) :: w4, t, u !!$ ! Set molecular mass number !!$ ! !!$ call lbl_molmn( & !!$ & lparam % IMol, & !!$ & MolMassNum & !!$ & ) 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_edge2 = lparam % nline l_mp = ( l_edge1 + l_edge2 ) / 2 do l = l_mp if( lparam%wnlc1(l,k) >= wn(iwnsl) - lparam%wnco(l) ) 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 % wnco( l ) ) ) & 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 % wnco( l ) ) ) & 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 = lsf * lparam % renfac( 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 !************************************************************************** !************************************************************************** ! 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 ! use lbl_const_module, only : pi 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 lbl_const_module, 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 lbl_const_module, 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 ) use lbl_const_module 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