module opt_prop_kd use vtype_module implicit none private public :: OptPropKDInit public :: OptPropKDFinalize public :: OptPropKDGetNWN public :: OptPropKDGetNPress public :: OptPropKDGetNMol public :: OptPropKDGetNBand public :: OptPropKDGetDelWN public :: OptPropKDGetBandNum public :: OptPropKDGetBandWNBnds public :: OptPropKDGetBandBinIndexBnds public :: OptPropKDGetStrBandAvePF !!$ public :: OptPropKDGetTotalStrFlux public :: OptPropKDGetAbsCoefProf public :: OptPropKDGetRayScatCoef public :: OptPropKDGetPtclParam public :: OptPropKDGetPFIntedRatio public :: OptPropKDGetStrPFIntedRatio character(128) , save :: OptPropNcFNSave integer , parameter :: NBnd = 2 real(DP), allocatable, save :: ab_KDTblBandWaveNumBnds (:,:) integer , allocatable, save :: ab_KDTblBandBinIndexBnds (:,:) ! real(DP), allocatable, save :: c_DelWaveNum(:) ! integer , save :: KDTblIndexVMRMol1Num integer , save :: KDTblIndexVMRMol2Num integer , save :: KDTblIndexVMRMol3Num integer , save :: KDTblNBin integer , save :: KDTblNBand integer , save :: KDTblNPress integer , save :: KDTblNMol integer , save :: KDTblNTemp integer , save :: KDTblNVMRMol1 integer , save :: KDTblNVMRMol2 integer , save :: KDTblNVMRMol3 integer , allocatable, save :: c_KDTblBinNum (:) real(DP), allocatable, save :: p_KDTblPress (:) integer , allocatable, save :: m_KDTblMolNum (:) real(DP), allocatable, save :: t_KDTblTemp (:) real(DP), allocatable, save :: u_KDTblVMRMol1(:) real(DP), allocatable, save :: v_KDTblVMRMol2(:) real(DP), allocatable, save :: w_KDTblVMRMol3(:) real(DP), allocatable, save :: ct_AtmPFRatio (:,:) real(DP), allocatable, save :: c_StrPFRatio (:) !!$ real(DP), allocatable, save :: b_BandAveStrPF (:) real(DP), allocatable, save :: b_RatioBandIntFluxtoTotFlux(:) real(DP), allocatable, save :: pmctuvw_AbsCoef(:,:,:,:,:,:,:) real(DP), allocatable, save :: b_BandAveRayScatCoefNonRadAct(:) real(DP), allocatable, save :: bm_BandAveRayScatCoef(:,:) !!$ integer :: NChar !!$ integer , allocatable :: a_CharNum(:) integer :: KDTblNPtcl integer :: KDTblPtclNRadiusMax integer , allocatable :: a_KDTblPtclNum (:) integer , allocatable :: r_KDTblPtclRadiusNum(:) character(stdstr), allocatable :: a_KDTblPtclName (:) integer , allocatable :: a_KDTblPtclNRadius (:) real(DP) , allocatable :: ra_KDTblPtclRadius (:,:) real(DP) , allocatable :: bra_KDTblBandAveQExt(:,:,:) real(DP) , allocatable :: bra_KDTblBandAveSSA (:,:,:) real(DP) , allocatable :: bra_KDTblBandAveAF (:,:,:) integer , allocatable, save :: mc_IDAbsCoefType(:,:) character(128), save :: ModuleName = 'opt_prop_kd' logical, save :: FlagInited = .false. contains !---------------------------------------------------------------------------- subroutine OptPropKDInit( & & OptPropNcFN & & ) use ni3_module character(*), intent(in) :: OptPropNcFN ! local variables ! integer :: InNcID character(stdstr) :: Mode character(128) :: PressAxisName character(128) :: MolAxisName character(128) :: TempAxisName character(128) :: VMRMol1AxisName character(128) :: VMRMol2AxisName character(128) :: VMRMol3AxisName character(128) :: BinAxisName character(128) :: AbsCoefVarName character(128) :: AtmPFRatioVarName character(128) :: StrPFRatioVarName real(DP), allocatable :: cp_StrPFRatio(:,:) !!$ real(DP), allocatable :: bp_BandAveStrPF(:,:) real(DP), allocatable :: bp_RatioBandIntFluxtoTotFlux(:,:) integer :: VMRMol1Num integer :: VMRMol2Num integer :: VMRMol3Num integer :: iBin integer :: iBand integer :: m !!$ if ( .not. FlagInited ) then !!$ write( 6, * ) trim( ModuleName ), " is not initialized." !!$ stop !!$ end if KDTblNBin = -1 KDTblNBand = -1 KDTblNPress = -1 KDTblNMol = -1 KDTblNTemp = -1 KDTblNVMRMol1 = -1 KDTblNVMRMol2 = -1 KDTblNVMRMol3 = -1 OptPropNcFNSave = OptPropNcFN PressAxisName = "Press" MolAxisName = "MolNum" TempAxisName = "Temp" VMRMol1AxisName = "VMRMol1" VMRMol2AxisName = "VMRMol2" VMRMol3AxisName = "VMRMol3" BinAxisName = 'BinNum' AbsCoefVarName = 'BinAbsCoef' AtmPFRatioVarName = 'BinAtmPFRatio' StrPFRatioVarName = 'BinStrPFRatio' Mode = "read" call ni3_open( OptPropNcFNSave, Mode, InNcID ) call ni3_inq_dimlen( InNcID, BinAxisName , KDTblNBin ) call ni3_inq_dimlen( InNcID, 'BandNum' , KDTblNBand ) call ni3_inq_dimlen( InNcID, PressAxisName , KDTblNPress ) call ni3_inq_dimlen( InNcID, MolAxisName , KDTblNMol ) call ni3_inq_dimlen( InNcID, TempAxisName , KDTblNTemp ) call ni3_inq_dimlen( InNcID, VMRMol1AxisName , KDTblNVMRMol1 ) call ni3_inq_dimlen( InNcID, VMRMol2AxisName , KDTblNVMRMol2 ) call ni3_inq_dimlen( InNcID, VMRMol3AxisName , KDTblNVMRMol3 ) allocate( ab_KDTblBandWaveNumBnds (NBnd,KDTblNBand) ) allocate( ab_KDTblBandBinIndexBnds (NBnd,KDTblNBand) ) call ni3_get_var( InNcID, 'BandWaveNumBnds' , ab_KDTblBandWaveNumBnds ) call ni3_get_var( InNcID, 'BandBinIndexBnds', ab_KDTblBandBinIndexBnds ) allocate( c_DelWaveNum(KDTblNBin) ) call ni3_get_var( InNcID, "BinWeight", c_DelWaveNum ) ! call ni3_get_var( InNcID, "VMRMol1Num", VMRMol1Num ) call ni3_get_var( InNcID, "VMRMol2Num", VMRMol2Num ) call ni3_get_var( InNcID, "VMRMol3Num", VMRMol3Num ) ! allocate( c_KDTblBinNum (KDTblNBin ) ) allocate( p_KDTblPress (KDTblNPress ) ) allocate( m_KDTblMolNum (KDTblNMol ) ) allocate( t_KDTblTemp (KDTblNTemp ) ) allocate( u_KDTblVMRMol1(KDTblNVMRMol1) ) allocate( v_KDTblVMRMol2(KDTblNVMRMol2) ) allocate( w_KDTblVMRMol3(KDTblNVMRMol3) ) ! call ni3_get_var( InNcID, BinAxisName , c_KDTblBinNum ) call ni3_get_var( InNcID, PressAxisName , p_KDTblPress ) call ni3_get_var( InNcID, MolAxisName , m_KDTblMolNum ) call ni3_get_var( InNcID, TempAxisName , t_KDTblTemp ) call ni3_get_var( InNcID, VMRMol1AxisName, u_KDTblVMRMol1 ) call ni3_get_var( InNcID, VMRMol2AxisName, v_KDTblVMRMol2 ) call ni3_get_var( InNcID, VMRMol3AxisName, w_KDTblVMRMol3 ) ! allocate( pmctuvw_AbsCoef(KDTblNPress,KDTblNMol,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) ) call ni3_get_var( InNcID, AbsCoefVarName, pmctuvw_AbsCoef ) ! allocate( mc_IDAbsCoefType(KDTblNMol,KDTblNBin) ) call ni3_get_var( InNcID, "BinIDAbsCoefType", mc_IDAbsCoefType ) ! allocate( ct_AtmPFRatio(KDTblNBin,KDTblNTemp) ) call ni3_get_var( InNcID, AtmPFRatioVarName, ct_AtmPFRatio ) ! allocate( cp_StrPFRatio (KDTblNBin,KDTblNPress) ) call ni3_get_var( InNcID, StrPFRatioVarName, cp_StrPFRatio ) allocate( c_StrPFRatio (KDTblNBin) ) do iBin = 1, KDTblNBin c_StrPFRatio(iBin) = cp_StrPFRatio(iBin,1) end do deallocate( cp_StrPFRatio ) ! !!$ allocate( bp_BandAveStrPF(KDTblNBand,KDTblNPress) ) !!$ call ni3_get_var( InNcID, StrPFRatioVarName, bp_BandAveStrPF ) !!$ allocate( b_BandAveStrPF(KDTblNBand) ) !!$ do iBand = 1, KDTblNBand !!$ b_BandAveStrPF(iBand) = bp_BandAveStrPF(iBand,1) !!$ end do !!$ deallocate( bp_BandAveStrPF ) ! allocate( b_BandAveRayScatCoefNonRadAct(KDTblNBand) ) allocate( bm_BandAveRayScatCoef (KDTblNBand,KDTblNMol) ) call ni3_get_var( InNcID, "BandAveRayScatCoefNonRadAct", b_BandAveRayScatCoefNonRadAct ) call ni3_get_var( InNcID, "BandAveRayScatCoef" , bm_BandAveRayScatCoef ) !!$ call ni3_inq_dimlen( InNcID, "CharNum" , NChar ) call ni3_inq_dimlen( InNcID, "PtclNum" , KDTblNPtcl ) call ni3_inq_dimlen( InNcID, "PtclRadiusNum", KDTblPtclNRadiusMax ) !!$ allocate( a_CharNum (NChar ) ) allocate( a_KDTblPtclNum (KDTblNPtcl ) ) allocate( r_KDTblPtclRadiusNum(KDTblPtclNRadiusMax) ) allocate( a_KDTblPtclName (KDTblNPtcl ) ) allocate( a_KDTblPtclNRadius (KDTblNPtcl ) ) allocate( ra_KDTblPtclRadius (KDTblPtclNRadiusMax,KDTblNPtcl) ) !!$ call ni3_get_var( InNcID, "CharNum" , a_CharNum ) call ni3_get_var( InNcID, "PtclNum" , a_KDTblPtclNum ) call ni3_get_var( InNcID, "PtclRadiusNum", r_KDTblPtclRadiusNum ) call ni3_get_var( InNcID, "PtclName" , a_KDTblPtclName ) call ni3_get_var( InNcID, "PtclNRadius" , a_KDTblPtclNRadius ) call ni3_get_var( InNcID, "PtclRadius" , ra_KDTblPtclRadius ) allocate( bra_KDTblBandAveQExt (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) ) allocate( bra_KDTblBandAveSSA (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) ) allocate( bra_KDTblBandAveAF (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) ) call ni3_get_var( InNcID, "BandAveQExt" , bra_KDTblBandAveQExt ) call ni3_get_var( InNcID, "BandAveSSA" , bra_KDTblBandAveSSA ) call ni3_get_var( InNcID, "BandAveAF" , bra_KDTblBandAveAF ) allocate( bp_RatioBandIntFluxtoTotFlux(KDTblNBand,KDTblNPress) ) !!$ call ni3_get_var( InNcID, "RatioBandIntFluxtoTotFlux", pb_RatioBandIntFluxtoTotFlux ) !!$ allocate( b_RatioBandIntFluxtoTotFlux(KDTblNBand) ) !!$ do iBand = 1, KDTblNBand !!$ b_RatioBandIntFluxtoTotFlux(iBand) = bp_RatioBandIntFluxtoTotFlux(iBand,1) !!$ end do !!$ deallocate( bp_RatioBandIntFluxtoTotFlux ) !!$ call ni3_close( InNcID ) !!$ Mode = "read" !!$ call ni3_open( StrSpeNcFN, Mode, InNcID ) !!$ call ni3_inq_dimlen( InNcID, "KDTblNBand", TmpNum ) !!$ if ( TmpNum /= KDTblNBand ) then !!$ write( 6, * ) 'Array size for stellar spectrum is inappropriate: ', & !!$ TmpNum, KDTblNBand !!$ stop !!$ end if call ni3_get_var( InNcID, "RatioBandIntFluxtoTotFlux", bp_RatioBandIntFluxtoTotFlux ) allocate( b_RatioBandIntFluxtoTotFlux(KDTblNBand) ) do iBand = 1, KDTblNBand b_RatioBandIntFluxtoTotFlux(iBand) = bp_RatioBandIntFluxtoTotFlux(iBand,1) end do deallocate( bp_RatioBandIntFluxtoTotFlux ) call ni3_close( InNcID ) do m = 1, KDTblNMol if ( m_KDTblMolNum(m) == VMRMol1Num ) exit end do if ( m > KDTblNMol ) then stop 'Unexpected events when searching index for VMRMol1Num.' end if KDTblIndexVMRMol1Num = m ! do m = 1, KDTblNMol if ( m_KDTblMolNum(m) == VMRMol2Num ) exit end do if ( m > KDTblNMol ) then stop 'Unexpected events when searching index for VMRMol2Num.' end if KDTblIndexVMRMol2Num = m ! do m = 1, KDTblNMol if ( m_KDTblMolNum(m) == VMRMol3Num ) exit end do if ( m > KDTblNMol ) then stop 'Unexpected events when searching index for VMRMol3Num.' end if KDTblIndexVMRMol3Num = m ! check key and non-key molecules if ( mc_IDAbsCoefType(KDTblIndexVMRMol1Num,1) == 0 ) then if ( KDTblNVMRMol1 > 1 ) then write( 6, * ) & & 'VMR Axis 1 has an unexpected size nevertheless it is non-key molecule: ', & & KDTblNVMRMol1 stop end if end if if ( mc_IDAbsCoefType(KDTblIndexVMRMol2Num,2) == 0 ) then if ( KDTblNVMRMol2 > 1 ) then write( 6, * ) & & 'VMR Axis 2 has an unexpected size nevertheless it is non-key molecule: ', & & KDTblNVMRMol2 stop end if end if if ( mc_IDAbsCoefType(KDTblIndexVMRMol3Num,1) == 0 ) then if ( KDTblNVMRMol3 > 1 ) then write( 6, * ) & & 'VMR Axis 3 has an unexpected size nevertheless it is non-key molecule: ', & & KDTblNVMRMol3 stop end if end if FlagInited = .true. end subroutine OptPropKDInit !---------------------------------------------------------------------------- subroutine OptPropKDFinalize ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if KDTblNBin = -1 KDTblNBand = -1 KDTblNPress = -1 KDTblNMol = -1 KDTblNTemp = -1 KDTblNVMRMol1 = -1 KDTblNVMRMol2 = -1 KDTblNVMRMol3 = -1 deallocate( ab_KDTblBandWaveNumBnds ) deallocate( ab_KDTblBandBinIndexBnds ) ! deallocate( c_DelWaveNum ) ! deallocate( c_KDTblBinNum ) deallocate( p_KDTblPress ) deallocate( m_KDTblMolNum ) deallocate( t_KDTblTemp ) deallocate( u_KDTblVMRMol1 ) deallocate( v_KDTblVMRMol2 ) deallocate( w_KDTblVMRMol3 ) ! deallocate( pmctuvw_AbsCoef ) deallocate( ct_AtmPFRatio ) deallocate( c_StrPFRatio ) !!$ deallocate( b_BandAveStrPF ) deallocate( b_RatioBandIntFluxtoTotFlux ) deallocate( b_BandAveRayScatCoefNonRadAct ) deallocate( bm_BandAveRayScatCoef ) !!$ deallocate( a_CharNum ) deallocate( a_KDTblPtclNum ) deallocate( r_KDTblPtclRadiusNum ) deallocate( a_KDTblPtclName ) deallocate( a_KDTblPtclNRadius ) deallocate( ra_KDTblPtclRadius ) FlagInited = .false. end subroutine OptPropKDFinalize !---------------------------------------------------------------------------- function OptPropKDGetNWN() result( NWN ) integer :: NWN ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if NWN = KDTblNBin end function OptPropKDGetNWN !---------------------------------------------------------------------------- function OptPropKDGetNPress() result( NPress ) integer :: NPress ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if NPress = KDTblNPress end function OptPropKDGetNPress !---------------------------------------------------------------------------- function OptPropKDGetNMol() result( NMol ) integer :: NMol ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if NMol = KDTblNMol end function OptPropKDGetNMol !---------------------------------------------------------------------------- function OptPropKDGetNBand() result( NBand ) integer :: NBand ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if NBand = KDTblNBand end function OptPropKDGetNBand !---------------------------------------------------------------------------- function OptPropKDGetDelWN( iWaveNum ) result( DelWaveNum ) integer, intent(in) :: iWaveNum real(DP) :: DelWaveNum ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if DelWaveNum = c_DelWaveNum(iWaveNum) end function OptPropKDGetDelWN !---------------------------------------------------------------------------- function OptPropKDGetBandNum( iWaveNum ) result( BandNum ) integer, intent(in) :: iWaveNum integer :: BandNum ! local variables ! integer :: n if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if do n = 1, KDTblNBand if ( ( ab_KDTblBandBinIndexBnds(1,n) <= iWaveNum ) .and. & & ( iWaveNum <= ab_KDTblBandBinIndexBnds(2,n) ) ) & exit end do BandNum = n end function OptPropKDGetBandNum !---------------------------------------------------------------------------- subroutine OptPropKDGetBandWNBnds( & & NBand, & & aa_BandWNBnds & & ) integer , intent(in ) :: NBand real(DP), intent(out) :: aa_BandWNBnds(2,NBand) ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if ! check if ( NBand /= KDTblNBand ) then write( 6, * ) 'Number of bands is inappropriate: ', NBand, KDTblNBand stop end if aa_BandWNBnds = ab_KDTblBandWaveNumBnds end subroutine OptPropKDGetBandWNBnds !---------------------------------------------------------------------------- subroutine OptPropKDGetBandBinIndexBnds( & & NBand, & & ab_BandBinIndexBnds & & ) integer, intent(in ) :: NBand integer, intent(out) :: ab_BandBinIndexBnds(2,NBand) ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if ! check if ( NBand /= KDTblNBand ) then write( 6, * ) 'Number of bands is inappropriate: ', NBand, KDTblNBand stop end if ab_BandBinIndexBnds = ab_KDTblBandBinIndexBnds end subroutine OptPropKDGetBandBinIndexBnds !---------------------------------------------------------------------------- subroutine OptPropKDGetStrBandAvePF( & & iBand, StrFluxTOA, & & StrBandAvePF & & ) integer , intent(in ) :: iBand real(DP), intent(in ) :: StrFluxTOA real(DP), intent(out) :: StrBandAvePF ! local variables ! !!$ StrBandAvePF = b_BandAveStrPF( iBand ) StrBandAvePF = b_RatioBandIntFluxtoTotFlux( iBand ) StrBandAvePF = StrBandAvePF * StrFluxTOA StrBandAvePF = StrBandAvePF & & / ( ab_KDTblBandWaveNumBnds(2,iBand) - ab_KDTblBandWaveNumBnds(1,iBand) ) end subroutine OptPropKDGetStrBandAvePF !---------------------------------------------------------------------------- !!$ subroutine OptPropKDGetTotalStrFlux( & !!$ & TotalStrFlux & !!$ & ) !!$ !!$ real(DP), intent(out) :: TotalStrFlux !!$ !!$ !!$ ! local variables !!$ ! !!$ real(DP) :: DelWaveNum !!$ !!$ integer :: l !!$ !!$ !!$ TotalStrFlux = 0.0d0 !!$ do l = 1, KDTblNBand !!$ DelWaveNum = ab_KDTblBandWaveNumBnds(2,l) - ab_KDTblBandWaveNumBnds(1,l) !!$ TotalStrFlux = TotalStrFlux & !!$ & + b_BandAveStrPF( l ) * DelWaveNum !!$ end do !!$ !!$ !!$ end subroutine OptPropKDGetTotalStrFlux !---------------------------------------------------------------------------- subroutine OptPropKDGetAbsCoefProf( & & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, rm_VMR, & & rm_AbsCoef & & ) integer , intent(in ) :: iWaveNum integer , intent(in ) :: kmax integer , intent(in ) :: NMol integer , intent(in ) :: m_MolNum (1:NMol) real(DP) , intent(in ) :: r_Press (1:kmax) real(DP) , intent(in ) :: r_Temp (1:kmax) real(DP) , intent(in ) :: rm_VMR (1:kmax,1:NMol) real(DP) , intent(out) :: rm_AbsCoef(1:kmax,1:NMol) ! local variables ! integer :: r_KDTblPressIndex (1:kmax) integer :: r_KDTblTempIndex (1:kmax) integer :: r_KDTblVMRMol1Index(1:kmax) integer :: r_KDTblVMRMol2Index(1:kmax) integer :: r_KDTblVMRMol3Index(1:kmax) ! It should be noticed that NIntpol is defined twice in this module! integer , parameter :: NIntpol = 2 !!$ integer , parameter :: NIntpol = 3 real(DP) :: a_Grid (NIntpol) real(DP) :: a_Weight (NIntpol) real(DP) :: p_Weight1(NIntpol) real(DP) :: t_Weight2(NIntpol) real(DP) :: u_Weight3(NIntpol) real(DP) :: v_Weight4(NIntpol) real(DP) :: w_Weight5(NIntpol) real(DP) :: ptuvw_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol) real(DP) :: ptuvwm_Val (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,KDTblNMol) real(DP) :: m_Val (KDTblNMol) integer :: NIntpolMax integer :: NIntpol3 integer :: NIntpol4 integer :: NIntpol5 integer :: IndexVMRMol1Num integer :: IndexVMRMol2Num integer :: IndexVMRMol3Num real(DP) :: Val real(DP) :: Press real(DP) :: Temp real(DP) :: VMRMol1 real(DP) :: VMRMol2 real(DP) :: VMRMol3 integer :: iPress integer :: iTemp integer :: iVMRMol1 integer :: iVMRMol2 integer :: iVMRMol3 integer :: k integer :: l integer :: m integer :: n integer :: l1 integer :: l2 integer :: l3 integer :: l4 integer :: l5 integer :: n1 integer :: n2 integer :: n3 integer :: n4 integer :: n5 if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if ! check size of VMR of argument and number of molecules in table if ( NMol < KDTblNMol ) then write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol stop end if do m = 1, NMol if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then write( 6, * ) 'Given molecular number, ', m_MolNum(m), ', does not match with that in table.' write( 6, * ) ' Argument: ', m_MolNum write( 6, * ) ' Table : ', m_KDTblMolNum stop end if end do if ( KDTblNPress < NIntpol ) then write( 6, * ) 'Array size for Press is not enough for interpolation: ', & & KDTblNPress, NIntpol stop end if if ( KDTblNTemp < NIntpol ) then write( 6, * ) 'Array size for Temp is not enough for interpolation: ', & & KDTblNTemp, NIntpol stop end if if ( KDTblNVMRMol1 > 1 ) then if ( KDTblNVMRMol1 < NIntpol ) then write( 6, * ) 'Array size for VMRMol1 is not enough for interpolation: ', & & KDTblNVMRMol1, NIntpol stop end if end if if ( KDTblNVMRMol2 > 1 ) then if ( KDTblNVMRMol2 < NIntpol ) then write( 6, * ) 'Array size for VMRMol2 is not enough for interpolation: ', & & KDTblNVMRMol2, NIntpol stop end if end if if ( KDTblNVMRMol3 > 1 ) then if ( KDTblNVMRMol3 < NIntpol ) then write( 6, * ) 'Array size for VMRMol3 is not enough for interpolation: ', & & KDTblNVMRMol3, NIntpol stop end if end if do m = 1, NMol if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol1Num) ) exit end do if ( m > NMol ) then write( 6, * ) & & 'Unable to find a molecular number, ', & & m_KDTblMolNum(KDTblIndexVMRMol1Num), ', in MolNum in argument.' stop end if IndexVMRMol1Num = m ! do m = 1, NMol if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol2Num) ) exit end do if ( m > NMol ) then write( 6, * ) & & 'Unable to find a molecular number, ', & & m_KDTblMolNum(KDTblIndexVMRMol2Num), ', in MolNum in argument.' stop end if IndexVMRMol2Num = m ! do m = 1, NMol if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol3Num) ) exit end do if ( m > NMol ) then write( 6, * ) & & 'Unable to find a molecular number, ', & & m_KDTblMolNum(KDTblIndexVMRMol3Num), ', in MolNum in argument.' stop end if IndexVMRMol3Num = m call FindTblIndex( & & 'Press', & & KDTblNPress, p_KDTblPress, NIntpol, kmax, r_Press, & & r_KDTblPressIndex & & ) call FindTblIndex( & & 'Temp', & & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, & & r_KDTblTempIndex & & ) if ( KDTblNVMRMol1 > 1 ) then call FindTblIndex( & & 'VMRMol1', & & KDTblNVMRMol1, u_KDTblVMRMol1, NIntpol, kmax, rm_VMR(:,IndexVMRMol1Num), & & r_KDTblVMRMol1Index & & ) else r_KDTblVMRMol1Index = 1 end if if ( KDTblNVMRMol2 > 1 ) then call FindTblIndex( & & 'VMRMol2', & & KDTblNVMRMol2, v_KDTblVMRMol2, NIntpol, kmax, rm_VMR(:,IndexVMRMol2Num), & & r_KDTblVMRMol2Index & & ) else r_KDTblVMRMol2Index = 1 end if if ( KDTblNVMRMol3 > 1 ) then call FindTblIndex( & & 'VMRMol3', & & KDTblNVMRMol3, w_KDTblVMRMol3, NIntpol, kmax, rm_VMR(:,IndexVMRMol3Num), & & r_KDTblVMRMol3Index & & ) else r_KDTblVMRMol3Index = 1 end if do k = 1, kmax Press = r_Press (k) Temp = r_Temp (k) VMRMol1 = rm_VMR (k,IndexVMRMol1Num) VMRMol2 = rm_VMR (k,IndexVMRMol2Num) VMRMol3 = rm_VMR (k,IndexVMRMol3Num) iPress = r_KDTblPressIndex (k) iTemp = r_KDTblTempIndex (k) iVMRMol1 = r_KDTblVMRMol1Index(k) iVMRMol2 = r_KDTblVMRMol2Index(k) iVMRMol3 = r_KDTblVMRMol3Index(k) do n = 1, NIntpol a_Grid(n) = p_KDTblPress (iPress+n-1) end do Val = Press if ( a_Grid(1) <= 0.0d0 ) then a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do else a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * log( Val / a_Grid(l) ) / log( a_Grid(n) / a_Grid(l) ) end do end do end if p_Weight1 = a_Weight ! do n = 1, NIntpol a_Grid(n) = t_KDTblTemp(iTemp+n-1) end do Val = Temp a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do t_Weight2 = a_Weight ! if ( KDTblNVMRMol1 > 1 ) then do n = 1, NIntpol a_Grid(n) = u_KDTblVMRMol1(iVMRMol1+n-1) end do Val = VMRMol1 if ( a_Grid(1) <= 0.0d0 ) then a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do else a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * log( Val / a_Grid(l) ) / log( a_Grid(n) / a_Grid(l) ) end do end do end if NIntpolMax = NIntpol else a_Weight = 0.0d0 a_Weight(1) = 1.0d0 NIntpolMax = 1 end if u_Weight3 = a_Weight NIntpol3 = NIntpolMax ! if ( KDTblNVMRMol2 > 1 ) then do n = 1, NIntpol a_Grid(n) = v_KDTblVMRMol2(iVMRMol2+n-1) end do Val = VMRMol2 if ( a_Grid(1) <= 0.0d0 ) then a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do else a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * log( Val / a_Grid(l) ) / log( a_Grid(n) / a_Grid(l) ) end do end do end if NIntpolMax = NIntpol else a_Weight = 0.0d0 a_Weight(1) = 1.0d0 NIntpolMax = 1 end if v_Weight4 = a_Weight NIntpol4 = NIntpolMax ! if ( KDTblNVMRMol3 > 1 ) then do n = 1, NIntpol a_Grid(n) = w_KDTblVMRMol3(iVMRMol3+n-1) end do Val = VMRMol3 if ( a_Grid(1) <= 0.0d0 ) then a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do else a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * log( Val / a_Grid(l) ) / log( a_Grid(n) / a_Grid(l) ) end do end do end if NIntpolMax = NIntpol else a_Weight = 0.0d0 a_Weight(1) = 1.0d0 NIntpolMax = 1 end if w_Weight5 = a_Weight NIntpol5 = NIntpolMax ! ptuvw_Weight = 1.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol do n1 = 1, NIntpol ptuvw_Weight(n1,n2,n3,n4,n5) = ptuvw_Weight(n1,n2,n3,n4,n5) & & * p_Weight1(n1) * t_Weight2(n2) & & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5) end do end do end do end do end do !!$ aaa_Weight = 1.0d0 !!$ do n3 = 1, 2 !!$ do n2 = 1, 2 !!$ do n1 = 1, 2 !!$ ! when n1=1, l1=2 !!$ ! when n1=2, l1=1 !!$ l1 = 3-n1 !!$ l2 = 3-n2 !!$ l3 = 3-n3 !!$ aaa_Weight(n1,n2,n3) = aaa_Weight(n1,n2,n3) & !!$ & * ( Press - a_Grid1(l1) ) / ( a_Grid1(2) - a_Grid1(1) ) & !!$ & * ( Temp - a_Grid2(l2) ) / ( a_Grid2(2) - a_Grid2(1) ) & !!$ & * ( VMRH2O - a_Grid3(l3) ) / ( a_Grid3(2) - a_Grid3(1) ) !!$ end do !!$ end do !!$ end do !!$ aaa_Weight = abs( aaa_Weight ) do m = 1, KDTblNMol do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol do n1 = 1, NIntpol l1 = iPress +n1-1 l2 = iTemp +n2-1 l3 = iVMRMol1+n3-1 l4 = iVMRMol2+n4-1 l5 = iVMRMol3+n5-1 !!$ aaaam_Val(n1,n2,n3,n4,m) = pmctvv_AbsCoef(l1,m,iWaveNum,l2,l3,l4) ptuvwm_Val(n1,n2,n3,n4,n5,m) = & & log( max( pmctuvw_AbsCoef(l1,m,iWaveNum,l2,l3,l4,l5), 1.0d-100 ) ) end do end do end do end do end do end do do m = 1, KDTblNMol m_Val(m) = 0.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol do n1 = 1, NIntpol m_Val(m) = m_Val(m) & & + ptuvwm_Val(n1,n2,n3,n4,n5,m) * ptuvw_Weight(n1,n2,n3,n4,n5) end do end do end do end do end do end do m_Val = exp( m_Val ) do m = 1, min( NMol, KDTblNMol ) rm_AbsCoef(k,m) = m_Val(m) end do do m = min( NMol, KDTblNMol )+1, max( NMol, KDTblNMol ) rm_AbsCoef(k,m) = 0.0d0 end do end do do m = 1, NMol select case ( mc_IDAbsCoefType(m,iWaveNum) ) case ( -1 ) ! zero abs. coef. rm_AbsCoef(:,m) = 0.0d0 case ( 0 ) ! abs. coef. per unit abs. molecule rm_AbsCoef(:,m) = rm_AbsCoef(:,m) * rm_VMR(:,m) !!$ case ( 1 ) ! abs. coef. per unit atm. molecule !!$ rm_AbsCoef(:,m) = rm_AbsCoef(:,m) end select end do end subroutine OptPropKDGetAbsCoefProf !---------------------------------------------------------------------------- subroutine OptPropKDGetRayScatCoef( & & iWaveNum, & & NMol, m_MolNum, & & RayScatCoefNonRadAct, m_RayScatCoef & & ) integer , intent(in ) :: iWaveNum integer , intent(in ) :: NMol integer , intent(in ) :: m_MolNum(NMol) real(DP) , intent(out) :: RayScatCoefNonRadAct real(DP) , intent(out) :: m_RayScatCoef(NMol) ! local variables ! integer :: iBand integer :: m if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if ! check size of VMR of argument and number of molecules in table if ( NMol < KDTblNMol ) then write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol stop end if do m = 1, NMol if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then write( 6, * ) 'Given molecular number, ', m_MolNum(m), ', does not match with that in table.' write( 6, * ) ' Argument: ', m_MolNum write( 6, * ) ' Table : ', m_KDTblMolNum stop end if end do iBand = OptPropKDGetBandNum( iWaveNum ) RayScatCoefNonRadAct = b_BandAveRayScatCoefNonRadAct( iBand ) do m = 1, NMol m_RayScatCoef(m) = bm_BandAveRayScatCoef( iBand, m ) end do end subroutine OptPropKDGetRayScatCoef !---------------------------------------------------------------------------- subroutine OptPropKDGetPtclParam( & & iWaveNum, & & kmax, NPtcl, a_PtclName, za_PtclRadius, & & za_QExt, za_SSA, za_AF & & ) integer , intent(in ) :: iWaveNum integer , intent(in ) :: kmax integer , intent(in ) :: NPtcl character(*), intent(in ) :: a_PtclName(NPtcl) real(DP) , intent(in ) :: za_PtclRadius(kmax,NPtcl) real(DP) , intent(out) :: za_QExt(kmax,NPtcl) real(DP) , intent(out) :: za_SSA (kmax,NPtcl) real(DP) , intent(out) :: za_AF (kmax,NPtcl) ! local variables ! integer :: NRadius real(DP) :: Radius real(DP) :: a_X(2) real(DP) :: a_Val(2) integer :: iBand integer :: iPtcl integer :: i integer :: k integer :: l integer :: n if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if iBand = OptPropKDGetBandNum( iWaveNum ) i = iBand do n = 1, NPtcl do iPtcl = 1, KDTblNPtcl if ( trim( a_PtclName(n) ) == trim( a_KDTblPtclName(iPtcl) ) ) exit end do if ( iPtcl > KDTblNPtcl ) then write( 6, * ) trim( a_PtclName(n) ) // ' cannot found.' stop end if NRadius = a_KDTblPtclNRadius(iPtcl) if ( NRadius == 1 ) then l = 1 do k = 1, kmax za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl) za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl) za_AF (k,n) = bra_KDTblBandAveAF (i,l,iPtcl) end do else do k = 1, kmax Radius = za_PtclRadius(k,n) if ( Radius <= ra_KDTblPtclRadius(1,n) ) then l = 1 za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl) za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl) za_AF (k,n) = bra_KDTblBandAveAF (i,l,iPtcl) else if ( ra_KDTblPtclRadius(NRadius,n) <= Radius ) then l = NRadius za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl) za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl) za_AF (k,n) = bra_KDTblBandAveAF (i,l,iPtcl) else do l = 1+1, NRadius if ( Radius <= ra_KDTblPtclRadius(l,n) ) exit end do a_X(1) = ra_KDTblPtclRadius(l-1,iPtcl) a_X(2) = ra_KDTblPtclRadius(l ,iPtcl) ! a_Val(1) = bra_KDTblBandAveQExt(i,l-1,iPtcl) a_Val(2) = bra_KDTblBandAveQExt(i,l ,iPtcl) za_QExt(k,n) = & & ( a_Val(2) - a_Val(1) ) & & / ( a_X (2) - a_X (1) ) & & * ( Radius - a_X (1) ) & & + a_Val(1) ! a_Val(1) = bra_KDTblBandAveSSA (i,l-1,iPtcl) a_Val(2) = bra_KDTblBandAveSSA (i,l ,iPtcl) za_SSA (k,n) = & & ( a_Val(2) - a_Val(1) ) & & / ( a_X (2) - a_X (1) ) & & * ( Radius - a_X (1) ) & & + a_Val(1) ! a_Val(1) = bra_KDTblBandAveAF (i,l-1,iPtcl) a_Val(2) = bra_KDTblBandAveAF (i,l ,iPtcl) za_AF (k,n) = & & ( a_Val(2) - a_Val(1) ) & & / ( a_X (2) - a_X (1) ) & & * ( Radius - a_X (1) ) & & + a_Val(1) end if end do end if end do end subroutine OptPropKDGetPtclParam !---------------------------------------------------------------------------- subroutine OptPropKDGetPFIntedRatio( & & iWaveNum, kmax, r_Temp, SurfTemp, & & r_AtmPFRatio, SurfPFRatio & & ) integer , intent(in ) :: iWaveNum integer , intent(in ) :: kmax real(DP) , intent(in ) :: r_Temp (1:kmax) real(DP) , intent(in ) :: SurfTemp real(DP) , intent(out) :: r_AtmPFRatio(1:kmax) real(DP) , intent(out) :: SurfPFRatio ! local variables ! integer :: r_KDTblTempIndex (1:kmax) ! It should be noticed that NIntpol is defined twice in this module! integer , parameter :: NIntpol = 2 !!$ integer , parameter :: NIntpol = 3 real(DP) :: a_Grid (NIntpol) real(DP) :: a_Val (NIntpol) real(DP) :: a_Weight(NIntpol) real(DP) :: Val real(DP) :: a_SurfTemp (1) integer :: a_KDTblSurfTempIndex(1) real(DP) :: Temp integer :: iTemp integer :: k integer :: l integer :: n if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if call FindTblIndex( & & 'Temp', & & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, & & r_KDTblTempIndex & & ) a_SurfTemp = SurfTemp call FindTblIndex( & & 'SurfTemp', & & KDTblNTemp, t_KDTblTemp, NIntpol, 1, a_SurfTemp, & & a_KDTblSurfTempIndex & & ) do k = 0, kmax if ( k == 0 ) then Temp = a_SurfTemp (1) iTemp = a_KDTblSurfTempIndex(1) else Temp = r_Temp (k) iTemp = r_KDTblTempIndex(k) end if do n = 1, NIntpol a_Grid(n) = t_KDTblTemp(iTemp+n-1) end do a_Weight = 1.0d0 do n = 1, NIntpol do l = 1, NIntpol if ( l == n ) cycle a_Weight(n) = a_Weight(n) & & * ( Temp - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) ) end do end do do n = 1, NIntpol l = iTemp+n-1 a_Val(n) = ct_AtmPFRatio(iWaveNum,l) end do Val = 0.0d0 do n = 1, NIntpol Val = Val + a_Val(n) * a_Weight(n) end do if ( k == 0 ) then SurfPFRatio = Val else r_AtmPFRatio(k) = Val end if end do end subroutine OptPropKDGetPFIntedRatio !---------------------------------------------------------------------------- subroutine OptPropKDGetStrPFIntedRatio( & & iWaveNum, & & StrPFRatio & & ) integer , intent(in ) :: iWaveNum real(DP) , intent(out) :: StrPFRatio ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if StrPFRatio = c_StrPFRatio(iWaveNum) end subroutine OptPropKDGetStrPFIntedRatio !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- subroutine FindTblIndex( & & AxisName, & & NTblGrid, a_TblGrid, NIntpol, NGrid, a_Grid, & & a_Index & & ) character(*), intent(in ) :: AxisName integer , intent(in ) :: NTblGrid real(DP) , intent(in ) :: a_TblGrid (1:NTblGrid) integer , intent(in ) :: NIntpol integer , intent(in ) :: NGrid real(DP) , intent(in ) :: a_Grid (1:NGrid) integer , intent(out) :: a_Index (1:NGrid) ! local variables ! integer :: i integer :: j ! check if ( a_TblGrid(1) > a_TblGrid(2) ) then write( 6, * ) 'Unexpected order of axis ', trim( AxisName ), & & ', i.e., descending order.' stop end if if ( NTblGrid-NIntpol+1 <= 0 ) then write( 6, * ) & & 'ERROR: Number of points in table axis is insufficient: ', & & NTblGrid, NIntpol end if do i = 1, NGrid if ( ( a_Grid(i) < a_TblGrid(1) ) .or. & & ( a_TblGrid(NTblGrid) < a_Grid(i) ) ) then write( 6, * ) trim( AxisName ), " grid is out of range of a table." write( 6, * ) a_TblGrid(1), a_Grid(i), a_TblGrid(NTblGrid) stop end if do j = 1+1, NTblGrid if ( a_TblGrid(j) >= a_Grid(i) ) exit end do a_Index(i) = j a_Index(i) = a_Index(i) - 1 a_Index(i) = min( a_Index(i), NTblGrid-NIntpol+1 ) end do end subroutine FindTblIndex !---------------------------------------------------------------------------- end module opt_prop_kd