module ac_calc use vtype_module implicit none private public :: ACCalc !---------------------------------------------------------------------------- contains !---------------------------------------------------------------------------- subroutine ACCalc( & & FlagCalcLine, FlagCalcCont, & & DelWaveNum, WaveNumS, WaveNumE, LineCutOffWaveNum, & & HITRANDir, HITRANFileNameExt, & & OutputFileName, & & NLev, NMolNum, z_Press, m_MolNum, z_Temp, zm_VMR & & ) use lbl_const_module use fi_module use lbl_module use ac_io_module use mtckd_wrapper logical , intent(in ) :: FlagCalcLine logical , intent(in ) :: FlagCalcCont real(dp) , intent(in ) :: DelWaveNum real(dp) , intent(in ) :: WaveNumS real(dp) , intent(in ) :: WaveNumE real(dp) , intent(in ) :: LineCutOffWaveNum character(*), intent(in ) :: HITRANDir character(*), intent(in ) :: HITRANFileNameExt character(*), intent(in ) :: OutputFileName integer , intent(in ) :: NLev integer , intent(in ) :: NMolNum real(dp) , intent(in ) :: z_Press( NLev ) integer , intent(in ) :: m_MolNum( NMolNum ) real(dp) , intent(in ) :: z_Temp( NLev ) real(dp) , intent(in ) :: zm_VMR( NLev, NMolNum ) !!$ allocate( z_Press ( NLev ) ) !!$ allocate( m_MolNum( NMolNum ) ) !!$ allocate( z_Temp ( NLev ) ) !!$ allocate( zm_VMR ( NLev, NMolNum ) ) ! Local variables ! type( LineParam ) :: LParam integer :: NLine ! ! NWaveNumSubLoop : number of wavenumber sub loop ! : (inner most loop for optimization) ! integer :: NWaveNumSubLoop integer :: iwnsl integer :: iwnsls integer :: iwnsle real(dp) , allocatable :: w_WaveNum( : ) integer(i4b) :: NWaveNum character(len=extstr) :: FileName real(dp) :: DelWaveNumMod ! ! ach( NWaveNumSubLoop ) :: absorption coefficient ! real(dp) , allocatable :: w_AbsCoef(:) real(dp) , allocatable :: w_AbsCoefLine(:) real(dp) , allocatable :: w_AbsCoefCont(:) real(dp) :: VMRH2O real(dp) :: VMRCO2 real(dp) :: VMRO3 real(dp) :: VMRO2 real(dp) :: VMRN2 real(dp) :: VMROTHERS real(dp) :: NumDen real(dp) :: PathLength real(dp) :: ColDen integer :: iWaveNum integer :: k integer :: m NWaveNumSubLoop = 10 !!$ NWaveNumSubLoop = 4096 ! for pi-computer allocate( w_WaveNum ( NWaveNumSubLoop ) ) allocate( w_AbsCoef( NWaveNumSubLoop ) ) NWaveNum = int( ( WaveNumE - WaveNumS ) / DelWaveNum ) DelWaveNumMod = ( WaveNumE - WaveNumS ) / dble( NWaveNum ) allocate( w_AbsCoefLine( NWaveNumSubLoop ) ) allocate( w_AbsCoefCont( NWaveNumSubLoop ) ) call ac_io_output_init( & & OutputFileName, NMolNum, m_MolNum, NLev, z_Press, zm_VMR, z_Temp, NWaveNum & & ) do m = 1, NMolNum if ( FlagCalcLine ) then write( FileName, '(a,i2.2,a)' ) & & trim( HITRANDir ) // "/", m_MolNum( m ), trim( HITRANFileNameExt ) write( 6, '(a,i2)' ) 'Molecular number: ', m write( 6, '(a)' ) ' File: ' // trim( FileName ) NLine = fi_lc( FileName ) write( 6, '(a,i6)' ) ' Number of lines : ', NLine !!$ call lbl_lineparam_init( m_MolNum( m ), NLine, NLev, LParam ) call lbl_lineparam_init( m_MolNum( m ), NLine, 1, LParam ) !!$ call lbl_readfile ( FileName, LParam ) call lbl_readfile_hitran04( FileName, LParam ) ! ! wavenumber width over integration (m^-1) ! !!$ LParam % wnco( : ) = 3.0d2 !!$ LParam % wnco( : ) = 25.0d2 LParam % wnco( : ) = LineCutOffWaveNum end if do k = 1, NLev if( mod( k, NLev/10 ) == 0 ) then write( 6, '(2i6,a,i5,a)' ) k, NLev, & & " --- ", int( real( k ) / NLev * 100 ), "%" end if if ( FlagCalcLine ) then ! ! pressure and temperature correction of line parameters ! call lbl_lineparam_correct( & z_Press(k), z_Press(k) * zm_VMR(k,m), z_Temp(k), & DelWaveNumMod, & LParam % IMol , & LParam % NLine , & LParam % iiso , & LParam % wnlc0 , & LParam % stren0 , & LParam % lowene , & LParam % abwid0 , & LParam % sbwid0 , & LParam % tdep_abw , & LParam % tdep_sbw , & LParam % apshift , & LParam % wnco , & !!$ LParam % wnlc1 (:,k), & !!$ LParam % stren1 (:,k), & !!$ LParam % pbwid1 (:,k), & !!$ LParam % dowid1 (:,k), & !!$ LParam % renfac (:,k), & LParam % wnlc1 (:,1), & LParam % stren1 (:,1), & LParam % pbwid1 (:,1), & LParam % dowid1 (:,1), & LParam % renfac (:,1), & 1, LParam % NLine & & ) end if ! ! loop for wavenumber ! do iWaveNum = 1, NWaveNum, NWaveNumSubLoop !!$ if( mod( iWaveNum, NWaveNum/NWaveNumSubLoop/10 ) == 0 ) then !!$ write( 6, '(2i6,a,i5,a)' ) iWaveNum-1, NWaveNum, & !!$ & " --- ", int( real( iWaveNum-1 ) / NWaveNum * 100 ), "%" !!$ end if iwnsls = 1 iwnsle = min( NWaveNum - ( iWaveNum - 1 ), NWaveNumSubLoop ) do iwnsl = iwnsls, iwnsle w_WaveNum( iwnsl ) = & & WaveNumS + DelWaveNumMod * dble( iWaveNum - 1 + iwnsl - 1 ) end do if ( FlagCalcLine ) then call lbl_calc_ac( & & LParam, & & NWaveNumSubLoop, iwnsls, iwnsle, w_WaveNum, w_AbsCoefLine & & ) else w_AbsCoefLine = 0.0d0 end if !!$ call ac_io_output( m, k, NWaveNumSubLoop, w_WaveNum, w_AbsCoefLine, iwnsls, iwnsle, iWaveNum ) if ( FlagCalcCont ) then PathLength = 1.0d0 NumDen = z_Press(k) / ( Boltz * z_Temp(k) ) !!$ NumDen = NumDen * ( 1.0d0 - zm_VMR(k,1) ) * zm_VMR(k,m) NumDen = NumDen * zm_VMR(k,m) ColDen = NumDen * PathLength if ( NumDen > 0.0d0 ) then VMRH2O = zm_VMR(k,1) if ( NMolNum >= 2 ) then VMRCO2 = zm_VMR(k,2) else VMRCO2 = 0.0d0 end if if ( NMolNum >= 3 ) then VMRO3 = zm_VMR(k,3) else VMRO3 = 0.0d0 end if if ( NMolNum >= 7 ) then VMRO2 = zm_VMR(k,7) else VMRO2 = 0.0d0 end if if ( NMolNum >= 22 ) then VMRN2 = zm_VMR(k,22) else VMRN2 = 0.0d0 end if VMROTHERS = 1.0d0 - VMRH2O - VMRCO2 - VMRO3 - VMRO2 - VMRN2 call call_mtckd_wrapper( & & z_Press(k), z_Temp(k), VMRH2O, PathLength, & & m, & & VMRCO2, VMRO3, VMRO2, VMRN2, VMROTHERS, & !!$ & WaveNumS, DelWaveNumMod, & & WaveNumS + DelWaveNumMod * dble( iWaveNum - 1 + iwnsls - 1 ), DelWaveNumMod, & & NWaveNumSubLoop, & & w_AbsCoefCont & & ) ! multiplying a factor to convert definition of absorption ! coefficient (see p.6 of Ptashnik et al. (2011). ! Reference ! Ptashnik et al., 2011, JGR, 116, D16305, doi:10.1029/2011JD015603 w_AbsCoefCont = w_AbsCoefCont * ( 296.0d0 / z_Temp(k) ) ! dividing by column density of absorbing molecule w_AbsCoefCont = w_AbsCoefCont / ColDen else w_AbsCoefCont = 0.0d0 end if else w_AbsCoefCont = 0.0d0 end if !!$ call ac_io_output( m, k, NWaveNumSubLoop, & !!$ & "AbsCoefLine", w_WaveNum, w_AbsCoefLine, & !!$ & iwnsls, iwnsle, iWaveNum ) !!$ !!$ call ac_io_output( m, k, NWaveNumSubLoop, & !!$ & "AbsCoefCont", w_WaveNum, w_AbsCoefCont, & !!$ & iwnsls, iwnsle, iWaveNum ) call ac_io_output( m, k, NWaveNumSubLoop, & & "AbsCoef", w_WaveNum, w_AbsCoefLine+w_AbsCoefCont, & & iwnsls, iwnsle, iWaveNum ) end do end do if ( FlagCalcLine ) then call lbl_lineparam_finalize( LParam ) end if end do call ac_io_output_end deallocate( w_AbsCoefLine ) deallocate( w_AbsCoefCont ) end subroutine ACCalc !---------------------------------------------------------------------------- end module ac_calc