module SurfaceRadiationMod !------------------------------------------------------------------------------ !BOP ! ! !MODULE: SurfaceRadiationMod ! ! !DESCRIPTION: ! Calculate solar fluxes absorbed by vegetation and ground surface ! ! !USES: use shr_kind_mod, only: r8 => shr_kind_r8 use clm_varctl , only: iulog ! !PUBLIC TYPES: implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: SurfaceRadiation ! Solar fluxes absorbed by veg and ground surface ! ! !REVISION HISTORY: ! Created by Mariana Vertenstein ! 11/26/03, Peter Thornton: Added new routine for improved treatment of ! sunlit/shaded canopy radiation. ! 4/26/05, Peter Thornton: Adopted the sun/shade algorithm as the default, ! removed the old SurfaceRadiation(), and renamed SurfaceRadiationSunShade() ! as SurfaceRadiation(). ! !EOP !------------------------------------------------------------------------------ contains !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SurfaceRadiation ! ! !INTERFACE: subroutine SurfaceRadiation(lbp, ubp, num_nourbanp, filter_nourbanp) ! ! !DESCRIPTION: ! Solar fluxes absorbed by vegetation and ground surface ! Note possible problem when land is on different grid than atmosphere. ! Land may have sun above the horizon (coszen > 0) but atmosphere may ! have sun below the horizon (forc_solad = 0 and forc_solai = 0). This is okay ! because all fluxes (absorbed, reflected, transmitted) are multiplied ! by the incoming flux and all will equal zero. ! Atmosphere may have sun above horizon (forc_solad > 0 and forc_solai > 0) but ! land may have sun below horizon. This is okay because fabd, fabi, ! ftdd, ftid, and ftii all equal zero so that sabv=sabg=fsa=0. Also, ! albd and albi equal one so that fsr=forc_solad+forc_solai. In other words, all ! the radiation is reflected. NDVI should equal zero in this case. ! However, the way the code is currently implemented this is only true ! if (forc_solad+forc_solai)|vis = (forc_solad+forc_solai)|nir. ! Output variables are parsun,parsha,sabv,sabg,fsa,fsr,ndvi ! ! !USES: use clmtype use clm_atmlnd , only : clm_a2l use clm_varpar , only : numrad use clm_varcon , only : spval, istsoil, degpsec, isecspday use clm_varcon , only : istice_mec use clm_varcon , only : istcrop use clm_time_manager, only : get_curr_date, get_step_size use clm_varpar , only : nlevsno use SNICARMod , only : DO_SNO_OC use abortutils , only : endrun ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbp, ubp ! pft upper and lower bounds integer, intent(in) :: num_nourbanp ! number of pfts in non-urban points in pft filter integer, intent(in) :: filter_nourbanp(ubp-lbp+1) ! pft filter for non-urban points ! ! !CALLED FROM: ! subroutine Biogeophysics1 in module Biogeophysics1Mod ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 2/18/02, Peter Thornton: Migrated to new data structures. Added a pft loop. ! 6/05/03, Peter Thornton: Modified sunlit/shaded canopy treatment. Original code ! had all radiation being absorbed in the sunlit canopy, and now the sunlit and shaded ! canopies are each given the appropriate fluxes. There was also an inconsistency in ! the original code, where parsun was not being scaled by leaf area, and so represented ! the entire canopy flux. This goes into Stomata (in CanopyFluxes) where it is assumed ! to be a flux per unit leaf area. In addition, the fpsn flux coming out of Stomata was ! being scaled back up to the canopy by multiplying by lai, but the input radiation flux was ! for the entire canopy to begin with. Corrected this inconsistency in this version, so that ! the parsun and parsha fluxes going into canopy fluxes are per unit lai in the sunlit and ! shaded canopies. ! 6/9/03, Peter Thornton: Moved coszen from g%gps to c%cps to avoid problem ! with OpenMP threading over columns, where different columns hit the radiation ! time step at different times during execution. ! 6/10/03, Peter Thornton: Added constraint on negative tot_aid, instead of ! exiting with error. Appears to be happening only at roundoff level. ! 6/11/03, Peter Thornton: Moved calculation of ext inside if (coszen), ! and added check on laisun = 0 and laisha = 0 in calculation of sun_aperlai ! and sha_aperlai. ! 11/26/03, Peter Thornton: During migration to new vector code, created ! this as a new routine to handle sunlit/shaded canopy calculations. ! 03/28/08, Mark Flanner: Incorporated SNICAR, including absorbed solar radiation ! in each snow layer and top soil layer, and optional radiative forcing calculation ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments ! integer , pointer :: ivt(:) ! pft vegetation type integer , pointer :: pcolumn(:) ! pft's column index integer , pointer :: pgridcell(:) ! pft's gridcell index real(r8), pointer :: pwtgcell(:) ! pft's weight relative to corresponding gridcell real(r8), pointer :: elai(:) ! one-sided leaf area index with burying by snow real(r8), pointer :: esai(:) ! one-sided stem area index with burying by snow real(r8), pointer :: londeg(:) ! longitude (degrees) real(r8), pointer :: latdeg(:) ! latitude (degrees) real(r8), pointer :: slasun(:) ! specific leaf area for sunlit canopy, projected area basis (m^2/gC) real(r8), pointer :: slasha(:) ! specific leaf area for shaded canopy, projected area basis (m^2/gC) real(r8), pointer :: gdir(:) ! leaf projection in solar direction (0 to 1) real(r8), pointer :: omega(:,:) ! fraction of intercepted radiation that is scattered (0 to 1) real(r8), pointer :: coszen(:) ! cosine of solar zenith angle real(r8), pointer :: forc_solad(:,:) ! direct beam radiation (W/m**2) real(r8), pointer :: forc_solai(:,:) ! diffuse radiation (W/m**2) real(r8), pointer :: fabd(:,:) ! flux absorbed by veg per unit direct flux real(r8), pointer :: fabi(:,:) ! flux absorbed by veg per unit diffuse flux real(r8), pointer :: ftdd(:,:) ! down direct flux below veg per unit dir flx real(r8), pointer :: ftid(:,:) ! down diffuse flux below veg per unit dir flx real(r8), pointer :: ftii(:,:) ! down diffuse flux below veg per unit dif flx real(r8), pointer :: albgrd(:,:) ! ground albedo (direct) real(r8), pointer :: albgri(:,:) ! ground albedo (diffuse) real(r8), pointer :: albd(:,:) ! surface albedo (direct) real(r8), pointer :: albi(:,:) ! surface albedo (diffuse) real(r8), pointer :: slatop(:) ! specific leaf area at top of canopy, projected area basis [m^2/gC] real(r8), pointer :: dsladlai(:) ! dSLA/dLAI, projected area basis [m^2/gC] ! ! local pointers to original implicit out arguments ! real(r8), pointer :: fsun(:) ! sunlit fraction of canopy real(r8), pointer :: laisun(:) ! sunlit leaf area real(r8), pointer :: laisha(:) ! shaded leaf area real(r8), pointer :: sabg(:) ! solar radiation absorbed by ground (W/m**2) real(r8), pointer :: sabv(:) ! solar radiation absorbed by vegetation (W/m**2) real(r8), pointer :: fsa(:) ! solar radiation absorbed (total) (W/m**2) real(r8), pointer :: fsa_r(:) ! rural solar radiation absorbed (total) (W/m**2) integer , pointer :: ityplun(:) ! landunit type integer , pointer :: plandunit(:) ! index into landunit level quantities real(r8), pointer :: parsun(:) ! average absorbed PAR for sunlit leaves (W/m**2) real(r8), pointer :: parsha(:) ! average absorbed PAR for shaded leaves (W/m**2) real(r8), pointer :: fsr(:) ! solar radiation reflected (W/m**2) real(r8), pointer :: fsds_vis_d(:) ! incident direct beam vis solar radiation (W/m**2) real(r8), pointer :: fsds_nir_d(:) ! incident direct beam nir solar radiation (W/m**2) real(r8), pointer :: fsds_vis_i(:) ! incident diffuse vis solar radiation (W/m**2) real(r8), pointer :: fsds_nir_i(:) ! incident diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsr_vis_d(:) ! reflected direct beam vis solar radiation (W/m**2) real(r8), pointer :: fsr_nir_d(:) ! reflected direct beam nir solar radiation (W/m**2) real(r8), pointer :: fsr_vis_i(:) ! reflected diffuse vis solar radiation (W/m**2) real(r8), pointer :: fsr_nir_i(:) ! reflected diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsds_vis_d_ln(:) ! incident direct beam vis solar rad at local noon (W/m**2) real(r8), pointer :: fsds_nir_d_ln(:) ! incident direct beam nir solar rad at local noon (W/m**2) real(r8), pointer :: fsr_vis_d_ln(:) ! reflected direct beam vis solar rad at local noon (W/m**2) real(r8), pointer :: fsr_nir_d_ln(:) ! reflected direct beam nir solar rad at local noon (W/m**2) real(r8), pointer :: eff_kid(:,:) ! effective extinction coefficient for indirect from direct real(r8), pointer :: eff_kii(:,:) ! effective extinction coefficient for indirect from indirect real(r8), pointer :: sun_faid(:,:) ! fraction sun canopy absorbed indirect from direct real(r8), pointer :: sun_faii(:,:) ! fraction sun canopy absorbed indirect from indirect real(r8), pointer :: sha_faid(:,:) ! fraction shade canopy absorbed indirect from direct real(r8), pointer :: sha_faii(:,:) ! fraction shade canopy absorbed indirect from indirect real(r8), pointer :: sun_add(:,:) ! sun canopy absorbed direct from direct (W/m**2) real(r8), pointer :: tot_aid(:,:) ! total canopy absorbed indirect from direct (W/m**2) real(r8), pointer :: sun_aid(:,:) ! sun canopy absorbed indirect from direct (W/m**2) real(r8), pointer :: sun_aii(:,:) ! sun canopy absorbed indirect from indirect (W/m**2) real(r8), pointer :: sha_aid(:,:) ! shade canopy absorbed indirect from direct (W/m**2) real(r8), pointer :: sha_aii(:,:) ! shade canopy absorbed indirect from indirect (W/m**2) real(r8), pointer :: sun_atot(:,:) ! sun canopy total absorbed (W/m**2) real(r8), pointer :: sha_atot(:,:) ! shade canopy total absorbed (W/m**2) real(r8), pointer :: sun_alf(:,:) ! sun canopy total absorbed by leaves (W/m**2) real(r8), pointer :: sha_alf(:,:) ! shade canopy total absored by leaves (W/m**2) real(r8), pointer :: sun_aperlai(:,:) ! sun canopy total absorbed per unit LAI (W/m**2) real(r8), pointer :: sha_aperlai(:,:) ! shade canopy total absorbed per unit LAI (W/m**2) real(r8), pointer :: flx_absdv(:,:) ! direct flux absorption factor (col,lyr): VIS [frc] real(r8), pointer :: flx_absdn(:,:) ! direct flux absorption factor (col,lyr): NIR [frc] real(r8), pointer :: flx_absiv(:,:) ! diffuse flux absorption factor (col,lyr): VIS [frc] real(r8), pointer :: flx_absin(:,:) ! diffuse flux absorption factor (col,lyr): NIR [frc] integer , pointer :: snl(:) ! negative number of snow layers [nbr] real(r8), pointer :: albgrd_pur(:,:) ! pure snow ground albedo (direct) real(r8), pointer :: albgri_pur(:,:) ! pure snow ground albedo (diffuse) real(r8), pointer :: albgrd_bc(:,:) ! ground albedo without BC (direct) (col,bnd) real(r8), pointer :: albgri_bc(:,:) ! ground albedo without BC (diffuse) (col,bnd) real(r8), pointer :: albgrd_oc(:,:) ! ground albedo without OC (direct) (col,bnd) real(r8), pointer :: albgri_oc(:,:) ! ground albedo without OC (diffuse) (col,bnd) real(r8), pointer :: albgrd_dst(:,:) ! ground albedo without dust (direct) (col,bnd) real(r8), pointer :: albgri_dst(:,:) ! ground albedo without dust (diffuse) (col,bnd) real(r8), pointer :: albsnd_hst(:,:) ! snow albedo, direct, for history files (col,bnd) [frc] real(r8), pointer :: albsni_hst(:,:) ! snow ground albedo, diffuse, for history files (col,bnd real(r8), pointer :: sabg_lyr(:,:) ! absorbed radiative flux (pft,lyr) [W/m2] real(r8), pointer :: sfc_frc_aer(:) ! surface forcing of snow with all aerosols (pft) [W/m2] real(r8), pointer :: sfc_frc_bc(:) ! surface forcing of snow with BC (pft) [W/m2] real(r8), pointer :: sfc_frc_oc(:) ! surface forcing of snow with OC (pft) [W/m2] real(r8), pointer :: sfc_frc_dst(:) ! surface forcing of snow with dust (pft) [W/m2] real(r8), pointer :: sfc_frc_aer_sno(:) ! surface forcing of snow with all aerosols, averaged only when snow is present (pft) [W/m2] real(r8), pointer :: sfc_frc_bc_sno(:) ! surface forcing of snow with BC, averaged only when snow is present (pft) [W/m2] real(r8), pointer :: sfc_frc_oc_sno(:) ! surface forcing of snow with OC, averaged only when snow is present (pft) [W/m2] real(r8), pointer :: sfc_frc_dst_sno(:) ! surface forcing of snow with dust, averaged only when snow is present (pft) [W/m2] real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: fsr_sno_vd(:) ! reflected visible, direct radiation from snow (for history files) (pft) [W/m2] real(r8), pointer :: fsr_sno_nd(:) ! reflected near-IR, direct radiation from snow (for history files) (pft) [W/m2] real(r8), pointer :: fsr_sno_vi(:) ! reflected visible, diffuse radiation from snow (for history files) (pft) [W/m2] real(r8), pointer :: fsr_sno_ni(:) ! reflected near-IR, diffuse radiation from snow (for history files) (pft) [W/m2] real(r8), pointer :: fsds_sno_vd(:) ! incident visible, direct radiation on snow (for history files) (pft) [W/m2] real(r8), pointer :: fsds_sno_nd(:) ! incident near-IR, direct radiation on snow (for history files) (pft) [W/m2] real(r8), pointer :: fsds_sno_vi(:) ! incident visible, diffuse radiation on snow (for history files) (pft) [W/m2] real(r8), pointer :: fsds_sno_ni(:) ! incident near-IR, diffuse radiation on snow (for history files) (pft) [W/m2] real(r8), pointer :: snowdp(:) ! snow height (m) ! ! ! !OTHER LOCAL VARIABLES: !EOP ! integer , parameter :: nband = numrad ! number of solar radiation waveband classes real(r8), parameter :: mpe = 1.e-06_r8 ! prevents overflow for division by zero integer :: fp ! non-urban filter pft index integer :: p ! pft index integer :: c ! column index integer :: l ! landunit index integer :: g ! grid cell index integer :: ib ! waveband number (1=vis, 2=nir) real(r8) :: absrad ! absorbed solar radiation (W/m**2) real(r8) :: rnir ! reflected solar radiation [nir] (W/m**2) real(r8) :: rvis ! reflected solar radiation [vis] (W/m**2) real(r8) :: laifra ! leaf area fraction of canopy real(r8) :: trd(lbp:ubp,numrad) ! transmitted solar radiation: direct (W/m**2) real(r8) :: tri(lbp:ubp,numrad) ! transmitted solar radiation: diffuse (W/m**2) real(r8) :: cad(lbp:ubp,numrad) ! direct beam absorbed by canopy (W/m**2) real(r8) :: cai(lbp:ubp,numrad) ! diffuse radiation absorbed by canopy (W/m**2) real(r8) :: vai(lbp:ubp) ! total leaf area index + stem area index, one sided real(r8) :: ext ! optical depth direct beam per unit LAI+SAI real(r8) :: t1, t2 ! temporary variables real(r8) :: cosz integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) integer :: year,month,day,secs ! calendar info for current time step integer :: i ! layer index [idx] real(r8) :: sabg_snl_sum ! temporary, absorbed energy in all active snow layers [W/m2] real(r8) :: absrad_pur ! temp: absorbed solar radiation by pure snow [W/m2] real(r8) :: absrad_bc ! temp: absorbed solar radiation without BC [W/m2] real(r8) :: absrad_oc ! temp: absorbed solar radiation without OC [W/m2] real(r8) :: absrad_dst ! temp: absorbed solar radiation without dust [W/m2] real(r8) :: sabg_pur(lbp:ubp) ! solar radiation absorbed by ground with pure snow [W/m2] real(r8) :: sabg_bc(lbp:ubp) ! solar radiation absorbed by ground without BC [W/m2] real(r8) :: sabg_oc(lbp:ubp) ! solar radiation absorbed by ground without OC [W/m2] real(r8) :: sabg_dst(lbp:ubp) ! solar radiation absorbed by ground without dust [W/m2] !------------------------------------------------------------------------------ ! Assign local pointers to multi-level derived type members (gridcell level) londeg => clm3%g%londeg latdeg => clm3%g%latdeg forc_solad => clm_a2l%forc_solad forc_solai => clm_a2l%forc_solai ! Assign local pointers to multi-level derived type members (landunit level) ityplun => clm3%g%l%itype ! Assign local pointers to multi-level derived type members (column level) albgrd => clm3%g%l%c%cps%albgrd albgri => clm3%g%l%c%cps%albgri coszen => clm3%g%l%c%cps%coszen ! Assign local pointers to derived type members (pft-level) plandunit => clm3%g%l%c%p%landunit ivt => clm3%g%l%c%p%itype pcolumn => clm3%g%l%c%p%column pgridcell => clm3%g%l%c%p%gridcell pwtgcell => clm3%g%l%c%p%wtgcell elai => clm3%g%l%c%p%pps%elai esai => clm3%g%l%c%p%pps%esai slasun => clm3%g%l%c%p%pps%slasun slasha => clm3%g%l%c%p%pps%slasha gdir => clm3%g%l%c%p%pps%gdir omega => clm3%g%l%c%p%pps%omega laisun => clm3%g%l%c%p%pps%laisun laisha => clm3%g%l%c%p%pps%laisha fabd => clm3%g%l%c%p%pps%fabd fabi => clm3%g%l%c%p%pps%fabi ftdd => clm3%g%l%c%p%pps%ftdd ftid => clm3%g%l%c%p%pps%ftid ftii => clm3%g%l%c%p%pps%ftii albd => clm3%g%l%c%p%pps%albd albi => clm3%g%l%c%p%pps%albi fsun => clm3%g%l%c%p%pps%fsun sabg => clm3%g%l%c%p%pef%sabg sabv => clm3%g%l%c%p%pef%sabv snowdp => clm3%g%l%c%cps%snowdp fsa => clm3%g%l%c%p%pef%fsa fsa_r => clm3%g%l%c%p%pef%fsa_r fsr => clm3%g%l%c%p%pef%fsr parsun => clm3%g%l%c%p%pef%parsun parsha => clm3%g%l%c%p%pef%parsha fsds_vis_d => clm3%g%l%c%p%pef%fsds_vis_d fsds_nir_d => clm3%g%l%c%p%pef%fsds_nir_d fsds_vis_i => clm3%g%l%c%p%pef%fsds_vis_i fsds_nir_i => clm3%g%l%c%p%pef%fsds_nir_i fsr_vis_d => clm3%g%l%c%p%pef%fsr_vis_d fsr_nir_d => clm3%g%l%c%p%pef%fsr_nir_d fsr_vis_i => clm3%g%l%c%p%pef%fsr_vis_i fsr_nir_i => clm3%g%l%c%p%pef%fsr_nir_i fsds_vis_d_ln => clm3%g%l%c%p%pef%fsds_vis_d_ln fsds_nir_d_ln => clm3%g%l%c%p%pef%fsds_nir_d_ln fsr_vis_d_ln => clm3%g%l%c%p%pef%fsr_vis_d_ln fsr_nir_d_ln => clm3%g%l%c%p%pef%fsr_nir_d_ln eff_kid => clm3%g%l%c%p%pps%eff_kid eff_kii => clm3%g%l%c%p%pps%eff_kii sun_faid => clm3%g%l%c%p%pps%sun_faid sun_faii => clm3%g%l%c%p%pps%sun_faii sha_faid => clm3%g%l%c%p%pps%sha_faid sha_faii => clm3%g%l%c%p%pps%sha_faii sun_add => clm3%g%l%c%p%pef%sun_add tot_aid => clm3%g%l%c%p%pef%tot_aid sun_aid => clm3%g%l%c%p%pef%sun_aid sun_aii => clm3%g%l%c%p%pef%sun_aii sha_aid => clm3%g%l%c%p%pef%sha_aid sha_aii => clm3%g%l%c%p%pef%sha_aii sun_atot => clm3%g%l%c%p%pef%sun_atot sha_atot => clm3%g%l%c%p%pef%sha_atot sun_alf => clm3%g%l%c%p%pef%sun_alf sha_alf => clm3%g%l%c%p%pef%sha_alf sun_aperlai => clm3%g%l%c%p%pef%sun_aperlai sha_aperlai => clm3%g%l%c%p%pef%sha_aperlai ! Assign local pointers to derived type members (ecophysiological) slatop => pftcon%slatop dsladlai => pftcon%dsladlai frac_sno => clm3%g%l%c%cps%frac_sno flx_absdv => clm3%g%l%c%cps%flx_absdv flx_absdn => clm3%g%l%c%cps%flx_absdn flx_absiv => clm3%g%l%c%cps%flx_absiv flx_absin => clm3%g%l%c%cps%flx_absin sabg_lyr => clm3%g%l%c%p%pef%sabg_lyr snl => clm3%g%l%c%cps%snl sfc_frc_aer => clm3%g%l%c%p%pef%sfc_frc_aer sfc_frc_aer_sno => clm3%g%l%c%p%pef%sfc_frc_aer_sno albgrd_pur => clm3%g%l%c%cps%albgrd_pur albgri_pur => clm3%g%l%c%cps%albgri_pur sfc_frc_bc => clm3%g%l%c%p%pef%sfc_frc_bc sfc_frc_bc_sno => clm3%g%l%c%p%pef%sfc_frc_bc_sno albgrd_bc => clm3%g%l%c%cps%albgrd_bc albgri_bc => clm3%g%l%c%cps%albgri_bc sfc_frc_oc => clm3%g%l%c%p%pef%sfc_frc_oc sfc_frc_oc_sno => clm3%g%l%c%p%pef%sfc_frc_oc_sno albgrd_oc => clm3%g%l%c%cps%albgrd_oc albgri_oc => clm3%g%l%c%cps%albgri_oc sfc_frc_dst => clm3%g%l%c%p%pef%sfc_frc_dst sfc_frc_dst_sno => clm3%g%l%c%p%pef%sfc_frc_dst_sno albgrd_dst => clm3%g%l%c%cps%albgrd_dst albgri_dst => clm3%g%l%c%cps%albgri_dst albsnd_hst => clm3%g%l%c%cps%albsnd_hst albsni_hst => clm3%g%l%c%cps%albsni_hst fsr_sno_vd => clm3%g%l%c%p%pef%fsr_sno_vd fsr_sno_nd => clm3%g%l%c%p%pef%fsr_sno_nd fsr_sno_vi => clm3%g%l%c%p%pef%fsr_sno_vi fsr_sno_ni => clm3%g%l%c%p%pef%fsr_sno_ni fsds_sno_vd => clm3%g%l%c%p%pef%fsds_sno_vd fsds_sno_nd => clm3%g%l%c%p%pef%fsds_sno_nd fsds_sno_vi => clm3%g%l%c%p%pef%fsds_sno_vi fsds_sno_ni => clm3%g%l%c%p%pef%fsds_sno_ni ! Determine seconds off current time step dtime = get_step_size() call get_curr_date (year, month, day, secs) ! Determine fluxes do fp = 1,num_nourbanp p = filter_nourbanp(fp) ! was redundant b/c filter already included wt>0; ! not redundant anymore with chg in filter definition l = plandunit(p) !Note: Some glacier_mec pfts may have zero weight if (pwtgcell(p)>0._r8 .or. ityplun(l)==istice_mec) then sabg(p) = 0._r8 sabv(p) = 0._r8 fsa(p) = 0._r8 l = plandunit(p) if (ityplun(l)==istsoil .or. ityplun(l)==istcrop) then fsa_r(p) = 0._r8 end if sabg_lyr(p,:) = 0._r8 sabg_pur(p) = 0._r8 sabg_bc(p) = 0._r8 sabg_oc(p) = 0._r8 sabg_dst(p) = 0._r8 end if end do ! Loop over pfts to calculate fsun, etc do fp = 1,num_nourbanp p = filter_nourbanp(fp) l = plandunit(p) if (pwtgcell(p)>0._r8 .or. ityplun(l)==istice_mec) then c = pcolumn(p) g = pgridcell(p) vai(p) = elai(p) + esai(p) if (coszen(c) > 0._r8 .and. elai(p) > 0._r8 .and. gdir(p) > 0._r8) then cosz = max(0.001_r8, coszen(c)) ext = gdir(p)/cosz t1 = min(ext*elai(p), 40.0_r8) t2 = exp(-t1) fsun(p) = (1._r8-t2)/t1 ! new control on low lai, to avoid numerical problems in ! calculation of slasun, slasha ! PET: 2/29/04 if (elai(p) > 0.01_r8) then laisun(p) = elai(p)*fsun(p) laisha(p) = elai(p)*(1._r8-fsun(p)) ! calculate the average specific leaf area for sunlit and shaded ! canopies, when effective LAI > 0 slasun(p) = (t2*dsladlai(ivt(p))*ext*elai(p) + & t2*dsladlai(ivt(p)) + & t2*slatop(ivt(p))*ext - & dsladlai(ivt(p)) - & slatop(ivt(p))*ext) / & (ext*(t2-1._r8)) slasha(p) = ((slatop(ivt(p)) + & (dsladlai(ivt(p)) * elai(p)/2.0_r8)) * elai(p) - & laisun(p)*slasun(p)) / laisha(p) else ! special case for low elai fsun(p) = 1._r8 laisun(p) = elai(p) laisha(p) = 0._r8 slasun(p) = slatop(ivt(p)) slasha(p) = 0._r8 end if else fsun(p) = 0._r8 laisun(p) = 0._r8 laisha(p) = elai(p) slasun(p) = 0._r8 slasha(p) = 0._r8 end if end if end do ! Loop over nband wavebands do ib = 1, nband do fp = 1,num_nourbanp p = filter_nourbanp(fp) l = plandunit(p) if (pwtgcell(p)>0._r8 .or. ityplun(l)==istice_mec) then c = pcolumn(p) g = pgridcell(p) ! Absorbed by canopy cad(p,ib) = forc_solad(g,ib)*fabd(p,ib) cai(p,ib) = forc_solai(g,ib)*fabi(p,ib) sabv(p) = sabv(p) + cad(p,ib) + cai(p,ib) fsa(p) = fsa(p) + cad(p,ib) + cai(p,ib) l = plandunit(p) if (ityplun(l)==istsoil .or. ityplun(l)==istcrop) then fsa_r(p) = fsa_r(p) + cad(p,ib) + cai(p,ib) end if ! Transmitted = solar fluxes incident on ground trd(p,ib) = forc_solad(g,ib)*ftdd(p,ib) tri(p,ib) = forc_solad(g,ib)*ftid(p,ib) + forc_solai(g,ib)*ftii(p,ib) ! Solar radiation absorbed by ground surface absrad = trd(p,ib)*(1._r8-albgrd(c,ib)) + tri(p,ib)*(1._r8-albgri(c,ib)) sabg(p) = sabg(p) + absrad fsa(p) = fsa(p) + absrad if (ityplun(l)==istsoil .or. ityplun(l)==istcrop) then fsa_r(p) = fsa_r(p) + absrad end if #if (defined SNICAR_FRC) ! Solar radiation absorbed by ground surface without BC absrad_bc = trd(p,ib)*(1._r8-albgrd_bc(c,ib)) + tri(p,ib)*(1._r8-albgri_bc(c,ib)) sabg_bc(p) = sabg_bc(p) + absrad_bc ! Solar radiation absorbed by ground surface without OC absrad_oc = trd(p,ib)*(1._r8-albgrd_oc(c,ib)) + tri(p,ib)*(1._r8-albgri_oc(c,ib)) sabg_oc(p) = sabg_oc(p) + absrad_oc ! Solar radiation absorbed by ground surface without dust absrad_dst = trd(p,ib)*(1._r8-albgrd_dst(c,ib)) + tri(p,ib)*(1._r8-albgri_dst(c,ib)) sabg_dst(p) = sabg_dst(p) + absrad_dst ! Solar radiation absorbed by ground surface without any aerosols absrad_pur = trd(p,ib)*(1._r8-albgrd_pur(c,ib)) + tri(p,ib)*(1._r8-albgri_pur(c,ib)) sabg_pur(p) = sabg_pur(p) + absrad_pur #endif ! New sunlit.shaded canopy algorithm if (coszen(c) > 0._r8 .and. elai(p) > 0._r8 .and. gdir(p) > 0._r8 ) then ! 1. calculate flux of direct beam radiation absorbed in the ! sunlit canopy as direct (sun_add), and the flux of direct ! beam radiation absorbed in the total canopy as indirect sun_add(p,ib) = forc_solad(g,ib) * (1._r8-ftdd(p,ib)) * (1._r8-omega(p,ib)) tot_aid(p,ib) = (forc_solad(g,ib) * fabd(p,ib)) - sun_add(p,ib) ! the following constraint set to catch round-off level errors ! that can cause negative tot_aid tot_aid(p,ib) = max(tot_aid(p,ib), 0._r8) ! 2. calculate the effective extinction coefficients for indirect ! transmission originating from direct and indirect streams, ! using ftid and ftii !eff_kid(p,ib) = -(log(ftid(p,ib)))/vai(p) !eff_kii(p,ib) = -(log(ftii(p,ib)))/vai(p) ! 3. calculate the fraction of indirect radiation being absorbed ! in the sunlit and shaded canopy fraction. Some of this indirect originates in ! the direct beam and some originates in the indirect beam. !sun_faid(p,ib) = 1.-exp(-eff_kid(p,ib) * vaisun(p)) !sun_faii(p,ib) = 1.-exp(-eff_kii(p,ib) * vaisun(p)) sun_faid(p,ib) = fsun(p) sun_faii(p,ib) = fsun(p) sha_faid(p,ib) = 1._r8-sun_faid(p,ib) sha_faii(p,ib) = 1._r8-sun_faii(p,ib) ! 4. calculate the total indirect flux absorbed by the sunlit ! and shaded canopy based on these fractions and the fabd and ! fabi from surface albedo calculations sun_aid(p,ib) = tot_aid(p,ib) * sun_faid(p,ib) sun_aii(p,ib) = forc_solai(g,ib)*fabi(p,ib)*sun_faii(p,ib) sha_aid(p,ib) = tot_aid(p,ib) * sha_faid(p,ib) sha_aii(p,ib) = forc_solai(g,ib)*fabi(p,ib)*sha_faii(p,ib) ! 5. calculate the total flux absorbed in the sunlit and shaded ! canopy as the sum of these terms sun_atot(p,ib) = sun_add(p,ib) + sun_aid(p,ib) + sun_aii(p,ib) sha_atot(p,ib) = sha_aid(p,ib) + sha_aii(p,ib) ! 6. calculate the total flux absorbed by leaves in the sunlit ! and shaded canopies laifra = elai(p)/vai(p) sun_alf(p,ib) = sun_atot(p,ib) * laifra sha_alf(p,ib) = sha_atot(p,ib) * laifra ! 7. calculate the fluxes per unit lai in the sunlit and shaded ! canopies if (laisun(p) > 0._r8) then sun_aperlai(p,ib) = sun_alf(p,ib)/laisun(p) else sun_aperlai(p,ib) = 0._r8 endif if (laisha(p) > 0._r8) then sha_aperlai(p,ib) = sha_alf(p,ib)/laisha(p) else sha_aperlai(p,ib) = 0._r8 endif else ! coszen = 0 or elai = 0 sun_add(p,ib) = 0._r8 tot_aid(p,ib) = 0._r8 eff_kid(p,ib) = 0._r8 eff_kii(p,ib) = 0._r8 sun_faid(p,ib) = 0._r8 sun_faii(p,ib) = 0._r8 sha_faid(p,ib) = 0._r8 sha_faii(p,ib) = 0._r8 sun_aid(p,ib) = 0._r8 sun_aii(p,ib) = 0._r8 sha_aid(p,ib) = 0._r8 sha_aii(p,ib) = 0._r8 sun_atot(p,ib) = 0._r8 sha_atot(p,ib) = 0._r8 sun_alf(p,ib) = 0._r8 sha_alf(p,ib) = 0._r8 sun_aperlai(p,ib) = 0._r8 sha_aperlai(p,ib) = 0._r8 end if end if end do ! end of pft loop end do ! end nbands loop ! compute absorbed flux in each snow layer and top soil layer, ! based on flux factors computed in the radiative transfer portion of SNICAR. do fp = 1,num_nourbanp p = filter_nourbanp(fp) l = plandunit(p) if (pwtgcell(p)>0._r8 .or. ityplun(l)==istice_mec) then c = pcolumn(p) sabg_snl_sum = 0._r8 ! CASE1: No snow layers: all energy is absorbed in top soil layer if (snl(c) == 0) then sabg_lyr(p,:) = 0._r8 sabg_lyr(p,1) = sabg(p) sabg_snl_sum = sabg_lyr(p,1) ! CASE 2: Snow layers present: absorbed radiation is scaled according to ! flux factors computed by SNICAR else do i = -nlevsno+1,1,1 sabg_lyr(p,i) = flx_absdv(c,i)*trd(p,1) + flx_absdn(c,i)*trd(p,2) + & flx_absiv(c,i)*tri(p,1) + flx_absin(c,i)*tri(p,2) ! summed radiation in active snow layers: if (i >= snl(c)+1) then sabg_snl_sum = sabg_snl_sum + sabg_lyr(p,i) endif enddo ! Error handling: The situation below can occur when solar radiation is ! NOT computed every timestep. ! When the number of snow layers has changed in between computations of the ! absorbed solar energy in each layer, we must redistribute the absorbed energy ! to avoid physically unrealistic conditions. The assumptions made below are ! somewhat arbitrary, but this situation does not arise very frequently. ! This error handling is implemented to accomodate any value of the ! radiation frequency. if (abs(sabg_snl_sum-sabg(p)) > 0.00001_r8) then if (snl(c) == 0) then sabg_lyr(p,-4:0) = 0._r8 sabg_lyr(p,1) = sabg(p) elseif (snl(c) == -1) then sabg_lyr(p,-4:-1) = 0._r8 sabg_lyr(p,0) = sabg(p)*0.6_r8 sabg_lyr(p,1) = sabg(p)*0.4_r8 else sabg_lyr(p,:) = 0._r8 sabg_lyr(p,snl(c)+1) = sabg(p)*0.75_r8 sabg_lyr(p,snl(c)+2) = sabg(p)*0.25_r8 endif endif ! If shallow snow depth, all solar radiation absorbed in top or top two snow layers ! to prevent unrealistic timestep soil warming if (snowdp(c) < 0.10_r8) then if (snl(c) == 0) then sabg_lyr(p,-4:0) = 0._r8 sabg_lyr(p,1) = sabg(p) elseif (snl(c) == -1) then sabg_lyr(p,-4:-1) = 0._r8 sabg_lyr(p,0) = sabg(p) sabg_lyr(p,1) = 0._r8 else sabg_lyr(p,:) = 0._r8 sabg_lyr(p,snl(c)+1) = sabg(p)*0.75_r8 sabg_lyr(p,snl(c)+2) = sabg(p)*0.25_r8 endif endif endif ! This situation should not happen: if (abs(sum(sabg_lyr(p,:))-sabg(p)) > 0.00001_r8) then write(iulog,*) "SNICAR ERROR: Absorbed ground radiation not equal to summed snow layer radiation. pft = ", & p," Col= ", c, " Diff= ",sum(sabg_lyr(p,:))-sabg(p), " sabg(p)= ", sabg(p), " sabg_sum(p)= ", & sum(sabg_lyr(p,:)), " snl(c)= ", snl(c) write(iulog,*) "flx_absdv1= ", trd(p,1)*(1.-albgrd(c,1)), "flx_absdv2= ", sum(flx_absdv(c,:))*trd(p,1) write(iulog,*) "flx_absiv1= ", tri(p,1)*(1.-albgri(c,1))," flx_absiv2= ", sum(flx_absiv(c,:))*tri(p,1) write(iulog,*) "flx_absdn1= ", trd(p,2)*(1.-albgrd(c,2))," flx_absdn2= ", sum(flx_absdn(c,:))*trd(p,2) write(iulog,*) "flx_absin1= ", tri(p,2)*(1.-albgri(c,2))," flx_absin2= ", sum(flx_absin(c,:))*tri(p,2) write(iulog,*) "albgrd_nir= ", albgrd(c,2) write(iulog,*) "coszen= ", coszen(c) call endrun() endif #if (defined SNICAR_FRC) ! BC aerosol forcing (pft-level): sfc_frc_bc(p) = sabg(p) - sabg_bc(p) ! OC aerosol forcing (pft-level): if (DO_SNO_OC) then sfc_frc_oc(p) = sabg(p) - sabg_oc(p) else sfc_frc_oc(p) = 0._r8 endif ! dust aerosol forcing (pft-level): sfc_frc_dst(p) = sabg(p) - sabg_dst(p) ! all-aerosol forcing (pft-level): sfc_frc_aer(p) = sabg(p) - sabg_pur(p) ! forcings averaged only over snow: if (frac_sno(c) > 0._r8) then sfc_frc_bc_sno(p) = sfc_frc_bc(p)/frac_sno(c) sfc_frc_oc_sno(p) = sfc_frc_oc(p)/frac_sno(c) sfc_frc_dst_sno(p) = sfc_frc_dst(p)/frac_sno(c) sfc_frc_aer_sno(p) = sfc_frc_aer(p)/frac_sno(c) else sfc_frc_bc_sno(p) = spval sfc_frc_oc_sno(p) = spval sfc_frc_dst_sno(p) = spval sfc_frc_aer_sno(p) = spval endif #endif endif enddo do fp = 1,num_nourbanp p = filter_nourbanp(fp) l = plandunit(p) if (pwtgcell(p)>0._r8 .or. ityplun(l)==istice_mec) then g = pgridcell(p) ! Final step of new sunlit/shaded canopy algorithm ! 8. calculate the total and per-unit-lai fluxes for PAR in the ! sunlit and shaded canopy leaf fractions parsun(p) = sun_aperlai(p,1) parsha(p) = sha_aperlai(p,1) ! The following code is duplicated from SurfaceRadiation ! NDVI and reflected solar radiation rvis = albd(p,1)*forc_solad(g,1) + albi(p,1)*forc_solai(g,1) rnir = albd(p,2)*forc_solad(g,2) + albi(p,2)*forc_solai(g,2) fsr(p) = rvis + rnir fsds_vis_d(p) = forc_solad(g,1) fsds_nir_d(p) = forc_solad(g,2) fsds_vis_i(p) = forc_solai(g,1) fsds_nir_i(p) = forc_solai(g,2) fsr_vis_d(p) = albd(p,1)*forc_solad(g,1) fsr_nir_d(p) = albd(p,2)*forc_solad(g,2) fsr_vis_i(p) = albi(p,1)*forc_solai(g,1) fsr_nir_i(p) = albi(p,2)*forc_solai(g,2) local_secp1 = secs + nint((londeg(g)/degpsec)/dtime)*dtime local_secp1 = mod(local_secp1,isecspday) if (local_secp1 == isecspday/2) then fsds_vis_d_ln(p) = forc_solad(g,1) fsds_nir_d_ln(p) = forc_solad(g,2) fsr_vis_d_ln(p) = albd(p,1)*forc_solad(g,1) fsr_nir_d_ln(p) = albd(p,2)*forc_solad(g,2) else fsds_vis_d_ln(p) = spval fsds_nir_d_ln(p) = spval fsr_vis_d_ln(p) = spval fsr_nir_d_ln(p) = spval end if ! diagnostic variables (downwelling and absorbed radiation partitioning) for history files ! (OPTIONAL) c = pcolumn(p) if (snl(c) < 0) then fsds_sno_vd(p) = forc_solad(g,1) fsds_sno_nd(p) = forc_solad(g,2) fsds_sno_vi(p) = forc_solai(g,1) fsds_sno_ni(p) = forc_solai(g,2) fsr_sno_vd(p) = fsds_vis_d(p)*albsnd_hst(c,1) fsr_sno_nd(p) = fsds_nir_d(p)*albsnd_hst(c,2) fsr_sno_vi(p) = fsds_vis_i(p)*albsni_hst(c,1) fsr_sno_ni(p) = fsds_nir_i(p)*albsni_hst(c,2) else fsds_sno_vd(p) = spval fsds_sno_nd(p) = spval fsds_sno_vi(p) = spval fsds_sno_ni(p) = spval fsr_sno_vd(p) = spval fsr_sno_nd(p) = spval fsr_sno_vi(p) = spval fsr_sno_ni(p) = spval endif end if end do end subroutine SurfaceRadiation end module SurfaceRadiationMod