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 public :: OptPropKDPFTblGetPFDPFDT interface OptPropKDPFTblGetPFDPFDT module procedure & OptPropKDPFTblGetPFDPFDT2DV2, & OptPropKDPFTblGetPFDPFDT1DV2, & OptPropKDPFTblGetPFDPFDT1D, & OptPropKDPFTblGetPFDPFDT0D end interface OptPropKDPFTblGetPFDPFDT 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 :: pctuvw_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 :: pmctuvw_LogAbsCoef(:,:,:,:,:,:,:) real(DP), allocatable, save :: b_BandAveRayScatCoefNonRadAct(:) real(DP), allocatable, save :: bm_BandAveRayScatCoef(:,:) !!$ integer :: NChar !!$ integer , allocatable :: a_CharNum(:) integer , save :: KDTblNPtcl integer , save :: KDTblPtclNRadiusMax integer , allocatable, save :: a_KDTblPtclNum (:) integer , allocatable, save :: r_KDTblPtclRadiusNum(:) character(stdstr), allocatable, save :: a_KDTblPtclName (:) integer , allocatable, save :: a_KDTblPtclNRadius (:) real(DP) , allocatable, save :: ra_KDTblPtclRadius (:,:) real(DP) , allocatable, save :: bra_KDTblBandAveQExt(:,:,:) real(DP) , allocatable, save :: bra_KDTblBandAveSSA (:,:,:) real(DP) , allocatable, save :: bra_KDTblBandAveAF (:,:,:) integer , allocatable, save :: mc_IDAbsCoefType(:,:) integer , save :: PFTblNTemp integer , save :: PFTblNBand real(DP), allocatable, save :: t_PFTblTemp (:) !!$ real(DP), allocatable, save :: b_PFTblBandNum(:) real(DP), allocatable, save :: tb_PFTblPF (:,:) real(DP), allocatable, save :: tb_PFTblDPFDT (:,:) 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( pmctuvw_LogAbsCoef(KDTblNPress,KDTblNMol,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) ) call ni3_get_var( InNcID, AbsCoefVarName, pmctuvw_LogAbsCoef ) pmctuvw_LogAbsCoef = log( max( pmctuvw_LogAbsCoef, 1.0d-100 ) ) ! 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( pctuvw_AtmPFRatio(KDTblNPress,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) ) call ni3_get_var( InNcID, AtmPFRatioVarName, pctuvw_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 ) ! Read Planck function table ! call ni3_inq_dimlen( InNcID, "PFTblTemp", PFTblNTemp ) PFTblNBand = KDTblNBand allocate( t_PFTblTemp (PFTblNTemp) ) !!$ allocate( b_PFTblBandNum (PFTblNBand) ) allocate( tb_PFTblPF (PFTblNTemp, KDTblNBand) ) allocate( tb_PFTblDPFDT (PFTblNTemp, KDTblNBand) ) call ni3_get_var( InNcID, "PFTblTemp" , t_PFTblTemp ) !!$ call ni3_get_var( InNcID, "PFTblBandNum", b_PFTblBandNum ) call ni3_get_var( InNcID, "PFTblPF" , tb_PFTblPF ) call ni3_get_var( InNcID, "PFTblDPFDT" , tb_PFTblDPFDT ) 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 ! Preparation of tables for Planck function ! !!$ call OptPropKDPFTblInit 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( pmctuvw_LogAbsCoef ) !!$ deallocate( ct_AtmPFRatio ) deallocate( pctuvw_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 ) ! Finalize tables for Planck function ! call OptPropKDPFTblFinalize 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 :: NIntpol1 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 ! call CalcWeight( & & KDTblNPress, p_KDTblPress, iPress, & & Press, & & NIntpol, & & p_Weight1, NIntpol1 & & ) ! 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 ! call CalcWeight( & & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, & & VMRMol1, & & NIntpol, & & u_Weight3, NIntpol3 & & ) ! !!$ 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 ! call CalcWeight( & & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, & & VMRMol2, & & NIntpol, & & v_Weight4, NIntpol4 & & ) ! !!$ 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 ! call CalcWeight( & & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, & & VMRMol3, & & NIntpol, & & w_Weight5, NIntpol5 & & ) ! ptuvw_Weight = 1.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol do n1 = 1, NIntpol1 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 ) ) ptuvwm_Val(n1,n2,n3,n4,n5,m) = & & pmctuvw_LogAbsCoef(l1,m,iWaveNum,l2,l3,l4,l5) 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 CalcWeight( & & KDTblNAxis, a_KDTblAxis, iAxis, & & TargetValue, & & NIntpol, & & a_Weight, NIntpolMax & & ) integer, intent(in ) :: KDTblNAxis real(8), intent(in ) :: a_KDTblAxis(KDTblNAxis) integer, intent(in ) :: iAxis real(8), intent(in ) :: TargetValue integer, intent(in ) :: NIntpol real(8), intent(out) :: a_Weight(NIntpol) integer, intent(out) :: NIntpolMax ! Local variables ! real(8) :: a_Grid(NIntpol) real(8) :: Val integer :: l integer :: n if ( KDTblNAxis > 1 ) then do n = 1, NIntpol a_Grid(n) = a_KDTblAxis(iAxis+n-1) end do Val = TargetValue 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 end subroutine CalcWeight !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- 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) ) // ' is not 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 OLDOptPropKDGetPFIntedRatio( & !!$ & 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 OLDOptPropKDGetPFIntedRatio !---------------------------------------------------------------------------- subroutine OptPropKDGetPFIntedRatio( & & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, SurfTemp, rm_VMR, & & r_AtmPFRatio, SurfPFRatio & & ) integer , intent(in ) :: iWaveNum integer , intent(in ) :: kmax integer , intent(in ) :: NMol integer , intent(in ) :: m_MolNum(NMol) real(DP) , intent(in ) :: r_Press (1:kmax) real(DP) , intent(in ) :: r_Temp (1:kmax) real(DP) , intent(in ) :: SurfTemp real(DP) , intent(in ) :: rm_VMR (1:kmax,1:NMol) real(DP) , intent(out) :: r_AtmPFRatio(1:kmax) real(DP) , intent(out) :: SurfPFRatio ! 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) :: ptuvw_Val (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol) integer :: NIntpol1 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) call CalcWeight( & & KDTblNPress, p_KDTblPress, iPress, & & Press, & & NIntpol, & & p_Weight1, NIntpol1 & & ) ! 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 ! call CalcWeight( & & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, & & VMRMol1, & & NIntpol, & & u_Weight3, NIntpol3 & & ) ! call CalcWeight( & & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, & & VMRMol2, & & NIntpol, & & v_Weight4, NIntpol4 & & ) ! call CalcWeight( & & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, & & VMRMol3, & & NIntpol, & & w_Weight5, NIntpol5 & & ) ! ptuvw_Weight = 1.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 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 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 l1 = iPress +n1-1 l2 = iTemp +n2-1 l3 = iVMRMol1+n3-1 l4 = iVMRMol2+n4-1 l5 = iVMRMol3+n5-1 ptuvw_Val(n1,n2,n3,n4,n5) = & & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5) end do end do end do end do end do Val = 0.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 Val = Val & & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5) end do end do end do end do end do r_AtmPFRatio(k) = Val end do call FindTblIndex( & & 'Temp', & & KDTblNTemp, t_KDTblTemp, NIntpol, 1, (/SurfTemp/), & & r_KDTblTempIndex(1:1) & & ) do k = 1, 1 Press = r_Press (k) Temp = SurfTemp 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) call CalcWeight( & & KDTblNPress, p_KDTblPress, iPress, & & Press, & & NIntpol, & & p_Weight1, NIntpol1 & & ) ! 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 ! call CalcWeight( & & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, & & VMRMol1, & & NIntpol, & & u_Weight3, NIntpol3 & & ) ! call CalcWeight( & & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, & & VMRMol2, & & NIntpol, & & v_Weight4, NIntpol4 & & ) ! call CalcWeight( & & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, & & VMRMol3, & & NIntpol, & & w_Weight5, NIntpol5 & & ) ! ptuvw_Weight = 1.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 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 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 l1 = iPress +n1-1 l2 = iTemp +n2-1 l3 = iVMRMol1+n3-1 l4 = iVMRMol2+n4-1 l5 = iVMRMol3+n5-1 ptuvw_Val(n1,n2,n3,n4,n5) = & & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5) end do end do end do end do end do Val = 0.0d0 do n5 = 1, NIntpol5 do n4 = 1, NIntpol4 do n3 = 1, NIntpol3 do n2 = 1, NIntpol !!$ do n1 = 1, NIntpol do n1 = 1, NIntpol1 Val = Val & & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5) end do end do end do end do end do SurfPFRatio = Val 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 OptPropKDPFTblInit use constants0, only : PI use rad_kd_utils, only : InqPFIntGQNum use planck_func_wrapper, only : & & Integ_PF_GQ_Array1D , & & Integ_DPFDT_GQ_Array1D ! Local variables ! real(dp), parameter :: DelTemp = 1.0d0 real(dp) :: WNs real(dp) :: WNe integer :: GQNum integer :: iBand integer :: iTemp write( 6, * ) 'Preparing Planck function table ...' ! Last +1 is preparation for quadratic interpolation. PFTblNTemp = int( ( t_KDTblTemp(KDTblNTemp) - t_KDTblTemp(1) ) / DelTemp ) + 1 + 1 PFTblNBand = KDTblNBand allocate( t_PFTblTemp (PFTblNTemp) ) !!$ allocate( b_PFTblBandNum(PFTblNBand) ) allocate( tb_PFTblPF (PFTblNTemp, PFTblNBand) ) allocate( tb_PFTblDPFDT (PFTblNTemp, PFTblNBand) ) do iTemp = 1, PFTblNTemp t_PFTblTemp(iTemp) = t_KDTblTemp(1) + DelTemp * (iTemp-1) end do do iBand = 1, PFTblNBand WNs = ab_KDTblBandWaveNumBnds(1,iBand) WNe = ab_KDTblBandWaveNumBnds(2,iBand) ! GQNum = InqPFIntGQNum( WNe - WNs ) ! call Integ_PF_GQ_Array1D( & & WNs, WNe, GQNum, & & 1, PFTblNTemp, & & t_PFTblTemp, & & tb_PFTblPF(:,iBand) & & ) tb_PFTblPF(:,iBand) = PI * tb_PFTblPF(:,iBand) / ( WNe - WNs ) call Integ_DPFDT_GQ_Array1D( & & WNs, WNe, GQNum, & ! (in ) & 1, PFTblNTemp, & & t_PFTblTemp, & & tb_PFTblDPFDT(:,iBand) & & ) tb_PFTblDPFDT(:,iBand) = PI * tb_PFTblDPFDT(:,iBand) / ( WNe - WNs ) end do end subroutine OptPropKDPFTblInit !---------------------------------------------------------------------------- subroutine OptPropKDPFTblFinalize ! Local variables ! deallocate( t_PFTblTemp ) !!$ deallocate( b_PFTblBandNum ) deallocate( tb_PFTblPF ) deallocate( tb_PFTblDPFDT ) end subroutine OptPropKDPFTblFinalize !---------------------------------------------------------------------------- subroutine OptPropKDPFTblGetPFDPFDT2DV2( & & iBand, & & imax, & & is, ie, & & NLev, xa_Temp, & & xa_PF, & & FlagDPFDT & & ) integer , intent(in ) :: iBand integer , intent(in ) :: imax integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: NLev real(dp), intent(in ) :: xa_Temp(imax,NLev) real(dp), intent(out) :: xa_PF (imax,NLev) logical , intent(in ), optional :: FlagDPFDT ! Local variables ! real(dp) :: a_Temp(NLev) real(dp) :: a_PF (NLev) real(dp) :: DelTemp integer :: iTemp integer :: a_iTemp(NLev) real(dp) :: XVal real(dp) :: a_X(3) real(dp) :: a_Y(3) integer :: i integer :: k do i = is, ie do k = 1, NLev a_Temp(k) = xa_Temp(i,k) end do DelTemp = t_PFTblTemp(2) - t_PFTblTemp(1) do k = 1, NLev if ( ( a_Temp(k) < t_PFTblTemp(1) ) & & .or. ( t_PFTblTemp(PFTblNTemp) < a_Temp(k) ) ) then write( 6, * ) 'PF Interpolation: Temperature is out of range:' write( 6, * ) ' ', k, a_Temp(k) write( 6, * ) ' ', t_PFTblTemp(1), t_PFTblTemp(PFTblNTemp) stop end if end do iTemp = int( ( a_Temp(k) - t_PFTblTemp(1) ) / DelTemp ) + 1 + 1 iTemp = max( iTemp, 2 ) iTemp = min( iTemp, PFTblNTemp-1 ) if ( ( a_Temp(k) < t_PFTblTemp(iTemp-1) ) & & .or. ( t_PFTblTemp(iTemp+1) < a_Temp(k) ) ) then write( 6, * ) 'PF Interpolation: Deduced position is out of range:' write( 6, * ) ' ', k, a_Temp(k) write( 6, * ) ' ', t_PFTblTemp(iTemp-1), t_PFTblTemp(iTemp+1) stop end if a_iTemp(k) = iTemp if ( present( FlagDPFDT ) .and. FlagDPFDT ) then do k = 1, NLev iTemp = a_iTemp(k) !!$ a_PF(k) = & !!$ & ( tb_PFTblDPFDT(iTemp,iBand) - tb_PFTblDPFDT(iTemp-1,iBand) ) & !!$ & / ( t_PFTblTemp(iTemp) - t_PFTblTemp(iTemp-1) ) & !!$ & * ( a_Temp(k) - t_PFTblTemp(iTemp-1) ) & !!$ & + tb_PFTblDPFDT(iTemp-1,iBand) XVal = a_Temp(k) a_X(1) = t_PFTblTemp(iTemp-1) a_X(2) = t_PFTblTemp(iTemp ) a_X(3) = t_PFTblTemp(iTemp+1) a_Y(1) = tb_PFTblDPFDT(iTemp-1,iBand) a_Y(2) = tb_PFTblDPFDT(iTemp ,iBand) a_Y(3) = tb_PFTblDPFDT(iTemp+1,iBand) a_PF(k) = & & ( ( XVal - a_X(2) ) * ( XVal - a_X(3) ) ) & & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) & & * a_Y(1) & & + ( ( XVal - a_X(3) ) * ( XVal - a_X(1) ) ) & & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) & & * a_Y(2) & & + ( ( XVal - a_X(1) ) * ( XVal - a_X(2) ) ) & & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) & & * a_Y(3) end do else do k = 1, NLev iTemp = a_iTemp(k) !!$ a_PF(k) = & !!$ & ( tb_PFTblPF(iTemp,iBand) - tb_PFTblPF(iTemp-1,iBand) ) & !!$ & / ( t_PFTblTemp(iTemp) - t_PFTblTemp(iTemp-1) ) & !!$ & * ( a_Temp(k) - t_PFTblTemp(iTemp-1) ) & !!$ & + tb_PFTblPF(iTemp-1,iBand) XVal = a_Temp(k) a_X(1) = t_PFTblTemp(iTemp-1) a_X(2) = t_PFTblTemp(iTemp ) a_X(3) = t_PFTblTemp(iTemp+1) a_Y(1) = tb_PFTblPF(iTemp-1,iBand) a_Y(2) = tb_PFTblPF(iTemp ,iBand) a_Y(3) = tb_PFTblPF(iTemp+1,iBand) a_PF(k) = & & ( ( XVal - a_X(2) ) * ( XVal - a_X(3) ) ) & & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) & & * a_Y(1) & & + ( ( XVal - a_X(3) ) * ( XVal - a_X(1) ) ) & & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) & & * a_Y(2) & & + ( ( XVal - a_X(1) ) * ( XVal - a_X(2) ) ) & & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) & & * a_Y(3) end do end if do k = 1, NLev xa_PF(i,k) = a_PF(k) end do end do end subroutine OptPropKDPFTblGetPFDPFDT2DV2 !---------------------------------------------------------------------------- subroutine OptPropKDPFTblGetPFDPFDT1DV2( & & iBand, & & imax, & & is, ie, & & x_Temp, & & x_PF, & & FlagDPFDT & & ) integer , intent(in ) :: iBand integer , intent(in ) :: imax integer , intent(in ) :: is integer , intent(in ) :: ie real(dp), intent(in ) :: x_Temp(imax) real(dp), intent(out) :: x_PF (imax) logical , intent(in ), optional :: FlagDPFDT ! Local variables ! real(dp) :: xa_Temp(imax,1) real(dp) :: xa_PF (imax,1) integer :: i do i = is, ie xa_Temp(i,1) = x_Temp(i) end do call OptPropKDPFTblGetPFDPFDT( & & iBand, & & imax, & & is, ie, & & 1, xa_Temp, & & xa_PF, & & FlagDPFDT & & ) do i = is, ie x_PF(i) = xa_PF(i,1) end do end subroutine OptPropKDPFTblGetPFDPFDT1DV2 !---------------------------------------------------------------------------- subroutine OptPropKDPFTblGetPFDPFDT1D( & & iBand, NLev, a_Temp, & & a_PF, & & FlagDPFDT & & ) integer , intent(in ) :: iBand integer , intent(in ) :: NLev real(dp), intent(in ) :: a_Temp(NLev) real(dp), intent(out) :: a_PF (NLev) logical , intent(in ), optional :: FlagDPFDT ! Local variables ! real(dp) :: DelTemp integer :: iTemp integer :: a_iTemp(NLev) real(dp) :: XVal real(dp) :: a_X(3) real(dp) :: a_Y(3) integer :: k DelTemp = t_PFTblTemp(2) - t_PFTblTemp(1) do k = 1, NLev if ( ( a_Temp(k) < t_PFTblTemp(1) ) & & .or. ( t_PFTblTemp(PFTblNTemp) < a_Temp(k) ) ) then write( 6, * ) 'PF Interpolation: Temperature is out of range:' write( 6, * ) ' ', k, a_Temp(k) write( 6, * ) ' ', t_PFTblTemp(1), t_PFTblTemp(PFTblNTemp) stop end if iTemp = int( ( a_Temp(k) - t_PFTblTemp(1) ) / DelTemp ) + 1 + 1 iTemp = max( iTemp, 2 ) iTemp = min( iTemp, PFTblNTemp-1 ) if ( ( a_Temp(k) < t_PFTblTemp(iTemp-1) ) & & .or. ( t_PFTblTemp(iTemp+1) < a_Temp(k) ) ) then write( 6, * ) 'PF Interpolation: Deduced position is out of range:' write( 6, * ) ' ', k, a_Temp(k) write( 6, * ) ' ', t_PFTblTemp(iTemp-1), t_PFTblTemp(iTemp+1) stop end if a_iTemp(k) = iTemp end do if ( present( FlagDPFDT ) .and. FlagDPFDT ) then do k = 1, NLev iTemp = a_iTemp(k) !!$ a_PF(k) = & !!$ & ( tb_PFTblDPFDT(iTemp,iBand) - tb_PFTblDPFDT(iTemp-1,iBand) ) & !!$ & / ( t_PFTblTemp(iTemp) - t_PFTblTemp(iTemp-1) ) & !!$ & * ( a_Temp(k) - t_PFTblTemp(iTemp-1) ) & !!$ & + tb_PFTblDPFDT(iTemp-1,iBand) XVal = a_Temp(k) a_X(1) = t_PFTblTemp(iTemp-1) a_X(2) = t_PFTblTemp(iTemp ) a_X(3) = t_PFTblTemp(iTemp+1) a_Y(1) = tb_PFTblDPFDT(iTemp-1,iBand) a_Y(2) = tb_PFTblDPFDT(iTemp ,iBand) a_Y(3) = tb_PFTblDPFDT(iTemp+1,iBand) a_PF(k) = & & ( ( XVal - a_X(2) ) * ( XVal - a_X(3) ) ) & & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) & & * a_Y(1) & & + ( ( XVal - a_X(3) ) * ( XVal - a_X(1) ) ) & & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) & & * a_Y(2) & & + ( ( XVal - a_X(1) ) * ( XVal - a_X(2) ) ) & & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) & & * a_Y(3) end do else do k = 1, NLev iTemp = a_iTemp(k) !!$ a_PF(k) = & !!$ & ( tb_PFTblPF(iTemp,iBand) - tb_PFTblPF(iTemp-1,iBand) ) & !!$ & / ( t_PFTblTemp(iTemp) - t_PFTblTemp(iTemp-1) ) & !!$ & * ( a_Temp(k) - t_PFTblTemp(iTemp-1) ) & !!$ & + tb_PFTblPF(iTemp-1,iBand) XVal = a_Temp(k) a_X(1) = t_PFTblTemp(iTemp-1) a_X(2) = t_PFTblTemp(iTemp ) a_X(3) = t_PFTblTemp(iTemp+1) a_Y(1) = tb_PFTblPF(iTemp-1,iBand) a_Y(2) = tb_PFTblPF(iTemp ,iBand) a_Y(3) = tb_PFTblPF(iTemp+1,iBand) a_PF(k) = & & ( ( XVal - a_X(2) ) * ( XVal - a_X(3) ) ) & & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) & & * a_Y(1) & & + ( ( XVal - a_X(3) ) * ( XVal - a_X(1) ) ) & & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) & & * a_Y(2) & & + ( ( XVal - a_X(1) ) * ( XVal - a_X(2) ) ) & & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) & & * a_Y(3) end do end if end subroutine OptPropKDPFTblGetPFDPFDT1D !---------------------------------------------------------------------------- subroutine OptPropKDPFTblGetPFDPFDT0D( & & iBand, Temp, & & PF, & & FlagDPFDT & & ) integer , intent(in ) :: iBand real(dp), intent(in ) :: Temp real(dp), intent(out) :: PF logical , intent(in ), optional :: FlagDPFDT ! Local variables ! integer , parameter :: NLev = 1 real(dp) :: a_Temp(NLev) real(dp) :: a_PF (NLev) a_Temp = Temp call OptPropKDPFTblGetPFDPFDT1D( & & iBand, NLev, a_Temp, & & a_PF, & & FlagDPFDT & & ) PF = a_PF(1) end subroutine OptPropKDPFTblGetPFDPFDT0D !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- 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