module rad_flux use vtype_module implicit none private public :: RadFluxCalcFlux public :: RadFluxInit public :: RadFluxFinalize public :: RadFluxOutputInit public :: RadFluxOutputFinalize logical, parameter :: FlagRayleigh = .true. !!$ logical, parameter :: FlagRayleigh = .false. integer, save :: IDMethod integer, parameter :: IDMethodKD = 0 integer, parameter :: IDMethodKDDevelop = 1 integer, parameter :: IDMethodLBL = 2 integer, parameter :: IDMethodLBLConstPF = 3 real(8), save :: Grav integer , save :: NBand real(DP), allocatable, save :: aa_BandWNBnds (:,:) integer , allocatable, save :: aa_BandBinIndexBnds(:,:) logical , save :: FlagOutputInited = .false. integer , save :: NcIDFlux real(DP) , save :: StrFluxTOA !!$ real(DP) , save :: FluxRatTOAtoStelSfc integer , save :: iWaveNumMaxforAlbedoLW integer , save :: iWaveNumMaxforMethodLW character(128), save :: ModuleName = 'flux' logical, save :: FlagInited = .false. !---------------------------------------------------------------------------- contains !---------------------------------------------------------------------------- subroutine RadFluxCalcFlux( & & kmax, NMolNum, & & r_Press, & & SurfTemp, r_Temp, rm_VMR, & & NPtcl, a_PtclName, za_PtclEffRad, za_DelPtclNum, & & AlbedoSW, AlbedoLW, SZADeg, & & r_MeanMolWt, & & r_RadUwFlux, r_RadDwFlux, & & ra_DelRadUwFlux, ra_DelRadDwFlux, & & WNIntegStart, WNIntegEnd & & ) use ni3_module use constants0, only : & & PI, & & AtomMassConst, & & AvogConst use rad_kd_utils, only : InqPFIntGQNum use rad_rte_two_stream_app_1d_wrapper, only : & & RadRTETwoStreamApp1DLWWrapper, & & RadRTETwoStreamApp1DSWWrapper, & & RadRTETwoStreamApp1DWrapper, & & IDScatApproxEddington, & & IDScatApproxHemiMean use opt_prop, only : & & OptPropGetNMol, & & OptPropGetNWN, & & OptPropGetNBand, & & OptPropGetDelWN, & & OptPropGetWN, & & OptPropGetBandNum, & & OptPropGetAbsCoefProf, & & OptPropGetRayScatCoef, & & OptPropGetPtclParam, & & OptPropGetPFInted, & & OptPropGetStrPFInted, & & OptPropGetStrBandAvePF, & & OptPropGetPFIntedRatio, & & OptPropGetStrPFIntedRatio use planck_func_wrapper integer , intent(in ) :: kmax integer , intent(in ) :: NMolNum real(DP), intent(in ) :: r_Press (0:kmax) real(DP), intent(in ) :: SurfTemp real(DP), intent(in ) :: r_Temp (0:kmax) real(DP), intent(in ) :: rm_VMR (0:kmax,1:NMolNum) integer , intent(in ) :: NPtcl character(*), intent(in ) :: a_PtclName (NPtcl) real(8) , intent(in ) :: za_PtclEffRad(1:kmax,NPtcl) real(8) , intent(in ) :: za_DelPtclNum(1:kmax,NPtcl) real(DP), intent(in ) :: AlbedoSW real(DP), intent(in ) :: AlbedoLW real(DP), intent(in ) :: SZADeg real(DP), intent(in ) :: r_MeanMolWt (0:kmax) real(DP), intent(out) :: r_RadUwFlux (0:kmax) real(DP), intent(out) :: r_RadDwFlux (0:kmax) real(DP), intent(out) :: ra_DelRadUwFlux(0:kmax,0:1) real(DP), intent(out) :: ra_DelRadDwFlux(0:kmax,0:1) real(DP), intent(in ) :: WNIntegStart real(DP), intent(in ) :: WNIntegEnd ! local variables ! real(DP), allocatable :: a_BandAveSurfPF (:) real(DP), allocatable :: ra_BandAvePF (:,:) real(DP), allocatable :: a_BandAveSurfDPFDT(:) real(DP), allocatable :: a_BandAveStrFluxTOA(:) integer :: NMolTmp integer :: NWaveNum integer :: iBand integer :: iBandStart integer :: iBandEnd integer :: iWaveNum integer :: iWaveNumStart integer :: iWaveNumEnd real(8) :: WaveNum real(8) :: DelWaveNum real(8), allocatable :: pm_AbsCoefPerAtmMol(:,:) real(DP) :: RayScatCoefNonRadAct real(DP) :: m_RayScatCoef(NMolNum) real(DP) :: GasAbsCoef real(DP) :: RayScatCoef real(DP) :: DelGasAbsOptDep real(DP) :: DelRayScatOptDep real(DP) :: a_DelPtclExtOptDep(NPtcl) real(dp) :: za_PtclQExt(1:kmax, NPtcl) real(dp) :: za_PtclSSA (1:kmax, NPtcl) real(dp) :: za_PtclAF (1:kmax, NPtcl) real(8) :: z_DelAtmMols (1:kmax) real(8) :: DelOptDep real(DP) :: z_SSA (1:kmax) real(DP) :: z_AF (1:kmax) real(DP) :: r_OptDep (0:kmax) real(DP) :: r_PFInted (0:kmax) real(DP) :: SurfPFInted real(DP) :: SurfDPFDTInted real(DP) :: r_RadUwFluxEach (0:kmax) real(DP) :: r_RadDwFluxEach (0:kmax) real(DP) :: ra_DelRadUwFluxEach(0:kmax,0:1) real(DP) :: ra_DelRadDwFluxEach(0:kmax,0:1) real(DP) :: r_RadUwFluxEachLW (0:kmax) real(DP) :: r_RadDwFluxEachLW (0:kmax) real(DP) :: ra_DelRadUwFluxEachLW(0:kmax,0:1) real(DP) :: ra_DelRadDwFluxEachLW(0:kmax,0:1) real(DP) :: r_RadUwFluxEachSW (0:kmax) real(DP) :: r_RadDwFluxEachSW (0:kmax) real(DP) :: ra_DelRadUwFluxEachSW(0:kmax,0:1) real(DP) :: ra_DelRadDwFluxEachSW(0:kmax,0:1) real(DP) :: z_DRadFluxDPressEach(1:kmax) real(DP) :: r_PFIntedRatio (0:kmax) real(DP) :: SurfPFIntedRatio real(DP) :: StrPFRatio integer :: IDScatApprox logical :: FlagSrcFuncTech real(DP) :: StrFluxTOAEach real(DP) :: InAngle real(DP) :: SurfAlbedo !!$ real(DP) :: TempMin real(DP) :: VMRNonRadAct real(DP) :: WNs real(DP) :: WNe integer :: GQNum integer :: k integer :: m integer :: n if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if select case ( IDMethod ) case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD ) allocate( a_BandAveSurfPF (NBand) ) allocate( ra_BandAvePF (0:kmax,NBand) ) allocate( a_BandAveSurfDPFDT (NBand) ) allocate( a_BandAveStrFluxTOA(NBand) ) do n = 1, NBand WNs = aa_BandWNBnds(1,n) WNe = aa_BandWNBnds(2,n) ! GQNum = InqPFIntGQNum( WNe - WNs ) ! call Integ_PF_GQ_Array0D( & & WNs, WNe, GQNum, & & SurfTemp, & & a_BandAveSurfPF(n) & & ) a_BandAveSurfPF(n) = PI * a_BandAveSurfPF(n) / ( WNe - WNs ) call Integ_PF_GQ_Array1D( & & WNs, WNe, GQNum, & & 0, kmax, & & r_Temp, & & ra_BandAvePF(:,n) & & ) ra_BandAvePF(:,n) = PI * ra_BandAvePF(:,n) / ( WNe - WNs ) call Integ_DPFDT_GQ_Array0D( & & WNs, WNe, GQNum, & ! (in ) & SurfTemp, & ! (in ) & a_BandAveSurfDPFDT(n) & & ) a_BandAveSurfDPFDT(n) = PI * a_BandAveSurfDPFDT(n) / ( WNe - WNs ) end do do n = 1, NBand call OptPropGetStrBandAvePF( & & n, StrFluxTOA, & & a_BandAveStrFluxTOA(n) & & ) end do end select do k = 1, kmax !!$ z_DelAtmMols(k) = & !!$ & ( r_Press(k-1) - r_Press(k) ) / Grav z_DelAtmMols(k) = & & ( r_Press(k-1) - r_Press(k) ) / Grav & & / ( ( r_MeanMolWt(k-1) + r_MeanMolWt(k) ) / 2.0d0 & & / AvogConst ) end do NMolTmp = OptPropGetNMol() if ( NMolTmp /= NMolNum ) then write( 6, * ) 'NMolNum in argument does not agree with NMol in data.' write( 6, * ) NMolNum, NMolTmp stop end if NWaveNum = OptPropGetNWN() allocate( pm_AbsCoefPerAtmMol(0:kmax, NMolNum) ) r_RadUwFlux = 0.0d0 r_RadDwFlux = 0.0d0 ra_DelRadUwFlux = 0.0d0 ra_DelRadDwFlux = 0.0d0 select case ( IDMethod ) case ( IDMethodKDDevelop, IDMethodKD ) if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then iWaveNumStart = 1 iWaveNumEnd = NWaveNum else do iBand = 1, NBand if ( aa_BandWNBnds(1,iBand) >= WNIntegStart ) exit end do iBandStart = iBand iBandStart = min( iBandStart, NBand ) do iBand = iBandStart, NBand if ( aa_BandWNBnds(2,iBand) >= WNIntegEnd ) exit end do iBandEnd = min( iBand, NBand ) iWaveNumStart = aa_BandBinIndexBnds(1,iBandStart) iWaveNumEnd = aa_BandBinIndexBnds(2,iBandEnd ) end if case ( IDMethodLBL, IDMethodLBLConstPF ) if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then iWaveNumStart = 1 iWaveNumEnd = NWaveNum else loop_search_wnstart: do iWaveNum = 1, NWaveNum if ( OptPropGetWN( iWaveNum ) >= WNIntegStart ) exit end do loop_search_wnstart iWaveNumStart = iWaveNum iWaveNumStart = min( iWaveNumStart, NWaveNum ) loop_search_wnend: do iWaveNum = iWaveNumStart, NWaveNum if ( OptPropGetWN( iWaveNum ) > WNIntegEnd ) exit end do loop_search_wnend iWaveNumEnd = iWaveNum - 1 end if case default write( 6, * ) 'In ', trim( ModuleName ) write( 6, * ) ' Unexpected case default.' stop end select !!$ do iWaveNum = 1, NWaveNum do iWaveNum = iWaveNumStart, iWaveNumEnd DelWaveNum = OptPropGetDelWN( iWaveNum ) call OptPropGetAbsCoefProf( & & iWaveNum, kmax+1, NMolNum, r_Press, r_Temp, rm_VMR, & & pm_AbsCoefPerAtmMol & & ) if ( FlagRayleigh ) then call OptPropGetRayScatCoef( & & iWaveNum, & & NMolNum, & & RayScatCoefNonRadAct, m_RayScatCoef & & ) else RayScatCoefNonRadAct = 0.0d0 m_RayScatCoef = 0.0d0 end if call OptPropGetPtclParam( & & iWaveNum, & & kmax, NPtcl, a_PtclName, za_PtclEffRad, & & za_PtclQExt, za_PtclSSA, za_PtclAF & & ) r_OptDep(kmax) = 0.0d0 do k = kmax, 1, -1 GasAbsCoef = 0.0d0 do m = 1, NMolNum GasAbsCoef = GasAbsCoef & & + ( pm_AbsCoefPerAtmMol(k-1,m) & & + pm_AbsCoefPerAtmMol(k ,m) ) / 2.0d0 end do DelGasAbsOptDep = GasAbsCoef * z_DelAtmMols(k) ! VMRNonRadAct = 1.0d0 do m = 1, NMolNum VMRNonRadAct = VMRNonRadAct - ( rm_VMR(k-1,m) + rm_VMR(k,m) ) / 2.0d0 end do VMRNonRadAct = max( VMRNonRadAct, 0.0d0 ) RayScatCoef = RayScatCoefNonRadAct * VMRNonRadAct do m = 1, NMolNum RayScatCoef = RayScatCoef & & + m_RayScatCoef(m) * ( rm_VMR(k-1,m) + rm_VMR(k,m) ) / 2.0d0 end do DelRayScatOptDep = RayScatCoef * z_DelAtmMols(k) ! do m = 1, NPtcl a_DelPtclExtOptDep(m) = & & za_PtclQExt(k,m) * PI * za_PtclEffRad(k,m) * za_DelPtclNum(k,m) end do DelOptDep = DelGasAbsOptDep + DelRayScatOptDep do m = 1, NPtcl DelOptDep = DelOptDep + a_DelPtclExtOptDep(m) end do r_OptDep(k-1) = r_OptDep(k) + DelOptDep ! if ( DelOptDep == 0.0d0 ) then z_SSA(k) = 0.0d0 else z_SSA(k) = 0.0d0 * DelGasAbsOptDep + 1.0d0 * DelRayScatOptDep do m = 1, NPtcl z_SSA(k) = z_SSA(k) + za_PtclSSA(k,m) * a_DelPtclExtOptDep(m) end do z_SSA(k) = z_SSA(k) / DelOptDep end if ! if ( z_SSA(k) == 0.0d0 ) then z_AF(k) = 0.0d0 else z_AF(k) = & & 0.0d0 * 0.0d0 * DelGasAbsOptDep & & + 0.0d0 * 1.0d0 * DelRayScatOptDep do m = 1, NPtcl z_AF(k) = z_AF(k) & & + za_PtclAF(k,m) * za_PtclSSA(k,m) * a_DelPtclExtOptDep(m) end do z_AF(k) = z_AF(k) / ( z_SSA(k) * DelOptDep ) end if end do z_SSA = min( z_SSA, 1.0d0 - 1.0d-10 ) z_AF = min( z_AF , 1.0d0 ) select case ( IDMethod ) case ( IDMethodLBL ) call OptPropGetPFInted( & & iWaveNum, & & kmax+1, & & r_Temp, SurfTemp, & & r_PFInted, SurfPFInted, SurfDPFDTInted & & ) ! call OptPropGetStrPFInted( & & iWaveNum, StrFluxTOA, & & StrFluxTOAEach & & ) case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD ) n = OptPropGetBandNum( iWaveNum ) if ( n >= NBand + 1 ) then write( 6, * ) 'In ', trim( ModuleName ) write( 6, * ) ' Unable to find appropriate band number: ', n stop end if r_PFInted = ra_BandAvePF(:,n) SurfPFInted = a_BandAveSurfPF(n) SurfDPFDTInted = a_BandAveSurfDPFDT(n) ! call OptPropGetPFIntedRatio( & & iWaveNum, kmax+1, r_Temp, SurfTemp, & & r_PFIntedRatio, SurfPFIntedRatio & & ) r_PFInted = r_PFInted * r_PFIntedRatio SurfPFInted = SurfPFInted * SurfPFIntedRatio SurfDPFDTInted = SurfDPFDTInted * SurfPFIntedRatio StrFluxTOAEach = a_BandAveStrFluxTOA(n) ! call OptPropGetStrPFIntedRatio( & & iWaveNum, & & StrPFRatio & & ) StrFluxTOAEach = StrFluxTOAEach * StrPFRatio r_PFInted = r_PFInted * DelWaveNum SurfPFInted = SurfPFInted * DelWaveNum SurfDPFDTInted = SurfDPFDTInted * DelWaveNum ! StrFluxTOAEach = StrFluxTOAEach * DelWaveNum case default write( 6, * ) 'In ', trim( ModuleName ) write( 6, * ) ' Unexpected case default' stop end select if ( iWaveNum <= iWaveNumMaxforAlbedoLW ) then SurfAlbedo = AlbedoLW else SurfAlbedo = AlbedoSW end if if ( iWaveNum <= iWaveNumMaxforMethodLW ) then IDScatApprox = IDScatApproxHemiMean FlagSrcFuncTech = .true. !!$ IDScatApprox = IDScatApproxEddington !!$ FlagSrcFuncTech = .false. else IDScatApprox = IDScatApproxEddington FlagSrcFuncTech = .false. end if !!$ write( 6, * ) iWaveNum, OptPropGetWN( iWaveNum ), SurfAlbedo !!$ write( 6, * ) iWaveNum, SurfAlbedo, iWaveNumMaxforAlbedoLW ! SurfAlbedo = 0.0d0 ! SurfAlbedo = AlbedoLW call RadRTETwoStreamApp1DLWWrapper( & & kmax, & ! (in) & z_SSA, z_AF, & ! (in) & r_OptDep, & ! (in) & SurfAlbedo, & ! (in) & r_PFInted, SurfPFInted, SurfDPFDTInted, & ! (in) ! & r_RadUwFluxEach, r_RadDwFluxEach, & ! (out) ! & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach & ! (out) & r_RadUwFluxEachLW, r_RadDwFluxEachLW, & ! (out) & ra_DelRadUwFluxEachLW, ra_DelRadDwFluxEachLW & ! (out) & ) InAngle = 1.0d0 / cos( SZADeg * PI / 180.0d0 ) call RadRTETwoStreamApp1DSWWrapper( & & kmax, & ! (in) & z_SSA, z_AF, & ! (in) & r_OptDep, & ! (in) & SurfAlbedo, & ! (in) & StrFluxTOAEach, InAngle, & ! (in) ! & r_RadUwFluxEach, r_RadDwFluxEach & ! (out) & r_RadUwFluxEachSW, r_RadDwFluxEachSW & ! (out) & ) ! ra_DelRadUwFluxEach = 0.0d0 ! ra_DelRadDwFluxEach = 0.0d0 ra_DelRadUwFluxEachSW = 0.0d0 ra_DelRadDwFluxEachSW = 0.0d0 r_RadUwFluxEach = r_RadUwFluxEachLW + r_RadUwFluxEachSW r_RadDwFluxEach = r_RadDwFluxEachLW + r_RadDwFluxEachSW ra_DelRadUwFluxEach = ra_DelRadUwFluxEachLW + ra_DelRadUwFluxEachSW ra_DelRadDwFluxEach = ra_DelRadDwFluxEachLW + ra_DelRadDwFluxEachSW !!$ InAngle = 1.0d0 / cos( SZADeg * PI / 180.0d0 ) !!$ !!$ call RadRTETwoStreamApp1DWrapper( & !!$ & IDScatApprox, FlagSrcFuncTech, & ! (in) !!$ & kmax, & ! (in) !!$ & z_SSA, z_AF, & ! (in) !!$ & r_OptDep, & ! (in) !!$ & SurfAlbedo, & ! (in) !!$! & r_RadUwFlux, r_RadDwFlux, & ! (out) !!$ & r_RadUwFluxEach, r_RadDwFluxEach, & ! (out) !!$! & SolarFluxTOA, InAngle, & ! (in) !!$ & StrFluxTOAEach, InAngle, & ! (in) !!$ & r_PFInted, SurfPFInted, SurfDPFDTInted, & ! (in) !!$! & ra_DelRadUwFlux, ra_DelRadDwFlux & ! (out) !!$ & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach & ! (out) !!$ & ) !!$ if ( iWaveNumEnd < 10000 ) then !!$ select case ( IDMethod ) !!$ case ( IDMethodLBLConstPF, IDMethodKDDevelop ) !!$ write( 60, * ) iWaveNum, r_OptDep(0), StrFluxTOAEach, r_RadDwFluxEach(0), pm_AbsCoefPerAtmMol(0,1), StrFluxTOAEach, InAngle !!$ write( 60, * ) iWaveNum, r_PFInted(0), SurfPFInted, SurfDPFDTInted, & !!$ & IDScatApprox, FlagSrcFuncTech !!$ end select !!$ end if if ( FlagOutputInited ) then call ni3_put_varss( NcIDFlux, "RadUwFlux", iWaveNum, r_RadUwFluxEach/DelWaveNum ) call ni3_put_varss( NcIDFlux, "RadDwFlux", iWaveNum, r_RadDwFluxEach/DelWaveNum ) call ni3_put_varss( NcIDFlux, "OptDep" , iWaveNum, r_OptDep ) do k = 1, kmax z_DRadFluxDPressEach(k) = & & ( ( r_RadUwFluxEach(k ) - r_RadDwFluxEach(k ) ) & & - ( r_RadUwFluxEach(k-1) - r_RadDwFluxEach(k-1) ) ) & & / ( r_Press(k) - r_Press(k-1) ) end do z_DRadFluxDPressEach = z_DRadFluxDPressEach / DelWaveNum call ni3_put_varss( NcIDFlux, "DRadFluxDPress", iWaveNum, z_DRadFluxDPressEach ) end if r_RadUwFlux = r_RadUwFlux + r_RadUwFluxEach r_RadDwFlux = r_RadDwFlux + r_RadDwFluxEach ra_DelRadUwFlux = ra_DelRadUwFlux + ra_DelRadUwFluxEach ra_DelRadDwFlux = ra_DelRadDwFlux + ra_DelRadDwFluxEach end do deallocate( pm_AbsCoefPerAtmMol ) select case ( IDMethod ) case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD ) deallocate( a_BandAveSurfPF ) deallocate( ra_BandAvePF ) deallocate( a_BandAveSurfDPFDT ) deallocate( a_BandAveStrFluxTOA ) end select end subroutine RadFluxCalcFlux !---------------------------------------------------------------------------- subroutine RadFluxInit( & & ArgGrav, ArgStrFluxTOA, & & WaveNumMethodSwitch, & & WaveNumAlbedoSwitch, & & OptPropNcFN, & & ArgRayScatCoefNcFN, ArgStrSpeNcFN, & & Method, & & a_ArgBandWNBnds & & ) use constants0, only : StB use opt_prop, only : & & OptPropInit, & & OptPropGetNBand, & & OptPropGetBandWNBnds, & & OptPropGetBandBinIndexBnds !!$ & OptPropGetTotalStrFlux use rad_rte_two_stream_app_1d_wrapper, only : RadRTETwoStreamApp1DWrapperInit real(8) , intent(in ) :: ArgGrav real(DP) , intent(in ) :: ArgStrFluxTOA real(DP) , intent(in ) :: WaveNumMethodSwitch real(DP) , intent(in ) :: WaveNumAlbedoSwitch character(*), intent(in ) :: OptPropNcFN character(*), intent(in ), optional :: ArgRayScatCoefNcFN character(*), intent(in ), optional :: ArgStrSpeNcFN character(*), intent(in ), optional :: Method real(DP) , intent(in ), optional :: a_ArgBandWNBnds(:) ! local variables ! character(extstr) :: RayScatCoefNcFN character(extstr) :: StrSpeNcFN integer :: n if ( FlagInited ) return if ( present( Method ) ) then select case ( Method ) case ( 'KD' ) IDMethod = IDMethodKD case ( 'KDDevelop' ) IDMethod = IDMethodKDDevelop case ( 'LBL' ) IDMethod = IDMethodLBL case ( 'LBLConstPF' ) IDMethod = IDMethodLBLConstPF case default write( 6, * ) "In module, ", trim( ModuleName ) write( 6, * ) " Unexpected method: ", trim( Method ) stop end select else IDMethod = IDMethodKD end if ! Check arguments ! select case ( IDMethod ) case ( IDMethodLBLConstPF ) if ( present( a_ArgBandWNBnds ) ) then else write( 6, * ) "In module, ", trim( ModuleName ) write( 6, * ) 'a_BandWNBnds has to be present when IDMethod is IDMethodLBLConstPF.' stop end if end select Grav = ArgGrav StrFluxTOA = ArgStrFluxTOA if ( present( ArgRayScatCoefNcFN ) ) then RayScatCoefNcFN = ArgRayScatCoefNcFN else RayScatCoefNcFN = "" end if if ( present( ArgStrSpeNcFN ) ) then StrSpeNcFN = ArgStrSpeNcFN else StrSpeNcFN = "" end if call OptPropInit( & & Method, & & OptPropNcFN, RayScatCoefNcFN, StrSpeNcFN & & ) call RadRTETwoStreamApp1DWrapperInit select case ( IDMethod ) case ( IDMethodKDDevelop, IDMethodKD ) NBand = OptPropGetNBand() allocate( aa_BandWNBnds (2,NBand) ) allocate( aa_BandBinIndexBnds(2,NBand) ) call OptPropGetBandWNBnds( & & NBand, & & aa_BandWNBnds & & ) call OptPropGetBandBinIndexBnds( & & NBand, & & aa_BandBinIndexBnds & & ) case ( IDMethodLBLConstPF ) NBand = size( a_ArgBandWNBnds ) - 1 allocate( aa_BandWNBnds(2,NBand) ) do n = 1, NBand aa_BandWNBnds(1,n) = a_ArgBandWNBnds(n ) aa_BandWNBnds(2,n) = a_ArgBandWNBnds(n+1) end do end select ! Check index for method switch from LW to SW iWaveNumMaxforMethodLW = SearchBoundaryIndex( WaveNumMethodSwitch ) ! Check index for albedo switch from LW to SW iWaveNumMaxforAlbedoLW = SearchBoundaryIndex( WaveNumAlbedoSwitch ) !!$ call OptPropGetTotalStrFlux( & !!$ & TotalStrFlux & !!$ & ) !!$ FluxRatTOAtoStelSfc = StrFluxTOA / TotalStrFlux FlagInited = .true. end subroutine RadFluxInit !---------------------------------------------------------------------------- subroutine RadFluxFinalize use opt_prop, only : OptPropFinalize ! local variables ! call OptPropFinalize if ( FlagOutputInited ) then call RadFluxOutputFinalize end if select case ( IDMethod ) case ( IDMethodKDDevelop ) deallocate( aa_BandWNBnds ) deallocate( aa_BandBinIndexBnds ) case ( IDMethodLBLConstPF ) deallocate( aa_BandWNBnds ) end select FlagInited = .false. end subroutine RadFluxFinalize !---------------------------------------------------------------------------- function SearchBoundaryIndex( WaveNumSwitch ) result( iWaveNum ) use opt_prop, only : & & OptPropGetNWN, & & OptPropGetWN real(DP), intent(in ) :: WaveNumSwitch integer :: iWaveNum ! Local variables ! integer :: NWaveNum integer :: iBand NWaveNum = OptPropGetNWN() ! Check index for switching from LW to SW select case ( IDMethod ) case ( IDMethodKDDevelop, IDMethodKD ) if ( WaveNumSwitch < aa_BandWNBnds(1,1) ) then ! ... for short wave is used in all bands. iWaveNum = 0 else if ( aa_BandWNBnds(2,NBand) < WaveNumSwitch ) then ! ... for long wave is used in all bands. iWaveNum = NWaveNum else do iBand = 1, NBand if ( ( aa_BandWNBnds(1,iBand) <= WaveNumSwitch ) .and. & & ( WaveNumSwitch <= aa_BandWNBnds(2,iBand) ) ) then exit end if end do if ( ( WaveNumSwitch - aa_BandWNBnds(1,iBand) ) & & < ( aa_BandWNBnds(2,iBand) - WaveNumSwitch ) ) then iBand = iBand - 1 end if ! obtain max bin number for ... for long wave if ( iBand >= 1 ) then iWaveNum = aa_BandBinIndexBnds(2,iBand) else iWaveNum = 0 end if end if case ( IDMethodLBL, IDMethodLBLConstPF ) if ( WaveNumSwitch < OptPropGetWN( 1 ) ) then ! ... for short wave is used in all bands. iWaveNum = 0 else if ( OptPropGetWN( NWaveNum ) < WaveNumSwitch ) then ! ... for long wave is used in all bands. iWaveNum = NWaveNum else do iWaveNum = 1, NWaveNum if ( WaveNumSwitch <= OptPropGetWN( iWaveNum ) ) exit end do iWaveNum = iWaveNum - 1 end if case default write( 6, * ) 'In ', trim( ModuleName ) write( 6, * ) ' Unexpected case default.' stop end select end function SearchBoundaryIndex !---------------------------------------------------------------------------- subroutine RadFluxOutputInit( & & OptPropNcFN, & & FluxOutNcFn, & & kmax, r_Press & & ) use ni3_module use netcdf character(*), intent(in) :: OptPropNcFN character(*), intent(in) :: FluxOutNcFn integer , intent(in) :: kmax real(8) , intent(in) :: r_Press (0:kmax) ! local variables ! real(8) :: z_Press (1:kmax) integer :: NcID character(stdstr) :: Mode character(stdstr) :: Name character(stdstr) :: StdName character(stdstr) :: LongName character(stdstr) :: Units character(stdstr) :: aa_DimNames(2) integer :: LBLNWN real(DP), allocatable :: w_WaveNum(:) integer :: k if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if do k = 1, kmax z_Press(k) = ( r_Press(k-1) + r_Press(k) ) / 2.0d0 end do Mode = "read" call ni3_open( OptPropNcFN, Mode, NcID ) call ni3_inq_dimlen( NcID, 'WaveNum', LBLNWN ) allocate( w_WaveNum( 1:LBLNWN ) ) call ni3_get_var( NcID, 'WaveNum', w_WaveNum ) call ni3_close( NcID ) ! Output Mode = "new" call ni3_open( FluxOutNcFN, Mode, NcIDFlux ) Name = "Press" StdName = "plev" Units = "Pa" call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, r_Press , & & stdname = stdname, units = units ) Name = "WaveNum" StdName = "wavenumber" Units = "m-1" call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, w_WaveNum , & & stdname = stdname, units = units ) Name = "PressLC" StdName = "plev" Units = "Pa" call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, z_Press , & & stdname = stdname, units = units ) aa_DimNames(1) = "Press" aa_DimNames(2) = "WaveNum" Name = "RadUwFlux" StdName = Name LongName = "radiative upward flux" Units = "W m-2 (m-1)-1" call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, & & stdname = StdName, longname = LongName, Units = Units ) Name = "RadDwFlux" StdName = Name LongName = "radiative downward flux" Units = "W m-2 (m-1)-1" call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, & & stdname = StdName, longname = LongName, Units = Units ) Name = "OptDep" StdName = Name LongName = "optical depth" Units = "1" call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, & & stdname = StdName, longname = LongName, Units = Units ) aa_DimNames(1) = "PressLC" aa_DimNames(2) = "WaveNum" Name = "DRadFluxDPress" StdName = Name LongName = "flux divergence" Units = "W m-2 Pa-1 (m-1)-1" call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, & & stdname = StdName, longname = LongName, Units = Units ) deallocate( w_WaveNum ) FlagOutputInited = .true. end subroutine RadFluxOutputInit !---------------------------------------------------------------------------- subroutine RadFluxOutputFinalize use ni3_module ! local variables ! if ( .not. FlagInited ) then write( 6, * ) trim( ModuleName ), " is not initialized." stop end if if ( FlagOutputInited ) then call ni3_close( NcIDFlux ) FlagOutputInited = .false. end if end subroutine RadFluxOutputFinalize !-------------------------------------------------------------------------- end module rad_flux