module SurfaceAlbedoMod !----------------------------------------------------------------------- !BOP ! ! !MODULE: SurfaceAlbedoMod ! ! !DESCRIPTION: ! Performs surface albedo calculations ! ! !PUBLIC TYPES: use clm_varcon , only : istsoil use clm_varpar , only : numrad use clm_varcon , only : istcrop use shr_kind_mod, only : r8 => shr_kind_r8 use clm_varpar , only : nlevsno use SNICARMod , only : sno_nbr_aer, SNICAR_RT, DO_SNO_AER, DO_SNO_OC implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: SurfaceAlbedo ! Surface albedo and two-stream fluxes ! ! !PUBLIC DATA MEMBERS: ! The CLM default albice values are too high. ! Full-spectral albedo for land ice is ~0.5 (Paterson, Physics of Glaciers, 1994, p. 59) ! This is the value used in CAM3 by Pritchard et al., GRL, 35, 2008. real(r8), public :: albice(numrad) = & ! albedo land ice by waveband (1=vis, 2=nir) (/ 0.80_r8, 0.55_r8 /) ! ! !PRIVATE MEMBER FUNCTIONS: private :: SoilAlbedo ! Determine ground surface albedo private :: TwoStream ! Two-stream fluxes for canopy radiative transfer ! ! !REVISION HISTORY: ! Created by Mariana Vertenstein ! !EOP !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: SurfaceAlbedo ! ! !INTERFACE: subroutine SurfaceAlbedo(lbg, ubg, lbc, ubc, lbp, ubp, & num_nourbanc, filter_nourbanc, & num_nourbanp, filter_nourbanp, & nextsw_cday, declinp1) ! ! !DESCRIPTION: ! Surface albedo and two-stream fluxes ! Surface albedos. Also fluxes (per unit incoming direct and diffuse ! radiation) reflected, transmitted, and absorbed by vegetation. ! Also sunlit fraction of the canopy. ! The calling sequence is: ! -> SurfaceAlbedo: albedos for next time step ! -> SoilAlbedo: soil/lake/glacier/wetland albedos ! -> SNICAR_RT: snow albedos: direct beam (SNICAR) ! -> SNICAR_RT: snow albedos: diffuse (SNICAR) ! -> TwoStream: absorbed, reflected, transmitted solar fluxes (vis dir,vis dif, nir dir, nir dif) ! ! !USES: use clmtype use shr_orb_mod use clm_time_manager, only : get_nstep ! ! !ARGUMENTS: implicit none integer , intent(in) :: lbg, ubg ! gridcell bounds integer , intent(in) :: lbc, ubc ! column bounds integer , intent(in) :: lbp, ubp ! pft bounds integer , intent(in) :: num_nourbanc ! number of columns in non-urban filter integer , intent(in) :: filter_nourbanc(ubc-lbc+1) ! column filter for non-urban points integer , intent(in) :: num_nourbanp ! number of pfts in non-urban filter integer , intent(in) :: filter_nourbanp(ubp-lbp+1) ! pft filter for non-urban points real(r8), intent(in) :: nextsw_cday ! calendar day at Greenwich (1.00, ..., days/year) real(r8), intent(in) :: declinp1 ! declination angle (radians) for next time step ! ! !CALLED FROM: ! subroutine clm_driver1 ! subroutine iniTimeVar ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 2/1/02, Peter Thornton: Migrate to new data structures ! 8/20/03, Mariana Vertenstein: Vectorized routine ! 11/3/03, Peter Thornton: added decl(c) output for use in CN code. ! 03/28/08, Mark Flanner: added SNICAR, which required reversing the ! order of calls to SNICAR_RT and SoilAlbedo and the location where ! ground albedo is calculated ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in arguments ! integer , pointer :: pgridcell(:) ! gridcell of corresponding pft integer , pointer :: plandunit(:) ! index into landunit level quantities integer , pointer :: itypelun(:) ! landunit type integer , pointer :: pcolumn(:) ! column of corresponding pft integer , pointer :: cgridcell(:) ! gridcell of corresponding column real(r8), pointer :: pwtgcell(:) ! weight of pft wrt corresponding gridcell real(r8), pointer :: lat(:) ! gridcell latitude (radians) real(r8), pointer :: lon(:) ! gridcell longitude (radians) 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 :: h2osno(:) ! snow water (mm H2O) real(r8), pointer :: rhol(:,:) ! leaf reflectance: 1=vis, 2=nir real(r8), pointer :: rhos(:,:) ! stem reflectance: 1=vis, 2=nir real(r8), pointer :: taul(:,:) ! leaf transmittance: 1=vis, 2=nir real(r8), pointer :: taus(:,:) ! stem transmittance: 1=vis, 2=nir integer , pointer :: ivt(:) ! pft vegetation type ! ! local pointers toimplicit out arguments ! real(r8), pointer :: coszen(:) ! cosine of solar zenith angle real(r8), pointer :: fsun(:) ! sunlit fraction of canopy 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 :: 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 :: decl(:) ! solar declination angle (radians) 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 :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water content (col,lyr) [kg/m2] real(r8), pointer :: h2osoi_ice(:,:) ! ice lens content (col,lyr) [kg/m2] real(r8), pointer :: mss_cnc_bcphi(:,:) ! mass concentration of hydrophilic BC (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_bcpho(:,:) ! mass concentration of hydrophobic BC (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_ocphi(:,:) ! mass concentration of hydrophilic OC (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_ocpho(:,:) ! mass concentration of hydrophobic OC (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_dst1(:,:) ! mass concentration of dust aerosol species 1 (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_dst2(:,:) ! mass concentration of dust aerosol species 2 (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_dst3(:,:) ! mass concentration of dust aerosol species 3 (col,lyr) [kg/kg] real(r8), pointer :: mss_cnc_dst4(:,:) ! mass concentration of dust aerosol species 4 (col,lyr) [kg/kg] real(r8), pointer :: albsod(:,:) ! direct-beam soil albedo (col,bnd) [frc] real(r8), pointer :: albsoi(:,:) ! diffuse soil albedo (col,bnd) [frc] 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] real(r8), pointer :: snw_rds(:,:) ! snow grain radius (col,lyr) [microns] 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) real(r8), pointer :: albgri_bc(:,:) ! ground albedo without BC (diffuse) real(r8), pointer :: albgrd_oc(:,:) ! ground albedo without OC (direct) real(r8), pointer :: albgri_oc(:,:) ! ground albedo without OC (diffuse) real(r8), pointer :: albgrd_dst(:,:) ! ground albedo without dust (direct) real(r8), pointer :: albgri_dst(:,:) ! ground albedo without dust (diffuse) 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) [frc] ! ! ! !OTHER LOCAL VARIABLES: !EOP ! real(r8), parameter :: mpe = 1.e-06_r8 ! prevents overflow for division by zero integer :: fp,fc,g,c,p ! indices integer :: ib ! band index integer :: ic ! 0=unit incoming direct; 1=unit incoming diffuse real(r8) :: wl(lbp:ubp) ! fraction of LAI+SAI that is LAI real(r8) :: ws(lbp:ubp) ! fraction of LAI+SAI that is SAI real(r8) :: vai(lbp:ubp) ! elai+esai real(r8) :: rho(lbp:ubp,numrad) ! leaf/stem refl weighted by fraction LAI and SAI real(r8) :: tau(lbp:ubp,numrad) ! leaf/stem tran weighted by fraction LAI and SAI real(r8) :: ftdi(lbp:ubp,numrad) ! down direct flux below veg per unit dif flux = 0 real(r8) :: albsnd(lbc:ubc,numrad) ! snow albedo (direct) real(r8) :: albsni(lbc:ubc,numrad) ! snow albedo (diffuse) real(r8) :: ext(lbp:ubp) ! optical depth direct beam per unit LAI+SAI real(r8) :: coszen_gcell(lbg:ubg) ! cosine solar zenith angle for next time step (gridcell level) real(r8) :: coszen_col(lbc:ubc) ! cosine solar zenith angle for next time step (pft level) real(r8) :: coszen_pft(lbp:ubp) ! cosine solar zenith angle for next time step (pft level) integer :: num_vegsol ! number of vegetated pfts where coszen>0 integer :: filter_vegsol(ubp-lbp+1) ! pft filter where vegetated and coszen>0 integer :: num_novegsol ! number of vegetated pfts where coszen>0 integer :: filter_novegsol(ubp-lbp+1) ! pft filter where vegetated and coszen>0 integer, parameter :: nband =numrad ! number of solar radiation waveband classes integer :: flg_slr ! flag for SNICAR (=1 if direct, =2 if diffuse) integer :: flg_snw_ice ! flag for SNICAR (=1 when called from CLM, =2 when called from sea-ice) real(r8) :: albsnd_pur(lbc:ubc,numrad) ! direct pure snow albedo (radiative forcing) real(r8) :: albsni_pur(lbc:ubc,numrad) ! diffuse pure snow albedo (radiative forcing) real(r8) :: albsnd_bc(lbc:ubc,numrad) ! direct snow albedo without BC (radiative forcing) real(r8) :: albsni_bc(lbc:ubc,numrad) ! diffuse snow albedo without BC (radiative forcing) real(r8) :: albsnd_oc(lbc:ubc,numrad) ! direct snow albedo without OC (radiative forcing) real(r8) :: albsni_oc(lbc:ubc,numrad) ! diffuse snow albedo without OC (radiative forcing) real(r8) :: albsnd_dst(lbc:ubc,numrad) ! direct snow albedo without dust (radiative forcing) real(r8) :: albsni_dst(lbc:ubc,numrad) ! diffuse snow albedo without dust (radiative forcing) integer :: i ! index for layers [idx] real(r8) :: flx_absd_snw(lbc:ubc,-nlevsno+1:1,numrad) ! flux absorption factor for just snow (direct) [frc] real(r8) :: flx_absi_snw(lbc:ubc,-nlevsno+1:1,numrad) ! flux absorption factor for just snow (diffuse) [frc] real(r8) :: foo_snw(lbc:ubc,-nlevsno+1:1,numrad) ! dummy array for forcing calls real(r8) :: albsfc(lbc:ubc,numrad) ! albedo of surface underneath snow (col,bnd) real(r8) :: h2osno_liq(lbc:ubc,-nlevsno+1:0) ! liquid snow content (col,lyr) [kg m-2] real(r8) :: h2osno_ice(lbc:ubc,-nlevsno+1:0) ! ice content in snow (col,lyr) [kg m-2] integer :: snw_rds_in(lbc:ubc,-nlevsno+1:0) ! snow grain size sent to SNICAR (col,lyr) [microns] real(r8) :: mss_cnc_aer_in_frc_pur(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for forcing calculation (zero) (col,lyr,aer) [kg kg-1] real(r8) :: mss_cnc_aer_in_frc_bc(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for BC forcing (col,lyr,aer) [kg kg-1] real(r8) :: mss_cnc_aer_in_frc_oc(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for OC forcing (col,lyr,aer) [kg kg-1] real(r8) :: mss_cnc_aer_in_frc_dst(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for dust forcing (col,lyr,aer) [kg kg-1] real(r8) :: mss_cnc_aer_in_fdb(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of all aerosol species for feedback calculation (col,lyr,aer) [kg kg-1] !----------------------------------------------------------------------- ! Assign local pointers to derived subtypes components (gridcell-level) lat => clm3%g%lat lon => clm3%g%lon ! Assign local pointers to derived subtypes components (landunit level) itypelun => clm3%g%l%itype ! Assign local pointers to derived subtypes components (column-level) cgridcell => clm3%g%l%c%gridcell h2osno => clm3%g%l%c%cws%h2osno albgrd => clm3%g%l%c%cps%albgrd albgri => clm3%g%l%c%cps%albgri decl => clm3%g%l%c%cps%decl coszen => clm3%g%l%c%cps%coszen albsod => clm3%g%l%c%cps%albsod albsoi => clm3%g%l%c%cps%albsoi 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 h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice snw_rds => clm3%g%l%c%cps%snw_rds albgrd_pur => clm3%g%l%c%cps%albgrd_pur albgri_pur => clm3%g%l%c%cps%albgri_pur albgrd_bc => clm3%g%l%c%cps%albgrd_bc albgri_bc => clm3%g%l%c%cps%albgri_bc albgrd_oc => clm3%g%l%c%cps%albgrd_oc albgri_oc => clm3%g%l%c%cps%albgri_oc albgrd_dst => clm3%g%l%c%cps%albgrd_dst albgri_dst => clm3%g%l%c%cps%albgri_dst mss_cnc_bcphi => clm3%g%l%c%cps%mss_cnc_bcphi mss_cnc_bcpho => clm3%g%l%c%cps%mss_cnc_bcpho mss_cnc_ocphi => clm3%g%l%c%cps%mss_cnc_ocphi mss_cnc_ocpho => clm3%g%l%c%cps%mss_cnc_ocpho mss_cnc_dst1 => clm3%g%l%c%cps%mss_cnc_dst1 mss_cnc_dst2 => clm3%g%l%c%cps%mss_cnc_dst2 mss_cnc_dst3 => clm3%g%l%c%cps%mss_cnc_dst3 mss_cnc_dst4 => clm3%g%l%c%cps%mss_cnc_dst4 albsnd_hst => clm3%g%l%c%cps%albsnd_hst albsni_hst => clm3%g%l%c%cps%albsni_hst ! Assign local pointers to derived subtypes components (pft-level) plandunit => clm3%g%l%c%p%landunit pgridcell => clm3%g%l%c%p%gridcell pcolumn => clm3%g%l%c%p%column pwtgcell => clm3%g%l%c%p%wtgcell albd => clm3%g%l%c%p%pps%albd albi => clm3%g%l%c%p%pps%albi 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 fsun => clm3%g%l%c%p%pps%fsun elai => clm3%g%l%c%p%pps%elai esai => clm3%g%l%c%p%pps%esai gdir => clm3%g%l%c%p%pps%gdir omega => clm3%g%l%c%p%pps%omega ivt => clm3%g%l%c%p%itype rhol => pftcon%rhol rhos => pftcon%rhos taul => pftcon%taul taus => pftcon%taus ! Cosine solar zenith angle for next time step do g = lbg, ubg coszen_gcell(g) = shr_orb_cosz (nextsw_cday, lat(g), lon(g), declinp1) end do ! Save coszen and declination values to clm3 data structures for ! use in other places in the CN and urban code do c = lbc,ubc g = cgridcell(c) coszen_col(c) = coszen_gcell(g) coszen(c) = coszen_col(c) decl(c) = declinp1 end do do fp = 1,num_nourbanp p = filter_nourbanp(fp) ! if (pwtgcell(p)>0._r8) then ! "if" added due to chg in filter definition g = pgridcell(p) coszen_pft(p) = coszen_gcell(g) ! end if ! then removed for CNDV (and dyn. landuse?) cases to work end do ! Initialize output because solar radiation only done if coszen > 0 do ib = 1, numrad do fc = 1,num_nourbanc c = filter_nourbanc(fc) albgrd(c,ib) = 0._r8 albgri(c,ib) = 0._r8 albgrd_pur(c,ib) = 0._r8 albgri_pur(c,ib) = 0._r8 albgrd_bc(c,ib) = 0._r8 albgri_bc(c,ib) = 0._r8 albgrd_oc(c,ib) = 0._r8 albgri_oc(c,ib) = 0._r8 albgrd_dst(c,ib) = 0._r8 albgri_dst(c,ib) = 0._r8 do i=-nlevsno+1,1,1 flx_absdv(c,i) = 0._r8 flx_absdn(c,i) = 0._r8 flx_absiv(c,i) = 0._r8 flx_absin(c,i) = 0._r8 enddo end do do fp = 1,num_nourbanp p = filter_nourbanp(fp) ! if (pwtgcell(p)>0._r8) then ! "if" added due to chg in filter definition albd(p,ib) = 1._r8 albi(p,ib) = 1._r8 fabd(p,ib) = 0._r8 fabi(p,ib) = 0._r8 ftdd(p,ib) = 0._r8 ftid(p,ib) = 0._r8 ftii(p,ib) = 0._r8 omega(p,ib)= 0._r8 if (ib==1) then gdir(p) = 0._r8 end if ! end if ! then removed for CNDV (and dyn. landuse?) cases to work end do end do ! SoilAlbedo called before SNICAR_RT ! so that reflectance of soil beneath snow column is known ! ahead of time for snow RT calculation. ! Snow albedos ! Note that snow albedo routine will only compute nonzero snow albedos ! where h2osno> 0 and coszen > 0 ! Ground surface albedos ! Note that ground albedo routine will only compute nonzero snow albedos ! where coszen > 0 call SoilAlbedo(lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, albsnd, albsni) ! set variables to pass to SNICAR. flg_snw_ice = 1 ! calling from CLM, not CSIM do c=lbc,ubc albsfc(c,:) = albsoi(c,:) h2osno_liq(c,:) = h2osoi_liq(c,-nlevsno+1:0) h2osno_ice(c,:) = h2osoi_ice(c,-nlevsno+1:0) snw_rds_in(c,:) = nint(snw_rds(c,:)) ! zero aerosol input arrays mss_cnc_aer_in_frc_pur(c,:,:) = 0._r8 mss_cnc_aer_in_frc_bc(c,:,:) = 0._r8 mss_cnc_aer_in_frc_oc(c,:,:) = 0._r8 mss_cnc_aer_in_frc_dst(c,:,:) = 0._r8 mss_cnc_aer_in_fdb(c,:,:) = 0._r8 end do ! Set aerosol input arrays ! feedback input arrays have been zeroed ! set soot and dust aerosol concentrations: if (DO_SNO_AER) then mss_cnc_aer_in_fdb(lbc:ubc,:,1) = mss_cnc_bcphi(lbc:ubc,:) mss_cnc_aer_in_fdb(lbc:ubc,:,2) = mss_cnc_bcpho(lbc:ubc,:) ! DO_SNO_OC is set in SNICAR_varpar. Default case is to ignore OC concentrations because: ! 1) Knowledge of their optical properties is primitive ! 2) When 'water-soluble' OPAC optical properties are applied to OC in snow, ! it has a negligible darkening effect. if (DO_SNO_OC) then mss_cnc_aer_in_fdb(lbc:ubc,:,3) = mss_cnc_ocphi(lbc:ubc,:) mss_cnc_aer_in_fdb(lbc:ubc,:,4) = mss_cnc_ocpho(lbc:ubc,:) endif mss_cnc_aer_in_fdb(lbc:ubc,:,5) = mss_cnc_dst1(lbc:ubc,:) mss_cnc_aer_in_fdb(lbc:ubc,:,6) = mss_cnc_dst2(lbc:ubc,:) mss_cnc_aer_in_fdb(lbc:ubc,:,7) = mss_cnc_dst3(lbc:ubc,:) mss_cnc_aer_in_fdb(lbc:ubc,:,8) = mss_cnc_dst4(lbc:ubc,:) endif ! If radiative forcing is being calculated, first estimate clean-snow albedo #if (defined SNICAR_FRC) ! 1. BC input array: ! set dust and (optionally) OC concentrations, so BC_FRC=[(BC+OC+dust)-(OC+dust)] mss_cnc_aer_in_frc_bc(lbc:ubc,:,5) = mss_cnc_dst1(lbc:ubc,:) mss_cnc_aer_in_frc_bc(lbc:ubc,:,6) = mss_cnc_dst2(lbc:ubc,:) mss_cnc_aer_in_frc_bc(lbc:ubc,:,7) = mss_cnc_dst3(lbc:ubc,:) mss_cnc_aer_in_frc_bc(lbc:ubc,:,8) = mss_cnc_dst4(lbc:ubc,:) if (DO_SNO_OC) then mss_cnc_aer_in_frc_bc(lbc:ubc,:,3) = mss_cnc_ocphi(lbc:ubc,:) mss_cnc_aer_in_frc_bc(lbc:ubc,:,4) = mss_cnc_ocpho(lbc:ubc,:) endif ! BC FORCING CALCULATIONS flg_slr = 1; ! direct-beam call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_bc, albsfc, albsnd_bc, foo_snw) flg_slr = 2; ! diffuse call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_bc, albsfc, albsni_bc, foo_snw) ! 2. OC input array: ! set BC and dust concentrations, so OC_FRC=[(BC+OC+dust)-(BC+dust)] if (DO_SNO_OC) then mss_cnc_aer_in_frc_oc(lbc:ubc,:,1) = mss_cnc_bcphi(lbc:ubc,:) mss_cnc_aer_in_frc_oc(lbc:ubc,:,2) = mss_cnc_bcpho(lbc:ubc,:) mss_cnc_aer_in_frc_oc(lbc:ubc,:,5) = mss_cnc_dst1(lbc:ubc,:) mss_cnc_aer_in_frc_oc(lbc:ubc,:,6) = mss_cnc_dst2(lbc:ubc,:) mss_cnc_aer_in_frc_oc(lbc:ubc,:,7) = mss_cnc_dst3(lbc:ubc,:) mss_cnc_aer_in_frc_oc(lbc:ubc,:,8) = mss_cnc_dst4(lbc:ubc,:) ! OC FORCING CALCULATIONS flg_slr = 1; ! direct-beam call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_oc, albsfc, albsnd_oc, foo_snw) flg_slr = 2; ! diffuse call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_oc, albsfc, albsni_oc, foo_snw) endif ! 3. DUST input array: ! set BC and OC concentrations, so DST_FRC=[(BC+OC+dust)-(BC+OC)] mss_cnc_aer_in_frc_dst(lbc:ubc,:,1) = mss_cnc_bcphi(lbc:ubc,:) mss_cnc_aer_in_frc_dst(lbc:ubc,:,2) = mss_cnc_bcpho(lbc:ubc,:) if (DO_SNO_OC) then mss_cnc_aer_in_frc_dst(lbc:ubc,:,3) = mss_cnc_ocphi(lbc:ubc,:) mss_cnc_aer_in_frc_dst(lbc:ubc,:,4) = mss_cnc_ocpho(lbc:ubc,:) endif ! DUST FORCING CALCULATIONS flg_slr = 1; ! direct-beam call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_dst, albsfc, albsnd_dst, foo_snw) flg_slr = 2; ! diffuse call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_dst, albsfc, albsni_dst, foo_snw) ! 4. ALL AEROSOL FORCING CALCULATION ! (pure snow albedo) flg_slr = 1; ! direct-beam call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_pur, albsfc, albsnd_pur, foo_snw) flg_slr = 2; ! diffuse call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_frc_pur, albsfc, albsni_pur, foo_snw) #endif ! CLIMATE FEEDBACK CALCULATIONS, ALL AEROSOLS: flg_slr = 1; ! direct-beam call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_fdb, albsfc, albsnd, flx_absd_snw) flg_slr = 2; ! diffuse call SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & coszen_col, flg_slr, h2osno_liq, h2osno_ice, snw_rds_in, & mss_cnc_aer_in_fdb, albsfc, albsni, flx_absi_snw) ! ground albedos and snow-fraction weighting of snow absorption factors do ib = 1, nband do fc = 1,num_nourbanc c = filter_nourbanc(fc) if (coszen(c) > 0._r8) then ! ground albedo was originally computed in SoilAlbedo, but is now computed here ! because the order of SoilAlbedo and SNICAR_RT was switched for SNICAR. albgrd(c,ib) = albsod(c,ib)*(1._r8-frac_sno(c)) + albsnd(c,ib)*frac_sno(c) albgri(c,ib) = albsoi(c,ib)*(1._r8-frac_sno(c)) + albsni(c,ib)*frac_sno(c) ! albedos for radiative forcing calculations: #if (defined SNICAR_FRC) ! BC forcing albedo albgrd_bc(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_bc(c,ib)*frac_sno(c) albgri_bc(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_bc(c,ib)*frac_sno(c) if (DO_SNO_OC) then ! OC forcing albedo albgrd_oc(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_oc(c,ib)*frac_sno(c) albgri_oc(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_oc(c,ib)*frac_sno(c) endif ! dust forcing albedo albgrd_dst(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_dst(c,ib)*frac_sno(c) albgri_dst(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_dst(c,ib)*frac_sno(c) ! pure snow albedo for all-aerosol radiative forcing albgrd_pur(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_pur(c,ib)*frac_sno(c) albgri_pur(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_pur(c,ib)*frac_sno(c) #endif ! also in this loop (but optionally in a different loop for vectorized code) ! weight snow layer radiative absorption factors based on snow fraction and soil albedo ! (NEEDED FOR ENERGY CONSERVATION) do i = -nlevsno+1,1,1 if (ib == 1) then flx_absdv(c,i) = flx_absd_snw(c,i,ib)*frac_sno(c) + & ((1.-frac_sno(c))*(1-albsod(c,ib))*(flx_absd_snw(c,i,ib)/(1.-albsnd(c,ib)))) flx_absiv(c,i) = flx_absi_snw(c,i,ib)*frac_sno(c) + & ((1.-frac_sno(c))*(1-albsoi(c,ib))*(flx_absi_snw(c,i,ib)/(1.-albsni(c,ib)))) elseif (ib == 2) then flx_absdn(c,i) = flx_absd_snw(c,i,ib)*frac_sno(c) + & ((1.-frac_sno(c))*(1-albsod(c,ib))*(flx_absd_snw(c,i,ib)/(1.-albsnd(c,ib)))) flx_absin(c,i) = flx_absi_snw(c,i,ib)*frac_sno(c) + & ((1.-frac_sno(c))*(1-albsoi(c,ib))*(flx_absi_snw(c,i,ib)/(1.-albsni(c,ib)))) endif enddo endif enddo enddo ! for diagnostics, set snow albedo to spval over non-snow points ! so that it is not averaged in history buffer ! (OPTIONAL) do ib = 1, nband do fc = 1,num_nourbanc c = filter_nourbanc(fc) if ((coszen(c) > 0._r8) .and. (h2osno(c) > 0._r8)) then albsnd_hst(c,ib) = albsnd(c,ib) albsni_hst(c,ib) = albsni(c,ib) else albsnd_hst(c,ib) = 0._r8 albsni_hst(c,ib) = 0._r8 endif enddo enddo ! Create solar-vegetated filter for the following calculations num_vegsol = 0 num_novegsol = 0 do fp = 1,num_nourbanp p = filter_nourbanp(fp) if (coszen_pft(p) > 0._r8) then if ((itypelun(plandunit(p)) == istsoil .or. & itypelun(plandunit(p)) == istcrop ) & .and. (elai(p) + esai(p)) > 0._r8 & .and. pwtgcell(p) > 0._r8) then num_vegsol = num_vegsol + 1 filter_vegsol(num_vegsol) = p else num_novegsol = num_novegsol + 1 filter_novegsol(num_novegsol) = p end if end if end do ! Weight reflectance/transmittance by lai and sai ! Only perform on vegetated pfts where coszen > 0 do fp = 1,num_vegsol p = filter_vegsol(fp) vai(p) = elai(p) + esai(p) wl(p) = elai(p) / max( vai(p), mpe ) ws(p) = esai(p) / max( vai(p), mpe ) end do do ib = 1, numrad do fp = 1,num_vegsol p = filter_vegsol(fp) rho(p,ib) = max( rhol(ivt(p),ib)*wl(p) + rhos(ivt(p),ib)*ws(p), mpe ) tau(p,ib) = max( taul(ivt(p),ib)*wl(p) + taus(ivt(p),ib)*ws(p), mpe ) end do end do ! Calculate surface albedos and fluxes ! Only perform on vegetated pfts where coszen > 0 call TwoStream (lbc, ubc, lbp, ubp, filter_vegsol, num_vegsol, & coszen_pft, vai, rho, tau) ! Determine values for non-vegetated pfts where coszen > 0 do ib = 1,numrad do fp = 1,num_novegsol p = filter_novegsol(fp) c = pcolumn(p) fabd(p,ib) = 0._r8 fabi(p,ib) = 0._r8 ftdd(p,ib) = 1._r8 ftid(p,ib) = 0._r8 ftii(p,ib) = 1._r8 albd(p,ib) = albgrd(c,ib) albi(p,ib) = albgri(c,ib) gdir(p) = 0._r8 end do end do end subroutine SurfaceAlbedo !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: SoilAlbedo ! ! !INTERFACE: subroutine SoilAlbedo (lbc, ubc, num_nourbanc, filter_nourbanc, coszen, albsnd, albsni) ! ! !DESCRIPTION: ! Determine ground surface albedo, accounting for snow ! ! !USES: use clmtype use clm_varpar, only : numrad use clm_varcon, only : albsat, albdry, alblak, tfrz, istice, istice_mec ! ! !ARGUMENTS: implicit none integer , intent(in) :: lbc, ubc ! column bounds integer , intent(in) :: num_nourbanc ! number of columns in non-urban points in column filter integer , intent(in) :: filter_nourbanc(ubc-lbc+1) ! column filter for non-urban points real(r8), intent(in) :: coszen(lbc:ubc) ! cos solar zenith angle next time step (column-level) real(r8), intent(in) :: albsnd(lbc:ubc,numrad) ! snow albedo (direct) real(r8), intent(in) :: albsni(lbc:ubc,numrad) ! snow albedo (diffuse) ! ! !CALLED FROM: ! subroutine SurfaceAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 2/5/02, Peter Thornton: Migrated to new data structures. ! 8/20/03, Mariana Vertenstein: Vectorized routine ! 03/28/08, Mark Flanner: changes for SNICAR ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments ! integer , pointer :: clandunit(:) ! landunit of corresponding column integer , pointer :: ltype(:) ! landunit type integer , pointer :: isoicol(:) ! soil color class real(r8), pointer :: t_grnd(:) ! ground temperature (Kelvin) real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water [m3/m3] ! ! local pointers to original implicit out arguments ! real(r8), pointer:: albgrd(:,:) ! ground albedo (direct) real(r8), pointer:: albgri(:,:) ! ground albedo (diffuse) ! albsod and albsoi are now clm_type variables so they can be used by SNICAR. real(r8), pointer :: albsod(:,:) ! soil albedo (direct) real(r8), pointer :: albsoi(:,:) ! soil albedo (diffuse) ! ! ! !OTHER LOCAL VARIABLES: !EOP ! integer, parameter :: nband =numrad ! number of solar radiation waveband classes integer :: fc ! non-urban filter column index integer :: c,l ! indices integer :: ib ! waveband number (1=vis, 2=nir) real(r8) :: inc ! soil water correction factor for soil albedo ! albsod and albsoi are now clm_type variables so they can be used by SNICAR. !real(r8) :: albsod ! soil albedo (direct) !real(r8) :: albsoi ! soil albedo (diffuse) integer :: soilcol ! soilcolor !----------------------------------------------------------------------- !dir$ inlinenever SoilAlbedo ! Assign local pointers to derived subtypes components (column-level) clandunit => clm3%g%l%c%landunit isoicol => clm3%g%l%c%cps%isoicol t_grnd => clm3%g%l%c%ces%t_grnd frac_sno => clm3%g%l%c%cps%frac_sno h2osoi_vol => clm3%g%l%c%cws%h2osoi_vol albgrd => clm3%g%l%c%cps%albgrd albgri => clm3%g%l%c%cps%albgri albsod => clm3%g%l%c%cps%albsod albsoi => clm3%g%l%c%cps%albsoi ! Assign local pointers to derived subtypes components (landunit-level) ltype => clm3%g%l%itype ! Compute soil albedos do ib = 1, nband do fc = 1,num_nourbanc c = filter_nourbanc(fc) if (coszen(c) > 0._r8) then l = clandunit(c) if (ltype(l) == istsoil .or. ltype(l) == istcrop) then ! soil inc = max(0.11_r8-0.40_r8*h2osoi_vol(c,1), 0._r8) soilcol = isoicol(c) ! changed from local variable to clm_type: !albsod = min(albsat(soilcol,ib)+inc, albdry(soilcol,ib)) !albsoi = albsod albsod(c,ib) = min(albsat(soilcol,ib)+inc, albdry(soilcol,ib)) albsoi(c,ib) = albsod(c,ib) else if (ltype(l) == istice .or. ltype(l) == istice_mec) then ! land ice ! changed from local variable to clm_type: !albsod = albice(ib) !albsoi = albsod albsod(c,ib) = albice(ib) albsoi(c,ib) = albsod(c,ib) else if (t_grnd(c) > tfrz) then ! unfrozen lake, wetland ! changed from local variable to clm_type: !albsod = 0.05_r8/(max(0.001_r8,coszen(c)) + 0.15_r8) !albsoi = albsod albsod(c,ib) = 0.05_r8/(max(0.001_r8,coszen(c)) + 0.15_r8) albsoi(c,ib) = albsod(c,ib) else ! frozen lake, wetland ! changed from local variable to clm_type: !albsod = alblak(ib) !albsoi = albsod albsod(c,ib) = alblak(ib) albsoi(c,ib) = albsod(c,ib) end if ! Weighting is done in SurfaceAlbedo, after the call to SNICAR_RT ! This had to be done, because SoilAlbedo is called before SNICAR_RT, so at ! this point, snow albedo is not yet known. end if end do end do end subroutine SoilAlbedo !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: TwoStream ! ! !INTERFACE: subroutine TwoStream (lbc, ubc, lbp, ubp, filter_vegsol, num_vegsol, & coszen, vai, rho, tau) ! ! !DESCRIPTION: ! Two-stream fluxes for canopy radiative transfer ! Use two-stream approximation of Dickinson (1983) Adv Geophysics ! 25:305-353 and Sellers (1985) Int J Remote Sensing 6:1335-1372 ! to calculate fluxes absorbed by vegetation, reflected by vegetation, ! and transmitted through vegetation for unit incoming direct or diffuse ! flux given an underlying surface with known albedo. ! ! !USES: use clmtype use clm_varpar, only : numrad use clm_varcon, only : omegas, tfrz, betads, betais ! ! !ARGUMENTS: implicit none integer , intent(in) :: lbc, ubc ! column bounds integer , intent(in) :: lbp, ubp ! pft bounds integer , intent(in) :: filter_vegsol(ubp-lbp+1) ! filter for vegetated pfts with coszen>0 integer , intent(in) :: num_vegsol ! number of vegetated pfts where coszen>0 real(r8), intent(in) :: coszen(lbp:ubp) ! cosine solar zenith angle for next time step real(r8), intent(in) :: vai(lbp:ubp) ! elai+esai real(r8), intent(in) :: rho(lbp:ubp,numrad) ! leaf/stem refl weighted by fraction LAI and SAI real(r8), intent(in) :: tau(lbp:ubp,numrad) ! leaf/stem tran weighted by fraction LAI and SAI ! ! !CALLED FROM: ! subroutine SurfaceAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! Modified for speedup: Mariana Vertenstein, 8/26/02 ! Vectorized routine: Mariana Vertenstein: 8/20/03 ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in scalars ! integer , pointer :: pcolumn(:) ! column of corresponding pft real(r8), pointer :: albgrd(:,:) ! ground albedo (direct) (column-level) real(r8), pointer :: albgri(:,:) ! ground albedo (diffuse)(column-level) real(r8), pointer :: t_veg(:) ! vegetation temperature (Kelvin) real(r8), pointer :: fwet(:) ! fraction of canopy that is wet (0 to 1) integer , pointer :: ivt(:) ! pft vegetation type real(r8), pointer :: xl(:) ! ecophys const - leaf/stem orientation index ! ! local pointers to implicit out scalars ! real(r8), pointer :: albd(:,:) ! surface albedo (direct) real(r8), pointer :: albi(:,:) ! surface albedo (diffuse) 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 :: gdir(:) ! leaf projection in solar direction (0 to 1) real(r8), pointer :: omega(:,:) ! fraction of intercepted radiation that is scattered (0 to 1) ! ! ! !OTHER LOCAL VARIABLES: !EOP ! integer :: fp,p,c ! array indices !integer :: ic ! 0=unit incoming direct; 1=unit incoming diffuse integer :: ib ! waveband number real(r8) :: cosz ! 0.001 <= coszen <= 1.000 real(r8) :: asu ! single scattering albedo real(r8) :: chil(lbp:ubp) ! -0.4 <= xl <= 0.6 real(r8) :: twostext(lbp:ubp)! optical depth of direct beam per unit leaf area real(r8) :: avmu(lbp:ubp) ! average diffuse optical depth real(r8) :: omegal ! omega for leaves real(r8) :: betai ! upscatter parameter for diffuse radiation real(r8) :: betail ! betai for leaves real(r8) :: betad ! upscatter parameter for direct beam radiation real(r8) :: betadl ! betad for leaves real(r8) :: tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9 ! temporary real(r8) :: p1,p2,p3,p4,s1,s2,u1,u2,u3 ! temporary real(r8) :: b,c1,d,d1,d2,f,h,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10 ! temporary real(r8) :: phi1,phi2,sigma ! temporary real(r8) :: temp0(lbp:ubp),temp1,temp2(lbp:ubp) ! temporary real(r8) :: t1 !----------------------------------------------------------------------- ! Assign local pointers to derived subtypes components (column-level) albgrd => clm3%g%l%c%cps%albgrd albgri => clm3%g%l%c%cps%albgri ! Assign local pointers to derived subtypes components (pft-level) pcolumn => clm3%g%l%c%p%column fwet => clm3%g%l%c%p%pps%fwet t_veg => clm3%g%l%c%p%pes%t_veg ivt => clm3%g%l%c%p%itype albd => clm3%g%l%c%p%pps%albd albi => clm3%g%l%c%p%pps%albi 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 gdir => clm3%g%l%c%p%pps%gdir omega => clm3%g%l%c%p%pps%omega xl => pftcon%xl ! Calculate two-stream parameters omega, betad, betai, avmu, gdir, twostext. ! Omega, betad, betai are adjusted for snow. Values for omega*betad ! and omega*betai are calculated and then divided by the new omega ! because the product omega*betai, omega*betad is used in solution. ! Also, the transmittances and reflectances (tau, rho) are linear ! weights of leaf and stem values. do fp = 1,num_vegsol p = filter_vegsol(fp) ! note that the following limit only acts on cosz values > 0 and less than ! 0.001, not on values cosz = 0, since these zero have already been filtered ! out in filter_vegsol cosz = max(0.001_r8, coszen(p)) chil(p) = min( max(xl(ivt(p)), -0.4_r8), 0.6_r8 ) if (abs(chil(p)) <= 0.01_r8) chil(p) = 0.01_r8 phi1 = 0.5_r8 - 0.633_r8*chil(p) - 0.330_r8*chil(p)*chil(p) phi2 = 0.877_r8 * (1._r8-2._r8*phi1) gdir(p) = phi1 + phi2*cosz twostext(p) = gdir(p)/cosz avmu(p) = ( 1._r8 - phi1/phi2 * log((phi1+phi2)/phi1) ) / phi2 temp0(p) = gdir(p) + phi2*cosz temp1 = phi1*cosz temp2(p) = ( 1._r8 - temp1/temp0(p) * log((temp1+temp0(p))/temp1) ) end do do ib = 1, numrad do fp = 1,num_vegsol p = filter_vegsol(fp) c = pcolumn(p) omegal = rho(p,ib) + tau(p,ib) asu = 0.5_r8*omegal*gdir(p)/temp0(p) *temp2(p) betadl = (1._r8+avmu(p)*twostext(p))/(omegal*avmu(p)*twostext(p))*asu betail = 0.5_r8 * ((rho(p,ib)+tau(p,ib)) + (rho(p,ib)-tau(p,ib)) & * ((1._r8+chil(p))/2._r8)**2) / omegal ! Adjust omega, betad, and betai for intercepted snow if (t_veg(p) > tfrz) then !no snow tmp0 = omegal tmp1 = betadl tmp2 = betail else tmp0 = (1._r8-fwet(p))*omegal + fwet(p)*omegas(ib) tmp1 = ( (1._r8-fwet(p))*omegal*betadl + fwet(p)*omegas(ib)*betads ) / tmp0 tmp2 = ( (1._r8-fwet(p))*omegal*betail + fwet(p)*omegas(ib)*betais ) / tmp0 end if omega(p,ib) = tmp0 betad = tmp1 betai = tmp2 ! Absorbed, reflected, transmitted fluxes per unit incoming radiation b = 1._r8 - omega(p,ib) + omega(p,ib)*betai c1 = omega(p,ib)*betai tmp0 = avmu(p)*twostext(p) d = tmp0 * omega(p,ib)*betad f = tmp0 * omega(p,ib)*(1._r8-betad) tmp1 = b*b - c1*c1 h = sqrt(tmp1) / avmu(p) sigma = tmp0*tmp0 - tmp1 p1 = b + avmu(p)*h p2 = b - avmu(p)*h p3 = b + tmp0 p4 = b - tmp0 ! PET, 03/01/04: added this test to avoid floating point errors in exp() ! EBK, 04/15/08: always do this for all modes -- not just CN t1 = min(h*vai(p), 40._r8) s1 = exp(-t1) t1 = min(twostext(p)*vai(p), 40._r8) s2 = exp(-t1) ! Determine fluxes for vegetated pft for unit incoming direct ! Loop over incoming direct and incoming diffuse ! 0=unit incoming direct; 1=unit incoming diffuse ! ic = 0 unit incoming direct flux ! ======================================== u1 = b - c1/albgrd(c,ib) u2 = b - c1*albgrd(c,ib) u3 = f + c1*albgrd(c,ib) tmp2 = u1 - avmu(p)*h tmp3 = u1 + avmu(p)*h d1 = p1*tmp2/s1 - p2*tmp3*s1 tmp4 = u2 + avmu(p)*h tmp5 = u2 - avmu(p)*h d2 = tmp4/s1 - tmp5*s1 h1 = -d*p4 - c1*f tmp6 = d - h1*p3/sigma tmp7 = ( d - c1 - h1/sigma*(u1+tmp0) ) * s2 h2 = ( tmp6*tmp2/s1 - p2*tmp7 ) / d1 h3 = - ( tmp6*tmp3*s1 - p1*tmp7 ) / d1 h4 = -f*p3 - c1*d tmp8 = h4/sigma tmp9 = ( u3 - tmp8*(u2-tmp0) ) * s2 h5 = - ( tmp8*tmp4/s1 + tmp9 ) / d2 h6 = ( tmp8*tmp5*s1 + tmp9 ) / d2 h7 = (c1*tmp2) / (d1*s1) h8 = (-c1*tmp3*s1) / d1 h9 = tmp4 / (d2*s1) h10 = (-tmp5*s1) / d2 ! Downward direct and diffuse fluxes below vegetation (ic = 0) ftdd(p,ib) = s2 ftid(p,ib) = h4*s2/sigma + h5*s1 + h6/s1 ! Flux reflected by vegetation (ic = 0) albd(p,ib) = h1/sigma + h2 + h3 ! Flux absorbed by vegetation (ic = 0) fabd(p,ib) = 1._r8 - albd(p,ib) & - (1._r8-albgrd(c,ib))*ftdd(p,ib) - (1._r8-albgri(c,ib))*ftid(p,ib) ! ic = 1 unit incoming diffuse ! ======================================== u1 = b - c1/albgri(c,ib) u2 = b - c1*albgri(c,ib) u3 = f + c1*albgri(c,ib) tmp2 = u1 - avmu(p)*h tmp3 = u1 + avmu(p)*h d1 = p1*tmp2/s1 - p2*tmp3*s1 tmp4 = u2 + avmu(p)*h tmp5 = u2 - avmu(p)*h d2 = tmp4/s1 - tmp5*s1 h1 = -d*p4 - c1*f tmp6 = d - h1*p3/sigma tmp7 = ( d - c1 - h1/sigma*(u1+tmp0) ) * s2 h2 = ( tmp6*tmp2/s1 - p2*tmp7 ) / d1 h3 = - ( tmp6*tmp3*s1 - p1*tmp7 ) / d1 h4 = -f*p3 - c1*d tmp8 = h4/sigma tmp9 = ( u3 - tmp8*(u2-tmp0) ) * s2 h5 = - ( tmp8*tmp4/s1 + tmp9 ) / d2 h6 = ( tmp8*tmp5*s1 + tmp9 ) / d2 h7 = (c1*tmp2) / (d1*s1) h8 = (-c1*tmp3*s1) / d1 h9 = tmp4 / (d2*s1) h10 = (-tmp5*s1) / d2 ! Downward direct and diffuse fluxes below vegetation ftii(p,ib) = h9*s1 + h10/s1 ! Flux reflected by vegetation albi(p,ib) = h7 + h8 ! Flux absorbed by vegetation fabi(p,ib) = 1._r8 - albi(p,ib) - (1._r8-albgri(c,ib))*ftii(p,ib) end do ! end of pft loop end do ! end of radiation band loop end subroutine TwoStream end module SurfaceAlbedoMod