diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 67206f73e1..8acb6e6c48 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1316,6 +1316,8 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l new_patch%frac_burnt = 0._r8 new_patch%total_tree_area = 0.0_r8 new_patch%NCL_p = 1 + + end subroutine create_patch @@ -1440,6 +1442,12 @@ subroutine zero_patch(cp_p) currentPatch%c_stomata = 0.0_r8 ! This is calculated immediately before use currentPatch%c_lblayer = 0.0_r8 + currentPatch%solar_zenith_flag = .false. + currentPatch%solar_zenith_angle = nan + + currentPatch%gnd_alb_dir(:) = nan + currentPatch%gnd_alb_dif(:) = nan + end subroutine zero_patch ! ============================================================================ diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 267d63d697..d9e4c6b3f2 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -14,6 +14,8 @@ module EDSurfaceRadiationMod use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : maxpft use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : pi_const use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : bc_out_type use FatesInterfaceMod , only : hlm_numSWb @@ -29,7 +31,8 @@ module EDSurfaceRadiationMod use EDTypesMod , only : ipar use EDCanopyStructureMod, only: calc_areaindex use FatesGlobals , only : fates_log - + use FatesGlobals, only : endrun => fates_endrun + ! CIME globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -37,6 +40,7 @@ module EDSurfaceRadiationMod private public :: ED_Norman_Radiation ! Surface albedo and two-stream fluxes + public :: PatchNormanRadiation public :: ED_SunShadeFracs logical :: debug = .false. ! for debugging this module @@ -66,84 +70,18 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) ! !LOCAL VARIABLES: - ! ============================================================================ - ! ED/NORMAN RADIATION DECS - ! ============================================================================ - type (ed_patch_type) , pointer :: currentPatch - integer :: radtype, L, ft, j, ifp - integer :: iter ! Iteration index - integer :: irep ! Flag to exit iteration loop - real(r8) :: sb - real(r8) :: error ! Error check - real(r8) :: down_rad, up_rad ! Iterative solution do Dif_dn and Dif_up - real(r8) :: ftweight(nclmax,maxpft,nlevleaf) - real(r8) :: k_dir(maxpft) ! Direct beam extinction coefficient - real(r8) :: tr_dir_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of direct beam radiation through a single layer - real(r8) :: tr_dif_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of diffuse radiation through a single layer - real(r8) :: forc_dir(maxPatchesPerSite,maxSWb) - real(r8) :: forc_dif(maxPatchesPerSite,maxSWb) - real(r8) :: weighted_dir_tr(nclmax) - real(r8) :: weighted_fsun(nclmax) - real(r8) :: weighted_dif_ratio(nclmax,maxSWb) - real(r8) :: weighted_dif_down(nclmax) - real(r8) :: weighted_dif_up(nclmax) - real(r8) :: refl_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation reflected by laye - real(r8) :: tran_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation transmitted by layer - real(r8) :: dif_ratio(nclmax,maxpft,nlevleaf,maxSWb) ! Ratio of upward to forward diffuse fluxes - real(r8) :: Dif_dn(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) - real(r8) :: Dif_up(nclmax,maxpft,nlevleaf) ! Upward diffuse flux above canopy layer J (W/m**2 ground area) - real(r8) :: lai_change(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) - real(r8) :: f_not_abs(maxpft,maxSWb) ! Fraction reflected + transmitted. 1-absorbtion. - real(r8) :: Abs_dir_z(maxpft,nlevleaf) - real(r8) :: Abs_dif_z(maxpft,nlevleaf) - real(r8) :: abs_rad(maxSWb) !radiation absorbed by soil - real(r8) :: tr_soili ! Radiation transmitted to the soil surface. - real(r8) :: tr_soild ! Radiation transmitted to the soil surface. - real(r8) :: phi1b(maxPatchesPerSite,maxpft) ! Radiation transmitted to the soil surface. - real(r8) :: phi2b(maxPatchesPerSite,maxpft) - real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer) - real(r8) :: angle - - real(r8),parameter :: tolerance = 0.000000001_r8 - real(r8), parameter :: pi = 3.141592654 ! PI - - - integer, parameter :: max_diag_nlevleaf = 4 - integer, parameter :: diag_nlevleaf = min(nlevleaf,max_diag_nlevleaf) ! for diagnostics, write a small number of leaf layers - - real(r8) :: denom - real(r8) :: lai_reduction(nclmax) - - integer :: fp,iv,s ! array indices - integer :: ib ! waveband number - real(r8) :: cosz ! 0.001 <= coszen <= 1.000 - real(r8) :: chil(maxPatchesPerSite) ! -0.4 <= xl <= 0.6 - real(r8) :: gdir(maxPatchesPerSite) ! leaf projection in solar direction (0 to 1) + integer :: s ! site loop counter + integer :: ifp ! patch loop counter + integer :: ib ! radiation broad band counter + type(ed_patch_type), pointer :: currentPatch ! patch pointer !----------------------------------------------------------------------- - - associate(& - rhol => EDPftvarcon_inst%rhol , & ! Input: [real(r8) (:) ] leaf reflectance: 1=vis, 2=nir - rhos => EDPftvarcon_inst%rhos , & ! Input: [real(r8) (:) ] stem reflectance: 1=vis, 2=nir - taul => EDPftvarcon_inst%taul , & ! Input: [real(r8) (:) ] leaf transmittance: 1=vis, 2=nir - taus => EDPftvarcon_inst%taus , & ! Input: [real(r8) (:) ] stem transmittance: 1=vis, 2=nir - xl => EDPftvarcon_inst%xl , & ! Input: [real(r8) (:) ] ecophys const - leaf/stem orientation index - clumping_index => EDPftvarcon_inst%clumping_index) ! Input: [real(r8) (:) ] ecophys const - leaf/stem clumping index - -! albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:) ] surface albedo (direct) (USED IN LND2ATM,BALANCE_CHECK) -! albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:) ] surface albedo (diffuse) (LND2ATM,BALANCE_CHECK) -! fabd => surfalb_inst%fabd_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit direct flux (BALANCE_CHECK) -! fabi => surfalb_inst%fabi_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit diffuse flux (BALANCE_CHECK) -! ftdd => surfalb_inst%ftdd_patch , & ! Output: [real(r8) (:,:) ] down direct flux below canopy per unit direct flx (BALANCE_CHECK) -! ftid => surfalb_inst%ftid_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit direct flx (BALANCE_CHECK) -! ftii => surfalb_inst%ftii_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flx (BALANCE_CHECK) - - ! ------------------------------------------------------------------------------- - ! TODO (mv, 2014-10-29) the filter here is different than below - ! this is needed to have the VOC's be bfb - this needs to be - ! re-examined int he future - ! RGK,2016-08-06: FATES is still incompatible with VOC emission module - ! ------------------------------------------------------------------------------- + ! ------------------------------------------------------------------------------- + ! TODO (mv, 2014-10-29) the filter here is different than below + ! this is needed to have the VOC's be bfb - this needs to be + ! re-examined int he future + ! RGK,2016-08-06: FATES is still incompatible with VOC emission module + ! ------------------------------------------------------------------------------- do s = 1, nsites @@ -167,24 +105,17 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) currentPatch%nrmlzd_parprof_dir_z(:,:,:) = 0._r8 currentPatch%nrmlzd_parprof_dif_z(:,:,:) = 0._r8 - if(bc_in(s)%filter_vegzen_pa(ifp))then + currentPatch%solar_zenith_flag = bc_in(s)%filter_vegzen_pa(ifp) + currentPatch%solar_zenith_angle = bc_in(s)%coszen_pa(ifp) + currentPatch%gnd_alb_dif(1:hlm_numSWb) = bc_in(s)%albgr_dif_rb(1:hlm_numSWb) + currentPatch%gnd_alb_dir(1:hlm_numSWb) = bc_in(s)%albgr_dir_rb(1:hlm_numSWb) + + if(currentPatch%solar_zenith_flag )then - weighted_dir_tr(:) = 0._r8 - weighted_dif_down(:) = 0._r8 - weighted_dif_up(:) = 0._r8 bc_out(s)%albd_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%albi_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%fabi_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%fabd_parb(ifp,:) = 0._r8 ! output HLM - tr_dir_z(:,:,:) = 0._r8 - tr_dif_z(:,:,:) = 0._r8 - ftweight(:,:,:) = 0._r8 - lai_change(:,:,:) = 0._r8 - Dif_up(:,:,:) = 0._r8 - Dif_dn(:,:,:) = 0._r8 - refl_dif(:,:,:,:) = 0.0_r8 - tran_dif(:,:,:,:) = 0.0_r8 - dif_ratio(:,:,:,:) = 0.0_r8 bc_out(s)%ftdd_parb(ifp,:) = 1._r8 ! output HLM bc_out(s)%ftid_parb(ifp,:) = 1._r8 ! output HLM bc_out(s)%ftii_parb(ifp,:) = 1._r8 ! output HLM @@ -196,755 +127,909 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) bc_out(s)%fabi_parb(ifp,:) = 0.0_r8 do ib = 1,hlm_numSWb bc_out(s)%albd_parb(ifp,ib) = bc_in(s)%albgr_dir_rb(ib) - bc_out(s)%albd_parb(ifp,ib) = bc_in(s)%albgr_dif_rb(ib) + bc_out(s)%albi_parb(ifp,ib) = bc_in(s)%albgr_dif_rb(ib) bc_out(s)%ftdd_parb(ifp,ib)= 1.0_r8 bc_out(s)%ftid_parb(ifp,ib)= 1.0_r8 bc_out(s)%ftii_parb(ifp,ib)= 1.0_r8 enddo + else + + call PatchNormanRadiation (currentPatch, & + bc_out(s)%albd_parb(ifp,:), & + bc_out(s)%albi_parb(ifp,:), & + bc_out(s)%fabd_parb(ifp,:), & + bc_out(s)%fabi_parb(ifp,:), & + bc_out(s)%ftdd_parb(ifp,:), & + bc_out(s)%ftid_parb(ifp,:), & + bc_out(s)%ftii_parb(ifp,:)) - ! Is this pft/canopy layer combination present in this patch? - do L = 1,nclmax - do ft = 1,numpft - currentPatch%canopy_mask(L,ft) = 0 - do iv = 1, currentPatch%nrad(L,ft) - if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then - currentPatch%canopy_mask(L,ft) = 1 - !I think 'present' is only used here... - endif - end do !iv - end do !ft - end do !L - - !do this once for one unit of diffuse, and once for one unit of direct radiation - do radtype = 1, n_rad_stream_types - do ib = 1,hlm_numSWb - if (radtype == idirect) then - ! Set the hypothetical driving radiation. We do this once for a single unit of direct and - ! once for a single unit of diffuse radiation. - forc_dir(ifp,ib) = 1.00_r8 - forc_dif(ifp,ib) = 0.00_r8 - else !dif - forc_dir(ifp,ib) = 0.00_r8 - forc_dif(ifp,ib) = 1.00_r8 - end if - end do !ib - - !Extract information that needs to be provided by ED into local array. - ftweight(:,:,:) = 0._r8 - do L = 1,currentPatch%NCL_p - do ft = 1,numpft - do iv = 1, currentPatch%nrad(L,ft) - !this is already corrected for area in CLAP - ftweight(L,ft,iv) = currentPatch%canopy_area_profile(L,ft,iv) - end do !iv - end do !ft1 - end do !L - if (sum(ftweight(1,:,1))<0.999_r8)then - write(fates_log(),*) 'canopy not full',ftweight(1,:,1) - endif - if (sum(ftweight(1,:,1))>1.0001_r8)then - write(fates_log(),*) 'canopy too full',ftweight(1,:,1) - endif - - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Direct beam extinction coefficient, k_dir. PFT specific. - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - cosz = max(0.001_r8, bc_in(s)%coszen_pa(ifp)) !copied from previous radiation code... - do ft = 1,numpft - sb = (90._r8 - (acos(cosz)*180/pi)) * (pi / 180._r8) - chil(ifp) = xl(ft) !min(max(xl(ft), -0.4_r8), 0.6_r8 ) - if (abs(chil(ifp)) <= 0.01_r8) then - chil(ifp) = 0.01_r8 - end if - phi1b(ifp,ft) = 0.5_r8 - 0.633_r8*chil(ifp) - 0.330_r8*chil(ifp)*chil(ifp) - phi2b(ifp,ft) = 0.877_r8 * (1._r8 - 2._r8*phi1b(ifp,ft)) !0 = horiz leaves, 1 - vert leaves. - gdir(ifp) = phi1b(ifp,ft) + phi2b(ifp,ft) * sin(sb) - !how much direct light penetrates a singleunit of lai? - k_dir(ft) = clumping_index(ft) * gdir(ifp) / sin(sb) - end do !FT - - do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.) - weighted_dir_tr(L) = 0.0_r8 - weighted_fsun(L) = 0._r8 - weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8 - !Each canopy layer (canopy, understorey) has multiple 'parallel' pft's - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves. - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. - ! Estimated do nine sky angles in increments of 10 degrees - ! PFT specific... - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - tr_dif_z(L,ft,:) = 0._r8 - do iv = 1,currentPatch%nrad(L,ft) - do j = 1,9 - angle = (5._r8 + (j - 1) * 10._r8) * 3.142 / 180._r8 - gdir(ifp) = phi1b(ifp,ft) + phi2b(ifp,ft) * sin(angle) - tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) + exp(-clumping_index(ft) * & - gdir(ifp) / sin(angle) * & - (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))) * & - sin(angle)*cos(angle) - end do - - tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) * 2._r8 * (10.00*pi/180._r8) - - end do - - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Direct beam transmittance, tr_dir_z, uses cumulative LAI above layer J to give - ! unscattered direct beam onto layer J. do each PFT section. - ! This is just an decay curve based on k_dir. (leaf & sun angle) - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - if (L==1)then - tr_dir_z(L,ft,1) = 1._r8 - else - tr_dir_z(L,ft,1) = weighted_dir_tr(L-1) - endif - laisum = 0.00_r8 - !total direct beam getting to the bottom of the top canopy. - do iv = 1,currentPatch%nrad(L,ft) - laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) - lai_change(L,ft,iv) = 0.0_r8 - if (( ftweight(L,ft,iv+1) > 0.0_r8 ) .and. ( ftweight(L,ft,iv+1) < ftweight(L,ft,iv) ))then - !where there is a partly empty leaf layer, some fluxes go straight through. - lai_change(L,ft,iv) = ftweight(L,ft,iv)-ftweight(L,ft,iv+1) - endif - if (ftweight(L,ft,iv+1) - ftweight(L,ft,iv) > 1.e-10_r8)then - write(fates_log(),*) 'lower layer has more coverage. This is wrong' , & - ftweight(L,ft,iv),ftweight(L,ft,iv+1),ftweight(L,ft,iv+1)-ftweight(L,ft,iv) - endif - - !n.b. in theory lai_change could be calculated daily in the ED code. - !This is light coming striaght through the canopy. - if (L==1)then - tr_dir_z(L,ft,iv+1) = exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - else - tr_dir_z(L,ft,iv+1) = weighted_dir_tr(L-1)*exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - endif - - if (iv == 1)then - !this is the top layer. - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & - ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) - else - !the lai_change(iv) affects the light incident on layer iv+2 not iv+1 - ! light coming from the layer above (iv-1) goes through iv and onto iv+1. - if (lai_change(L,ft,iv-1) > 0.0_r8)then - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv)* & - lai_change(L,ft,iv-1) / ftweight(L,ft,1) - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv-1)* & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) - else - !account fot the light that comes striaght down from unfilled layers above. - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & - ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) - endif - endif - end do - - !add up all the weighted contributions from the different PFT columns. - weighted_dir_tr(L) = weighted_dir_tr(L) + tr_dir_z(L,ft,currentPatch%nrad(L,ft)+1)*ftweight(L,ft,1) + + endif ! is there vegetation? + + end if ! if the vegetation and zenith filter is active + currentPatch => currentPatch%younger + end do ! Loop linked-list patches + enddo ! Loop Sites + + return + end subroutine ED_Norman_Radiation + - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Sunlit and shaded fraction of leaf layer - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - - !laisum = 0._r8 - do iv = 1,currentPatch%nrad(L,ft) - ! Cumulative leaf area. Original code uses cumulative lai do layer. - ! Now use cumulative lai at center of layer. - ! Same as tr_dir_z calcualtions, but in the middle of the layer? FIX(RF,032414)-WHY? - if (iv == 1) then - laisum = 0.5_r8 * (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)) - else - laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) - end if - - - if (L == 1)then !top canopy layer - currentPatch%f_sun(L,ft,iv) = exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - else - currentPatch%f_sun(L,ft,iv) = weighted_fsun(L-1)* exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - endif - - if ( iv > 1 ) then ! becasue we are looking at this layer (not the next) - ! we only ever add fluxes if iv>1 - if (lai_change(L,ft,iv-1) > 0.0_r8)then - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv) * & - lai_change(L,ft,iv-1)/ftweight(L,ft,1) - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv-1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) - else - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv-1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - endif - - end do !iv - weighted_fsun(L) = weighted_fsun(L) + currentPatch%f_sun(L,ft,currentPatch%nrad(L,ft))* & - ftweight(L,ft,1) - - ! instance where the first layer ftweight is used a proxy for the whole column. FTWA - ! this is possibly a source of slight error. If we use the ftweight at the top of the PFT column, - ! then we willl underestimate fsun, but if we use ftweight at the bottom of the column, we will - ! underestimate it. Really, we should be tracking the release of direct light from the column as it tapers - ! towards the ground. Is that necessary to get energy closure? It would be quite hard... - endif !present. - end do!pft loop - end do !L - - do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. - do ft = 1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !==============================================================================! - ! Iterative solution do scattering - !==============================================================================! + ! ====================================================================================== + + subroutine PatchNormanRadiation (currentPatch, & + albd_parb_out, & ! (ifp,ib) + albi_parb_out, & ! (ifp,ib) + fabd_parb_out, & ! (ifp,ib) + fabi_parb_out, & ! (ifp,ib) + ftdd_parb_out, & ! (ifp,ib) + ftid_parb_out, & ! (ifp,ib) + ftii_parb_out) ! (ifp,ib) + + ! ----------------------------------------------------------------------------------- + ! + ! This routine performs the Norman Radiation scattering for each patch. + ! + ! ----------------------------------------------------------------------------------- + + ! + ! !USES: + use EDPftvarcon , only : EDPftvarcon_inst + use EDtypesMod , only : ed_patch_type + + ! ----------------------------------------------------------------------------------- + ! !ARGUMENTS: + ! ----------------------------------------------------------------------------------- + + type(ed_patch_type), intent(inout), target :: currentPatch + real(r8), intent(inout) :: albd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: albi_parb_out(hlm_numSWb) + real(r8), intent(inout) :: fabd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: fabi_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftdd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftid_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftii_parb_out(hlm_numSWb) + + ! Locals + ! ----------------------------------------------------------------------------------- + + integer :: radtype, L, ft, j + integer :: iter ! Iteration index + integer :: irep ! Flag to exit iteration loop + real(r8) :: sb + real(r8) :: error ! Error check + real(r8) :: down_rad, up_rad ! Iterative solution do Dif_dn and Dif_up + real(r8) :: ftweight(nclmax,maxpft,nlevleaf) + real(r8) :: k_dir(maxpft) ! Direct beam extinction coefficient + real(r8) :: tr_dir_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of direct beam radiation through a single layer + real(r8) :: tr_dif_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of diffuse radiation through a single layer + real(r8) :: weighted_dir_tr(nclmax) + real(r8) :: weighted_fsun(nclmax) + real(r8) :: weighted_dif_ratio(nclmax,maxSWb) + real(r8) :: weighted_dif_down(nclmax) + real(r8) :: weighted_dif_up(nclmax) + real(r8) :: refl_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation reflected by laye + real(r8) :: tran_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation transmitted by layer + real(r8) :: dif_ratio(nclmax,maxpft,nlevleaf,maxSWb) ! Ratio of upward to forward diffuse fluxes + real(r8) :: Dif_dn(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) + real(r8) :: Dif_up(nclmax,maxpft,nlevleaf) ! Upward diffuse flux above canopy layer J (W/m**2 ground area) + real(r8) :: lai_change(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) + real(r8) :: f_not_abs(maxpft,maxSWb) ! Fraction reflected + transmitted. 1-absorbtion. + real(r8) :: Abs_dir_z(maxpft,nlevleaf) + real(r8) :: Abs_dif_z(maxpft,nlevleaf) + real(r8) :: abs_rad(maxSWb) !radiation absorbed by soil + real(r8) :: tr_soili ! Radiation transmitted to the soil surface. + real(r8) :: tr_soild ! Radiation transmitted to the soil surface. + real(r8) :: phi1b(maxpft) ! Radiation transmitted to the soil surface. + real(r8) :: phi2b(maxpft) + real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer) + real(r8) :: angle + + real(r8),parameter :: tolerance = 0.000000001_r8 + + + integer, parameter :: max_diag_nlevleaf = 4 + integer, parameter :: diag_nlevleaf = min(nlevleaf,max_diag_nlevleaf) ! for diagnostics, write a small number of leaf layers + + real(r8) :: denom + real(r8) :: lai_reduction(nclmax) + + integer :: fp,iv,s ! array indices + integer :: ib ! waveband number + real(r8) :: cosz ! 0.001 <= coszen <= 1.000 + real(r8) :: chil + real(r8) :: gdir + + + real(r8), parameter :: forc_dir(n_rad_stream_types) = (/ 1.0_r8, 0.0_r8 /) ! These are binary switches used + real(r8), parameter :: forc_dif(n_rad_stream_types) = (/ 0.0_r8, 1.0_r8 /) ! to turn off and on radiation streams + + + + associate(& + rhol => EDPftvarcon_inst%rhol , & ! Input: [real(r8) (:) ] leaf reflectance: 1=vis, 2=nir + rhos => EDPftvarcon_inst%rhos , & ! Input: [real(r8) (:) ] stem reflectance: 1=vis, 2=nir + taul => EDPftvarcon_inst%taul , & ! Input: [real(r8) (:) ] leaf transmittance: 1=vis, 2=nir + taus => EDPftvarcon_inst%taus , & ! Input: [real(r8) (:) ] stem transmittance: 1=vis, 2=nir + xl => EDPftvarcon_inst%xl , & ! Input: [real(r8) (:) ] ecophys const - leaf/stem orientation index + clumping_index => EDPftvarcon_inst%clumping_index) + + + ! Initialize local arrays + + weighted_dir_tr(:) = 0._r8 + weighted_dif_down(:) = 0._r8 + weighted_dif_up(:) = 0._r8 + + tr_dir_z(:,:,:) = 0._r8 + tr_dif_z(:,:,:) = 0._r8 + lai_change(:,:,:) = 0._r8 + Dif_up(:,:,:) = 0._r8 + Dif_dn(:,:,:) = 0._r8 + refl_dif(:,:,:,:) = 0.0_r8 + tran_dif(:,:,:,:) = 0.0_r8 + dif_ratio(:,:,:,:) = 0.0_r8 + + + ! Initialize the ouput arrays + ! --------------------------------------------------------------------------------- + albd_parb_out(1:hlm_numSWb) = 0.0_r8 + albi_parb_out(1:hlm_numSWb) = 0.0_r8 + fabd_parb_out(1:hlm_numSWb) = 0.0_r8 + fabi_parb_out(1:hlm_numSWb) = 0.0_r8 + ftdd_parb_out(1:hlm_numSWb) = 1.0_r8 + ftid_parb_out(1:hlm_numSWb) = 1.0_r8 + ftii_parb_out(1:hlm_numSWb) = 1.0_r8 + + ! Is this pft/canopy layer combination present in this patch? + + do L = 1,nclmax + do ft = 1,numpft + currentPatch%canopy_mask(L,ft) = 0 + do iv = 1, currentPatch%nrad(L,ft) + if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then + currentPatch%canopy_mask(L,ft) = 1 + !I think 'present' is only used here... + endif + end do !iv + end do !ft + end do !L + + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Direct beam extinction coefficient, k_dir. PFT specific. + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + cosz = max(0.001_r8, currentPatch%solar_zenith_angle ) !copied from previous radiation code... + do ft = 1,numpft + sb = (90._r8 - (acos(cosz)*180._r8/pi_const)) * (pi_const / 180._r8) + chil = xl(ft) !min(max(xl(ft), -0.4_r8), 0.6_r8 ) + if ( abs(chil) <= 0.01_r8) then + chil = 0.01_r8 + end if + phi1b(ft) = 0.5_r8 - 0.633_r8*chil - 0.330_r8*chil*chil + phi2b(ft) = 0.877_r8 * (1._r8 - 2._r8*phi1b(ft)) !0 = horiz leaves, 1 - vert leaves. + gdir = phi1b(ft) + phi2b(ft) * sin(sb) + !how much direct light penetrates a singleunit of lai? + k_dir(ft) = clumping_index(ft) * gdir / sin(sb) + end do !FT + + + + + !do this once for one unit of diffuse, and once for one unit of direct radiation + do radtype = 1, n_rad_stream_types + + ! Extract information that needs to be provided by ED into local array. + ! RGK: NOT SURE WHY WE NEED FTWEIGHT ... + ! ------------------------------------------------------------------------------ + + ftweight(:,:,:) = 0._r8 + do L = 1,currentPatch%NCL_p + do ft = 1,numpft + do iv = 1, currentPatch%nrad(L,ft) + !this is already corrected for area in CLAP + ftweight(L,ft,iv) = currentPatch%canopy_area_profile(L,ft,iv) + end do !iv + end do !ft1 + end do !L + if (sum(ftweight(1,:,1))<0.999_r8)then + write(fates_log(),*) 'canopy not full',ftweight(1,:,1) + endif + if (sum(ftweight(1,:,1))>1.0001_r8)then + write(fates_log(),*) 'canopy too full',ftweight(1,:,1) + endif + + do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.) + + weighted_dir_tr(L) = 0.0_r8 + weighted_fsun(L) = 0._r8 + weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8 + + !Each canopy layer (canopy, understorey) has multiple 'parallel' pft's + + do ft =1,numpft + + if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves. + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. + ! Estimated do nine sky angles in increments of 10 degrees + ! PFT specific... + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + tr_dif_z(L,ft,:) = 0._r8 + do iv = 1,currentPatch%nrad(L,ft) + do j = 1,9 + angle = (5._r8 + real(j - 1,r8) * 10._r8) * pi_const / 180._r8 + gdir = phi1b(ft) + phi2b(ft) * sin(angle) + tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) + exp(-clumping_index(ft) * & + gdir / sin(angle) * & + (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))) * & + sin(angle)*cos(angle) + end do + + tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) * 2._r8 * (10._r8 * pi_const / 180._r8) + + end do + + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Direct beam transmittance, tr_dir_z, uses cumulative LAI above layer J to give + ! unscattered direct beam onto layer J. do each PFT section. + ! This is just an decay curve based on k_dir. (leaf & sun angle) + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + if (L==1)then + tr_dir_z(L,ft,1) = 1._r8 + else + tr_dir_z(L,ft,1) = weighted_dir_tr(L-1) + endif + laisum = 0.00_r8 + !total direct beam getting to the bottom of the top canopy. + do iv = 1,currentPatch%nrad(L,ft) + laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) + lai_change(L,ft,iv) = 0.0_r8 + if (( ftweight(L,ft,iv+1) > 0.0_r8 ) .and. ( ftweight(L,ft,iv+1) < ftweight(L,ft,iv) ))then + !where there is a partly empty leaf layer, some fluxes go straight through. + lai_change(L,ft,iv) = ftweight(L,ft,iv)-ftweight(L,ft,iv+1) + endif + if (ftweight(L,ft,iv+1) - ftweight(L,ft,iv) > 1.e-10_r8)then + write(fates_log(),*) 'lower layer has more coverage. This is wrong' , & + ftweight(L,ft,iv),ftweight(L,ft,iv+1),ftweight(L,ft,iv+1)-ftweight(L,ft,iv) + endif + + !n.b. in theory lai_change could be calculated daily in the ED code. + !This is light coming striaght through the canopy. + if (L==1)then + tr_dir_z(L,ft,iv+1) = exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + else + tr_dir_z(L,ft,iv+1) = weighted_dir_tr(L-1)*exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + endif + + if (iv == 1)then + !this is the top layer. + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & + ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) + else + !the lai_change(iv) affects the light incident on layer iv+2 not iv+1 + ! light coming from the layer above (iv-1) goes through iv and onto iv+1. + if (lai_change(L,ft,iv-1) > 0.0_r8)then + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv)* & + lai_change(L,ft,iv-1) / ftweight(L,ft,1) + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv-1)* & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) + else + !account fot the light that comes striaght down from unfilled layers above. + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & + ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) + endif + endif + + end do + + !add up all the weighted contributions from the different PFT columns. + weighted_dir_tr(L) = weighted_dir_tr(L) + tr_dir_z(L,ft,currentPatch%nrad(L,ft)+1)*ftweight(L,ft,1) + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Sunlit and shaded fraction of leaf layer + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - do ib = 1,hlm_numSWb !vis, nir - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Leaf scattering coefficient and terms do diffuse radiation reflected - ! and transmitted by a layer - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - f_not_abs(ft,ib) = rhol(ft,ib) + taul(ft,ib) !leaf level fraction NOT absorbed. - !tr_dif_z is a term that uses the LAI in each layer, whereas rhol and taul do not, - !because they are properties of leaf surfaces and not of the leaf matrix. - do iv = 1,currentPatch%nrad(L,ft) - !How much diffuse light is intercepted and then reflected? - refl_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * rhol(ft,ib) - !How much diffuse light in this layer is transmitted? - tran_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * & - taul(ft,ib) + tr_dif_z(L,ft,iv) - end do + !laisum = 0._r8 + do iv = 1,currentPatch%nrad(L,ft) + ! Cumulative leaf area. Original code uses cumulative lai do layer. + ! Now use cumulative lai at center of layer. + ! Same as tr_dir_z calcualtions, but in the middle of the layer? FIX(RF,032414)-WHY? + if (iv == 1) then + laisum = 0.5_r8 * (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)) + else + laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) + end if + + + if (L == 1)then !top canopy layer + currentPatch%f_sun(L,ft,iv) = exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + else + currentPatch%f_sun(L,ft,iv) = weighted_fsun(L-1)* exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + endif + + if ( iv > 1 ) then ! becasue we are looking at this layer (not the next) + ! we only ever add fluxes if iv>1 + if (lai_change(L,ft,iv-1) > 0.0_r8)then + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv) * & + lai_change(L,ft,iv-1)/ftweight(L,ft,1) + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv-1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) + else + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv-1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + endif - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Ratio of upward to forward diffuse fluxes, dif_ratio - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Soil diffuse reflectance (ratio of down to up radiation). - iv = currentPatch%nrad(L,ft) + 1 - if (L == currentPatch%NCL_p)then !nearest the soil - dif_ratio(L,ft,iv,ib) = bc_in(s)%albgr_dif_rb(ib) - else - dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib) - end if - ! Canopy layers, working upwardfrom soil with dif_ratio(iv+1) known - ! FIX(RF,032414) ray tracing eqution - need to find derivation of this... - ! for each unit going down, there are x units going up. - do iv = currentPatch%nrad(L,ft),1, -1 - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv+1,ib) * & - tran_dif(L,ft,iv,ib)*tran_dif(L,ft,iv,ib) / & - (1._r8 - dif_ratio(L,ft,iv+1,ib) * refl_dif(L,ft,iv,ib)) & - + refl_dif(L,ft,iv,ib) - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) * & - ftweight(L,ft,iv)/ftweight(L,ft,1) - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) + dif_ratio(L,ft,iv+1,ib) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - end do - weighted_dif_ratio(L,ib) = weighted_dif_ratio(L,ib) + & - dif_ratio(L,ft,1,ib) * ftweight(L,ft,1) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - end do!hlm_numSWb - endif ! currentPatch%canopy_mask - end do!ft - end do!L + end do !iv + + weighted_fsun(L) = weighted_fsun(L) + currentPatch%f_sun(L,ft,currentPatch%nrad(L,ft))* & + ftweight(L,ft,1) + + ! instance where the first layer ftweight is used a proxy for the whole column. FTWA + ! this is possibly a source of slight error. If we use the ftweight at the top of the PFT column, + ! then we willl underestimate fsun, but if we use ftweight at the bottom of the column, we will + ! underestimate it. Really, we should be tracking the release of direct light from the column as it tapers + ! towards the ground. Is that necessary to get energy closure? It would be quite hard... + endif !present. + end do!pft loop + end do !L + + + do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. + do ft = 1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + + !==============================================================================! + ! Iterative solution do scattering + !==============================================================================! + + do ib = 1,hlm_numSWb !vis, nir + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Leaf scattering coefficient and terms do diffuse radiation reflected + ! and transmitted by a layer + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + f_not_abs(ft,ib) = rhol(ft,ib) + taul(ft,ib) !leaf level fraction NOT absorbed. + !tr_dif_z is a term that uses the LAI in each layer, whereas rhol and taul do not, + !because they are properties of leaf surfaces and not of the leaf matrix. + do iv = 1,currentPatch%nrad(L,ft) + !How much diffuse light is intercepted and then reflected? + refl_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * rhol(ft,ib) + !How much diffuse light in this layer is transmitted? + tran_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * & + taul(ft,ib) + tr_dif_z(L,ft,iv) + end do + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Ratio of upward to forward diffuse fluxes, dif_ratio + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Soil diffuse reflectance (ratio of down to up radiation). + iv = currentPatch%nrad(L,ft) + 1 + if (L == currentPatch%NCL_p)then !nearest the soil + dif_ratio(L,ft,iv,ib) = currentPatch%gnd_alb_dif(ib) !bc_in(s)%albgr_dif_rb(ib) + else + dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib) + end if + ! Canopy layers, working upwardfrom soil with dif_ratio(iv+1) known + ! FIX(RF,032414) ray tracing eqution - need to find derivation of this... + ! for each unit going down, there are x units going up. + do iv = currentPatch%nrad(L,ft),1, -1 + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv+1,ib) * & + tran_dif(L,ft,iv,ib)*tran_dif(L,ft,iv,ib) / & + (1._r8 - dif_ratio(L,ft,iv+1,ib) * refl_dif(L,ft,iv,ib)) & + + refl_dif(L,ft,iv,ib) + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) * & + ftweight(L,ft,iv)/ftweight(L,ft,1) + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) + dif_ratio(L,ft,iv+1,ib) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + end do + weighted_dif_ratio(L,ib) = weighted_dif_ratio(L,ib) + & + dif_ratio(L,ft,1,ib) * ftweight(L,ft,1) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + end do!hlm_numSWb + endif ! currentPatch%canopy_mask + end do!ft + end do!L - do ib = 1,hlm_numSWb - Dif_dn(:,:,:) = 0.00_r8 - Dif_up(:,:,:) = 0.00_r8 - do L = 1, currentPatch%NCL_p !work down from the top of the canopy. - weighted_dif_down(L) = 0._r8 - do ft = 1, numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! First estimates do downward and upward diffuse flux - ! - ! Dif_dn = forward diffuse flux onto layer J - ! Dif_up = Upward diffuse flux above layer J - ! - ! Solved here without direct beam radiation and using dif_ratio = Dif_up / Dif_dn - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! downward diffuse flux onto the top surface of the canopy - - if (L == 1)then - Dif_dn(L,ft,1) = forc_dif(ifp,ib) - else - Dif_dn(L,ft,1) = weighted_dif_down(L-1) - end if - ! forward diffuse flux within the canopy and at soil, working forward through canopy - do iv = 1,currentPatch%nrad(L,ft) - denom = refl_dif(L,ft,iv,ib) * dif_ratio(L,ft,iv,ib) - denom = 1._r8 - denom - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) / & - denom *ftweight(L,ft,iv)/ftweight(L,ft,1) - if (iv > 1)then - if (lai_change(L,ft,iv-1) > 0.0_r8)then - !here we are thinking about whether the layer above had an laichange, - !but calculating the flux onto the layer below. - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv)* & - lai_change(L,ft,iv-1)/ftweight(L,ft,1) - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv-1)* & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1)/ftweight(L,ft,1)) - else - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - else - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - end do - - weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & - ftweight(L,ft,1) - - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif !present - end do !ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey? - !Add on the radiation going through the canopy gaps. - weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif - end do !L - - do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. - weighted_dif_up(L) = 0._r8 - do ft = 1, numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !Bounce diffuse radiation off soil surface. - iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then !is this the bottom layer ? - Dif_up(L,ft,iv) =bc_in(s)%albgr_dif_rb(ib) * Dif_dn(L,ft,iv) - else - Dif_up(L,ft,iv) = weighted_dif_up(L+1) - end if - ! Upward diffuse flux within the canopy and above the canopy, working upward through canopy - - do iv = currentPatch%nrad(L,ft), 1, -1 - if (lai_change(L,ft,iv) > 0.0_r8)then - Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * & - ftweight(L,ft,iv) / ftweight(L,ft,1) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & - tran_dif(L,ft,iv,ib) * lai_change(L,ft,iv)/ftweight(L,ft,1) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - !nb is this the right constuction? - ! the radiation that hits the empty space is not reflected. - else - Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * ftweight(L,ft,iv) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * (1.0_r8-ftweight(L,ft,iv)) - endif - end do + + do ib = 1,hlm_numSWb + Dif_dn(:,:,:) = 0.00_r8 + Dif_up(:,:,:) = 0.00_r8 + do L = 1, currentPatch%NCL_p !work down from the top of the canopy. + weighted_dif_down(L) = 0._r8 + do ft = 1, numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! First estimates do downward and upward diffuse flux + ! + ! Dif_dn = forward diffuse flux onto layer J + ! Dif_up = Upward diffuse flux above layer J + ! + ! Solved here without direct beam radiation and using dif_ratio = Dif_up / Dif_dn + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! downward diffuse flux onto the top surface of the canopy + + if (L == 1)then + Dif_dn(L,ft,1) = forc_dif(radtype) + else + Dif_dn(L,ft,1) = weighted_dif_down(L-1) + end if + ! forward diffuse flux within the canopy and at soil, working forward through canopy + do iv = 1,currentPatch%nrad(L,ft) + denom = refl_dif(L,ft,iv,ib) * dif_ratio(L,ft,iv,ib) + denom = 1._r8 - denom + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) / & + denom *ftweight(L,ft,iv)/ftweight(L,ft,1) + if (iv > 1)then + if (lai_change(L,ft,iv-1) > 0.0_r8)then + !here we are thinking about whether the layer above had an laichange, + !but calculating the flux onto the layer below. + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv)* & + lai_change(L,ft,iv-1)/ftweight(L,ft,1) + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv-1)* & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1)/ftweight(L,ft,1)) + else + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + else + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + end do + + weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & + ftweight(L,ft,1) - weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif !present - end do !ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? - !Add on the radiation coming up through the canopy gaps. - !diffuse to diffuse - weighted_dif_up(L) = weighted_dif_up(L) +(1.0-sum(ftweight(L,1:numpft,1))) * & - weighted_dif_down(L-1) * bc_in(s)%albgr_dif_rb(ib) - !direct to diffuse - weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(ifp,ib) * & - weighted_dir_tr(L-1) * (1.0-sum(ftweight(L,1:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) - endif - end do !L - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! 3. Iterative calculation of forward and upward diffuse fluxes, iNCL_puding - ! scattered direct beam - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - - ! Flag to exit iteration loop: 0 = exit and 1 = iterate - irep = 1 - ! Iteration loop - iter = 0 - do while(irep ==1 .and. iter<50) - - iter = iter + 1 - irep = 0 - do L = 1,currentPatch%NCL_p !working from the top down - weighted_dif_down(L) = 0._r8 - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - ! forward diffuse flux within the canopy and at soil, working forward through canopy - ! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy. - ! Note: down = forward flux onto next layer - if (L == 1)then !is this the top layer? - Dif_dn(L,ft,1) = forc_dif(ifp,ib) - else - Dif_dn(L,ft,1) = weighted_dif_down(L-1) - end if - down_rad = 0._r8 - - do iv = 1, currentPatch%nrad(L,ft) - - down_rad = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) + & - Dif_up(L,ft,iv+1) * refl_dif(L,ft,iv,ib) + & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * (1.00_r8 - & - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & - currentPatch%esai_profile(L,ft,iv)))) * taul(ft,ib) - down_rad = down_rad *(ftweight(L,ft,iv)/ftweight(L,ft,1)) - - if (iv > 1)then - if (lai_change(L,ft,iv-1) > 0.0_r8)then - down_rad = down_rad + Dif_dn(L,ft,iv) * lai_change(L,ft,iv-1)/ftweight(L,ft,1) - down_rad = down_rad + Dif_dn(L,ft,iv-1) * (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ & - ftweight(L,ft,1) - else - down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & - ftweight(L,ft,1) - endif - else - down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & - ftweight(L,ft,1) - endif - - !this is just Dif down, plus refl up, plus dir intercepted and turned into dif... , - if (abs(down_rad - Dif_dn(L,ft,iv+1)) > tolerance)then - irep = 1 - end if - Dif_dn(L,ft,iv+1) = down_rad - - end do !iv - - weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & - ftweight(L,ft,1) - - endif !present - end do!ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? - weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1) * & - (1.0-sum(ftweight(L,1:numpft,1))) - end if - end do ! do L loop - - do L = 1, currentPatch%NCL_p ! working from the top down. - weighted_dif_up(L) = 0._r8 - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam) - iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil - Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) *bc_in(s)%albgr_dif_rb(ib) + & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) *bc_in(s)%albgr_dir_rb(ib) - else !In the other canopy layers, reflect off the underlying vegetation. - Dif_up(L,ft,iv) = weighted_dif_up(L+1) - end if - - ! Upward diffuse flux within and above the canopy, working upward through canopy - ! with Dif_dn from previous interation. Note: up = upward flux above current layer - do iv = currentPatch%nrad(L,ft),1,-1 - !this is radiation up, by layer transmittance, by - - !reflection of the lower layer, - up_rad = Dif_dn(L,ft,iv) * refl_dif(L,ft,iv,ib) - up_rad = up_rad + forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * (1.00_r8 - exp(-k_dir(ft) * & - (currentPatch%elai_profile(L,ft,iv) + currentPatch%esai_profile(L,ft,iv)))) * & - rhol(ft,ib) - up_rad = up_rad + Dif_up(L,ft,iv+1) * tran_dif(L,ft,iv,ib) - up_rad = up_rad * ftweight(L,ft,iv)/ftweight(L,ft,1) - up_rad = up_rad + Dif_up(L,ft,iv+1) *(ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - ! THE LOWER LAYER FLUX IS HOMOGENIZED, SO WE DON"T CONSIDER THE LAI_CHANGE HERE... + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif !present + end do !ft + if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey? + !Add on the radiation going through the canopy gaps. + weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif + end do !L + + do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. + weighted_dif_up(L) = 0._r8 + do ft = 1, numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + !Bounce diffuse radiation off soil surface. + iv = currentPatch%nrad(L,ft) + 1 + if (L==currentPatch%NCL_p)then !is this the bottom layer ? + Dif_up(L,ft,iv) = currentPatch%gnd_alb_dif(ib) * Dif_dn(L,ft,iv) + else + Dif_up(L,ft,iv) = weighted_dif_up(L+1) + end if + ! Upward diffuse flux within the canopy and above the canopy, working upward through canopy + + do iv = currentPatch%nrad(L,ft), 1, -1 + if (lai_change(L,ft,iv) > 0.0_r8)then + Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * & + ftweight(L,ft,iv) / ftweight(L,ft,1) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & + tran_dif(L,ft,iv,ib) * lai_change(L,ft,iv)/ftweight(L,ft,1) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + !nb is this the right constuction? + ! the radiation that hits the empty space is not reflected. + else + Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * ftweight(L,ft,iv) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * (1.0_r8-ftweight(L,ft,iv)) + endif + end do + + weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif !present + end do !ft + if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + !Add on the radiation coming up through the canopy gaps. + !diffuse to diffuse + weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & + weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) + !direct to diffuse + weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * & + weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) * currentPatch%gnd_alb_dir(ib) + endif + end do !L + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! 3. Iterative calculation of forward and upward diffuse fluxes, iNCL_puding + ! scattered direct beam + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + + ! Flag to exit iteration loop: 0 = exit and 1 = iterate + irep = 1 + ! Iteration loop + iter = 0 + do while(irep ==1 .and. iter<50) + + iter = iter + 1 + irep = 0 + do L = 1,currentPatch%NCL_p !working from the top down + weighted_dif_down(L) = 0._r8 + do ft =1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + ! forward diffuse flux within the canopy and at soil, working forward through canopy + ! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy. + ! Note: down = forward flux onto next layer + if (L == 1)then !is this the top layer? + Dif_dn(L,ft,1) = forc_dif(radtype) + else + Dif_dn(L,ft,1) = weighted_dif_down(L-1) + end if + down_rad = 0._r8 + + do iv = 1, currentPatch%nrad(L,ft) + + down_rad = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) + & + Dif_up(L,ft,iv+1) * refl_dif(L,ft,iv,ib) + & + forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - & + exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & + currentPatch%esai_profile(L,ft,iv)))) * taul(ft,ib) + down_rad = down_rad *(ftweight(L,ft,iv)/ftweight(L,ft,1)) - if (abs(up_rad - Dif_up(L,ft,iv)) > tolerance) then !are we close to the tolerance level? - irep = 1 - end if - Dif_up(L,ft,iv) = up_rad + if (iv > 1)then + if (lai_change(L,ft,iv-1) > 0.0_r8)then + down_rad = down_rad + Dif_dn(L,ft,iv) * lai_change(L,ft,iv-1)/ftweight(L,ft,1) + down_rad = down_rad + Dif_dn(L,ft,iv-1) * (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ & + ftweight(L,ft,1) + else + down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & + ftweight(L,ft,1) + endif + else + down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & + ftweight(L,ft,1) + endif - end do !iv - weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) - end if !present - end do!ft - - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? - !Add on the radiation coming up through the canopy gaps. - weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & - weighted_dif_down(L-1) * bc_in(s)%albgr_dif_rb(ib) - weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(ifp,ib) * & - weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) - end if - end do!L - end do ! do while over iter - - abs_rad(ib) = 0._r8 - tr_soili = 0._r8 - tr_soild = 0._r8 - - do L = 1, currentPatch%NCL_p !working from the top down. - abs_dir_z(:,:) = 0._r8 - abs_dif_z(:,:) = 0._r8 - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !==============================================================================! - ! Compute absorbed flux densities - !==============================================================================! - - ! Absorbed direct beam and diffuse do leaf layers - do iv = 1, currentPatch%nrad(L,ft) - Abs_dir_z(ft,iv) = ftweight(L,ft,iv)* forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * & - (1.00_r8 - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & - currentPatch%esai_profile(L,ft,iv)))) * (1.00_r8 - f_not_abs(ft,ib)) - Abs_dif_z(ft,iv) = ftweight(L,ft,iv)* ((Dif_dn(L,ft,iv) + & - Dif_up(L,ft,iv+1)) * (1.00_r8 - tr_dif_z(L,ft,iv)) * & - (1.00_r8 - f_not_abs(ft,ib))) - end do - - ! Absorbed direct beam and diffuse do soil - if (L == currentPatch%NCL_p)then - iv = currentPatch%nrad(L,ft) + 1 - Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 -bc_in(s)%albgr_dif_rb(ib)) - Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(ifp,ib) * & - tr_dir_z(L,ft,iv) * (1.0_r8 -bc_in(s)%albgr_dir_rb(ib)) - tr_soild = tr_soild + ftweight(L,ft,1)*forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) - tr_soili = tr_soili + ftweight(L,ft,1)*Dif_dn(L,ft,iv) - end if - ! Absorbed radiation, shaded and sunlit portions of leaf layers - !here we get one unit of diffuse radiation... how much of - !it is absorbed? - if (ib == ivis) then ! only set the absorbed PAR for the visible light band. - do iv = 1, currentPatch%nrad(L,ft) - if (radtype==idirect) then - if ( debug ) then - write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv) - write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), & - currentPatch%fabd_sun_z(L,ft,iv) - endif - currentPatch%fabd_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - (1._r8 - currentPatch%f_sun(L,ft,iv)) - currentPatch%fabd_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - currentPatch%f_sun(L,ft,iv) + & - Abs_dir_z(ft,iv) - else - currentPatch%fabi_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - (1._r8 - currentPatch%f_sun(L,ft,iv)) - currentPatch%fabi_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - currentPatch%f_sun(L,ft,iv) - endif - if ( debug ) then - write(fates_log(),*) 'EDsurfAlb 740 ', currentPatch%fabd_sha_z(L,ft,iv), & - currentPatch%fabd_sun_z(L,ft,iv) - endif - end do - endif ! ib - - !==============================================================================! - ! Sum fluxes - !==============================================================================! - ! Solar radiation absorbed by ground - iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then - abs_rad(ib) = abs_rad(ib) + (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv)) - end if - ! Solar radiation absorbed by vegetation and sunlit/shaded leaves - do iv = 1,currentPatch%nrad(L,ft) - if (radtype == idirect)then - currentPatch%fabd(ib) = currentPatch%fabd(ib) + & - Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv) - ! bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib) - else - currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv) - ! bc_out(s)%fabi_parb(ifp,ib) = currentPatch%fabi(ib) - endif - end do - ! Albefor - if (L==1)then !top canopy layer. - if (radtype == idirect)then - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + & - Dif_up(L,ft,1) * ftweight(L,ft,1) - else - bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + & - Dif_up(L,ft,1) * ftweight(L,ft,1) - end if - end if - - ! pass normalized PAR profiles for use in diagnostic averaging for history fields - if (ib == ivis) then ! only diagnose PAR profiles for the visible band - do iv = 1, currentPatch%nrad(L,ft) - currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) - currentPatch%nrmlzd_parprof_pft_dif_z(radtype,L,ft,iv) = & - Dif_dn(L,ft,iv) + Dif_up(L,ft,iv) - ! - currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) = & - currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) + & - (forc_dir(ifp,ib) * tr_dir_z(L,ft,iv)) * & - (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) - currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) = & - currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) + & - (Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)) * & - (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) - end do - end if ! ib = visible - end if ! present - end do !ft - if (radtype == idirect)then - bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib) - else - bc_out(s)%fabi_parb(ifp,ib) = currentPatch%fabi(ib) - endif - - - !radiation absorbed from fluxes through unfilled part of lower canopy. - if (currentPatch%NCL_p > 1.and.L == currentPatch%NCL_p)then - abs_rad(ib) = abs_rad(ib) + weighted_dif_down(L-1) * & - (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-bc_in(s)%albgr_dif_rb(ib)) - abs_rad(ib) = abs_rad(ib) + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * & - (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-bc_in(s)%albgr_dir_rb(ib)) - tr_soili = tr_soili + weighted_dif_down(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) - tr_soild = tr_soild + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) - endif - - if (radtype == idirect)then - currentPatch%tr_soil_dir(ib) = tr_soild - currentPatch%tr_soil_dir_dif(ib) = tr_soili - currentPatch%sabs_dir(ib) = abs_rad(ib) - bc_out(s)%ftdd_parb(ifp,ib) = tr_soild - bc_out(s)%ftid_parb(ifp,ib) = tr_soili - else - currentPatch%tr_soil_dif(ib) = tr_soili - currentPatch%sabs_dif(ib) = abs_rad(ib) - bc_out(s)%ftii_parb(ifp,ib) = tr_soili - end if - - end do!l - - - !==============================================================================! - ! Conservation check - !==============================================================================! - ! Total radiation balance: absorbed = incoming - outgoing - - if (radtype == idirect)then - error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)) + & - currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-bc_in(s)%albgr_dif_rb(ib)))) - if ( abs(error) > 0.0001)then - write(fates_log(),*)'dir ground absorption error',ifp,s,error,currentPatch%sabs_dir(ib), & - currentPatch%tr_soil_dir(ib)* & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) - write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)) - - do ft =1,3 - iv = currentPatch%nrad(1,ft) + 1 - write(fates_log(),*) 'abs soil fluxes', Abs_dir_z(ft,iv),Abs_dif_z(ft,iv) - end do - - end if - else - if ( abs(currentPatch%sabs_dif(ib)-(currentPatch%tr_soil_dif(ib) * & - (1.0_r8-bc_in(s)%albgr_dif_rb(ib)))) > 0.0001)then - write(fates_log(),*)'dif ground absorption error',ifp,s,currentPatch%sabs_dif(ib) , & - (currentPatch%tr_soil_dif(ib)* & - (1.0_r8-bc_in(s)%albgr_dif_rb(ib))),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) - endif - endif - - if (radtype == idirect)then - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib)) - else - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabi_parb(ifp,ib) + bc_out(s)%albi_parb(ifp,ib) + currentPatch%sabs_dif(ib)) - endif - lai_reduction(:) = 0.0_r8 - do L = 1, currentPatch%NCL_p - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - do iv = 1, currentPatch%nrad(L,ft) - if (lai_change(L,ft,iv) > 0.0_r8)then - lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv)) - endif - enddo - endif - enddo - enddo + !this is just Dif down, plus refl up, plus dir intercepted and turned into dif... , + if (abs(down_rad - Dif_dn(L,ft,iv+1)) > tolerance)then + irep = 1 + end if + Dif_dn(L,ft,iv+1) = down_rad + + end do !iv + + weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & + ftweight(L,ft,1) + + endif !present + end do!ft + if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1) * & + (1.0_r8-sum(ftweight(L,1:numpft,1))) + end if + end do ! do L loop + + do L = 1, currentPatch%NCL_p ! working from the top down. + weighted_dif_up(L) = 0._r8 + do ft =1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam) + iv = currentPatch%nrad(L,ft) + 1 + if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil + Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) * currentPatch%gnd_alb_dif(ib) + & + forc_dir(radtype) * tr_dir_z(L,ft,iv) * currentPatch%gnd_alb_dir(ib) + else !In the other canopy layers, reflect off the underlying vegetation. + Dif_up(L,ft,iv) = weighted_dif_up(L+1) + end if + + ! Upward diffuse flux within and above the canopy, working upward through canopy + ! with Dif_dn from previous interation. Note: up = upward flux above current layer + do iv = currentPatch%nrad(L,ft),1,-1 + !this is radiation up, by layer transmittance, by + + !reflection of the lower layer, + up_rad = Dif_dn(L,ft,iv) * refl_dif(L,ft,iv,ib) + up_rad = up_rad + forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - exp(-k_dir(ft) * & + (currentPatch%elai_profile(L,ft,iv) + currentPatch%esai_profile(L,ft,iv)))) * & + rhol(ft,ib) + up_rad = up_rad + Dif_up(L,ft,iv+1) * tran_dif(L,ft,iv,ib) + up_rad = up_rad * ftweight(L,ft,iv)/ftweight(L,ft,1) + up_rad = up_rad + Dif_up(L,ft,iv+1) *(ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + ! THE LOWER LAYER FLUX IS HOMOGENIZED, SO WE DON"T CONSIDER THE LAI_CHANGE HERE... + + if (abs(up_rad - Dif_up(L,ft,iv)) > tolerance) then !are we close to the tolerance level? + irep = 1 + end if + Dif_up(L,ft,iv) = up_rad + + end do !iv + weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) + end if !present + end do!ft + + if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + !Add on the radiation coming up through the canopy gaps. + weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & + weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) + weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * & + weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))*currentPatch%gnd_alb_dir(ib) + end if + end do!L + end do ! do while over iter - if (radtype == idirect)then - !here we are adding a within-ED radiation scheme tolerance, and then adding the diffrence onto the albedo - !it is important that the lower boundary for this is ~1000 times smaller than the tolerance in surface albedo. - if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + error - !this terms adds the error back on to the albedo. While this is partly inexcusable, it is - ! in the medium term a solution that - ! prevents the model from crashing with small and occasional energy balances issues. - ! These are extremely difficult to debug, many have been solved already, leading - ! to the complexity of this code, but where the system generates occasional errors, we - ! will deal with them for now. - end if - if (abs(error) > 0.15_r8)then - write(fates_log(),*) 'Large Dir Radn consvn error',error ,ifp,ib - write(fates_log(),*) 'diags', bc_out(s)%albd_parb(ifp,ib), bc_out(s)%ftdd_parb(ifp,ib), & - bc_out(s)%ftid_parb(ifp,ib), bc_out(s)%fabd_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'ftweight',ftweight(1,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno - write(fates_log(),*) 'bc_in(s)%albgr_dir_rb(ib)',bc_in(s)%albgr_dir_rb(ib) - - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + error - end if - else - - if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then - bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error - end if - - if (abs(error) > 0.15_r8)then - write(fates_log(),*) '>5% Dif Radn consvn error',error ,ifp,ib - write(fates_log(),*) 'diags', bc_out(s)%albi_parb(ifp,ib), bc_out(s)%ftii_parb(ifp,ib), & - bc_out(s)%fabi_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno - write(fates_log(),*) 'bc_in(s)%albgr_dif_rb(ib)',bc_in(s)%albgr_dif_rb(ib) - write(fates_log(),*) 'rhol',rhol(1:numpft,:) - write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1) - write(fates_log(),*) 'present',currentPatch%canopy_mask(1,1:numpft) - write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1) - - bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error - end if - - if (radtype == idirect)then - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib)) - else - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabi_parb(ifp,ib) + bc_out(s)%albi_parb(ifp,ib) + currentPatch%sabs_dif(ib)) - endif - - if (abs(error) > 0.00000001_r8)then - write(fates_log(),*) 'there is still error after correction',error ,ifp,ib - end if - - end if + abs_rad(ib) = 0._r8 + tr_soili = 0._r8 + tr_soild = 0._r8 + + do L = 1, currentPatch%NCL_p !working from the top down. + abs_dir_z(:,:) = 0._r8 + abs_dif_z(:,:) = 0._r8 + do ft =1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + !==============================================================================! + ! Compute absorbed flux densities + !==============================================================================! + + ! Absorbed direct beam and diffuse do leaf layers + do iv = 1, currentPatch%nrad(L,ft) + Abs_dir_z(ft,iv) = ftweight(L,ft,iv)* forc_dir(radtype) * tr_dir_z(L,ft,iv) * & + (1.00_r8 - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & + currentPatch%esai_profile(L,ft,iv)))) * (1.00_r8 - f_not_abs(ft,ib)) + Abs_dif_z(ft,iv) = ftweight(L,ft,iv)* ((Dif_dn(L,ft,iv) + & + Dif_up(L,ft,iv+1)) * (1.00_r8 - tr_dif_z(L,ft,iv)) * & + (1.00_r8 - f_not_abs(ft,ib))) + end do + + ! Absorbed direct beam and diffuse do soil + if (L == currentPatch%NCL_p)then + iv = currentPatch%nrad(L,ft) + 1 + Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dif(ib) ) + Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(radtype) * & + tr_dir_z(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dir(ib) ) + tr_soild = tr_soild + ftweight(L,ft,1)*forc_dir(radtype) * tr_dir_z(L,ft,iv) + tr_soili = tr_soili + ftweight(L,ft,1)*Dif_dn(L,ft,iv) + end if + + ! Absorbed radiation, shaded and sunlit portions of leaf layers + !here we get one unit of diffuse radiation... how much of + !it is absorbed? + if (ib == ivis) then ! only set the absorbed PAR for the visible light band. + do iv = 1, currentPatch%nrad(L,ft) + if (radtype==idirect) then + if ( debug ) then + write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv) + write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), & + currentPatch%fabd_sun_z(L,ft,iv) + endif + currentPatch%fabd_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + (1._r8 - currentPatch%f_sun(L,ft,iv)) + currentPatch%fabd_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + currentPatch%f_sun(L,ft,iv) + & + Abs_dir_z(ft,iv) + else + currentPatch%fabi_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + (1._r8 - currentPatch%f_sun(L,ft,iv)) + currentPatch%fabi_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + currentPatch%f_sun(L,ft,iv) + endif + if ( debug ) then + write(fates_log(),*) 'EDsurfAlb 740 ', currentPatch%fabd_sha_z(L,ft,iv), & + currentPatch%fabd_sun_z(L,ft,iv) + endif + end do + endif ! ib + + + !==============================================================================! + ! Sum fluxes + !==============================================================================! + ! Solar radiation absorbed by ground + iv = currentPatch%nrad(L,ft) + 1 + if (L==currentPatch%NCL_p)then + abs_rad(ib) = abs_rad(ib) + (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv)) + end if + ! Solar radiation absorbed by vegetation and sunlit/shaded leaves + do iv = 1,currentPatch%nrad(L,ft) + if (radtype == idirect)then + currentPatch%fabd(ib) = currentPatch%fabd(ib) + & + Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv) + ! bc_out(s)%fabd_parb_out(ib) = currentPatch%fabd(ib) + else + currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv) + ! bc_out(s)%fabi_parb_out(ib) = currentPatch%fabi(ib) + endif + end do + + ! Albefor + if (L==1)then !top canopy layer. + if (radtype == idirect)then + albd_parb_out(ib) = albd_parb_out(ib) + & + Dif_up(L,ft,1) * ftweight(L,ft,1) + else + albi_parb_out(ib) = albi_parb_out(ib) + & + Dif_up(L,ft,1) * ftweight(L,ft,1) + end if + end if + + ! pass normalized PAR profiles for use in diagnostic averaging for history fields + if (ib == ivis) then ! only diagnose PAR profiles for the visible band + do iv = 1, currentPatch%nrad(L,ft) + currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = & + forc_dir(radtype) * tr_dir_z(L,ft,iv) + currentPatch%nrmlzd_parprof_pft_dif_z(radtype,L,ft,iv) = & + Dif_dn(L,ft,iv) + Dif_up(L,ft,iv) + ! + currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) = & + currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) + & + (forc_dir(radtype) * tr_dir_z(L,ft,iv)) * & + (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) + currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) = & + currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) + & + (Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)) * & + (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) + end do + end if ! ib = visible + end if ! present + end do !ft + if (radtype == idirect)then + fabd_parb_out(ib) = currentPatch%fabd(ib) + else + fabi_parb_out(ib) = currentPatch%fabi(ib) + endif + + + !radiation absorbed from fluxes through unfilled part of lower canopy. + if (currentPatch%NCL_p > 1.and.L == currentPatch%NCL_p)then + abs_rad(ib) = abs_rad(ib) + weighted_dif_down(L-1) * & + (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dif(ib) ) + abs_rad(ib) = abs_rad(ib) + forc_dir(radtype) * weighted_dir_tr(L-1) * & + (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dir(ib) ) + tr_soili = tr_soili + weighted_dif_down(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) + tr_soild = tr_soild + forc_dir(radtype) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) + endif + + if (radtype == idirect)then + currentPatch%tr_soil_dir(ib) = tr_soild + currentPatch%tr_soil_dir_dif(ib) = tr_soili + currentPatch%sabs_dir(ib) = abs_rad(ib) + ftdd_parb_out(ib) = tr_soild + ftid_parb_out(ib) = tr_soili + else + currentPatch%tr_soil_dif(ib) = tr_soili + currentPatch%sabs_dif(ib) = abs_rad(ib) + ftii_parb_out(ib) = tr_soili + end if + + end do!l + - end do !hlm_numSWb - - enddo ! rad-type - endif ! is there vegetation? + !==============================================================================! + ! Conservation check + !==============================================================================! + ! Total radiation balance: absorbed = incoming - outgoing + + if (radtype == idirect)then + error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * & + (1.0_r8-currentPatch%gnd_alb_dir(ib) ) + & + currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-currentPatch%gnd_alb_dif(ib) ))) + if ( abs(error) > 0.0001)then + write(fates_log(),*)'dir ground absorption error',error,currentPatch%sabs_dir(ib), & + currentPatch%tr_soil_dir(ib)* & + (1.0_r8-currentPatch%gnd_alb_dir(ib) ),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) + write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), & + (1.0_r8-currentPatch%gnd_alb_dir(ib) ) + + do ft =1,3 + iv = currentPatch%nrad(1,ft) + 1 + write(fates_log(),*) 'abs soil fluxes', Abs_dir_z(ft,iv),Abs_dif_z(ft,iv) + end do + + end if + else + if ( abs(currentPatch%sabs_dif(ib)-(currentPatch%tr_soil_dif(ib) * & + (1.0_r8-currentPatch%gnd_alb_dif(ib) ))) > 0.0001_r8)then + write(fates_log(),*)'dif ground absorption error',currentPatch%sabs_dif(ib) , & + (currentPatch%tr_soil_dif(ib)* & + (1.0_r8-currentPatch%gnd_alb_dif(ib) )),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) + endif + endif + + if (radtype == idirect)then + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabd_parb_out(ib) + albd_parb_out(ib) + currentPatch%sabs_dir(ib)) + else + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabi_parb_out(ib) + albi_parb_out(ib) + currentPatch%sabs_dif(ib)) + endif + lai_reduction(:) = 0.0_r8 + do L = 1, currentPatch%NCL_p + do ft =1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + do iv = 1, currentPatch%nrad(L,ft) + if (lai_change(L,ft,iv) > 0.0_r8)then + lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv)) + endif + enddo + endif + enddo + enddo + + if (radtype == idirect)then + !here we are adding a within-ED radiation scheme tolerance, and then adding the diffrence onto the albedo + !it is important that the lower boundary for this is ~1000 times smaller than the tolerance in surface albedo. + if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then + albd_parb_out(ib) = albd_parb_out(ib) + error + !this terms adds the error back on to the albedo. While this is partly inexcusable, it is + ! in the medium term a solution that + ! prevents the model from crashing with small and occasional energy balances issues. + ! These are extremely difficult to debug, many have been solved already, leading + ! to the complexity of this code, but where the system generates occasional errors, we + ! will deal with them for now. + end if + if (abs(error) > 0.15_r8)then + write(fates_log(),*) 'Large Dir Radn consvn error',error ,ib + write(fates_log(),*) 'diags', albd_parb_out(ib), ftdd_parb_out(ib), & + ftid_parb_out(ib), fabd_parb_out(ib) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'ftweight',ftweight(1,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno + write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) + + albd_parb_out(ib) = albd_parb_out(ib) + error + end if + else + + if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then + albi_parb_out(ib) = albi_parb_out(ib) + error + end if + + if (abs(error) > 0.15_r8)then + write(fates_log(),*) '>5% Dif Radn consvn error',error ,ib + write(fates_log(),*) 'diags', albi_parb_out(ib), ftii_parb_out(ib), & + fabi_parb_out(ib) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno + write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) + write(fates_log(),*) 'rhol',rhol(1:numpft,:) + write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1) + write(fates_log(),*) 'present',currentPatch%canopy_mask(1,1:numpft) + write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1) + + albi_parb_out(ib) = albi_parb_out(ib) + error + end if + + if (radtype == idirect)then + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabd_parb_out(ib) + albd_parb_out(ib) + currentPatch%sabs_dir(ib)) + else + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabi_parb_out(ib) + albi_parb_out(ib) + currentPatch%sabs_dif(ib)) + endif + + if (abs(error) > 0.00000001_r8)then + write(fates_log(),*) 'there is still error after correction',error ,ib + end if + + end if + + end do !hlm_numSWb + + enddo ! rad-type + + + end associate + return + end subroutine PatchNormanRadiation - end if ! if the vegetation and zenith filter is active - currentPatch => currentPatch%younger - end do ! Loop linked-list patches - enddo ! Loop Sites - - end associate - return - end subroutine ED_Norman_Radiation - ! ====================================================================================== subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 06264b5229..ca08a686d2 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -344,6 +344,13 @@ module EDTypesMod integer :: ncan(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft !RADIATION FLUXES + + logical :: solar_zenith_flag ! integer flag specifying daylight (based on zenith angle) + real(r8) :: solar_zenith_angle ! solar zenith angle (radians) + + real(r8) :: gnd_alb_dif(maxSWb) ! ground albedo for diffuse rad, both bands (fraction) + real(r8) :: gnd_alb_dir(maxSWb) ! ground albedo for direct rad, both bands (fraction) + real(r8) :: fabd_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of direct light absorbed by each canopy ! layer, pft, and leaf layer:- real(r8) :: fabd_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of direct light absorbed by each canopy @@ -720,6 +727,10 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%total_canopy_area = ',cpatch%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',cpatch%total_tree_area write(fates_log(),*) 'pa%zstar = ',cpatch%zstar + write(fates_log(),*) 'pa%solar_zenith_flag = ',cpatch%solar_zenith_flag + write(fates_log(),*) 'pa%solar_zenith_angle = ',cpatch%solar_zenith_angle + write(fates_log(),*) 'pa%gnd_alb_dif = ',cpatch%gnd_alb_dif(:) + write(fates_log(),*) 'pa%gnd_alb_dir = ',cpatch%gnd_alb_dir(:) write(fates_log(),*) 'pa%c_stomata = ',cpatch%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',cpatch%c_lblayer write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index e0af48ddd3..c49831e918 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -18,6 +18,8 @@ module FatesInterfaceMod use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf use EDTypesMod , only : maxpft + use EDTypesMod , only : ncwd + use EDTypesMod , only : numWaterMem use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesGlobals , only : fates_global_verbose @@ -924,8 +926,14 @@ subroutine set_fates_global_elements(use_fates) ! These values are used to define the restart file allocations and general structure ! of memory for the cohort arrays - fates_maxElementsPerPatch = max(maxCohortsPerPatch, & - numpft * nclmax * nlevleaf) + fates_maxElementsPerPatch = max(maxCohortsPerPatch, numpft, ncwd ) + + if (maxPatchesPerSite * fates_maxElementsPerPatch < numWaterMem) then + write(fates_log(), *) 'By using such a tiny number of maximum patches and maximum cohorts' + write(fates_log(), *) ' this could create problems for indexing in restart files' + write(fates_log(), *) ' The multiple of the two has to be greater than numWaterMem' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if fates_maxElementsPerSite = maxPatchesPerSite * fates_maxElementsPerPatch diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index cfc929580d..a710fd922e 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -5,12 +5,15 @@ module FatesRestartInterfaceMod use FatesConstantsMod , only : fates_avg_flag_length use FatesConstantsMod , only : fates_short_string_length use FatesConstantsMod , only : fates_long_string_length + use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : ifalse use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use FatesIODimensionsMod, only : fates_io_dimension_type use FatesIOVariableKindMod, only : fates_io_variable_kind_type use FatesRestartVariableMod, only : fates_restart_variable_type use FatesInterfaceMod, only : bc_in_type + use FatesInterfaceMod, only : bc_out_type use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index use PRTGenericMod, only : prt_global @@ -102,6 +105,9 @@ module FatesRestartInterfaceMod integer, private :: ir_lmort_collateral_co integer, private :: ir_lmort_infra_co + ! Radiation + integer, private :: ir_solar_zenith_flag_pa + integer, private :: ir_solar_zenith_angle_pa integer, private :: ir_ddbhdt_co integer, private :: ir_resp_tstep_co @@ -110,6 +116,10 @@ module FatesRestartInterfaceMod integer, private :: ir_isnew_co integer, private :: ir_cwd_ag_pacw integer, private :: ir_cwd_bg_pacw + + integer, private :: ir_gnd_alb_dif_pasb + integer, private :: ir_gnd_alb_dir_pasb + integer, private :: ir_leaf_litter_paft integer, private :: ir_root_litter_paft integer, private :: ir_leaf_litter_in_paft @@ -119,11 +129,7 @@ module FatesRestartInterfaceMod integer, private :: ir_livegrass_pa integer, private :: ir_age_pa integer, private :: ir_area_pa - integer, private :: ir_fsun_paclftls - integer, private :: ir_fabd_sun_paclftls - integer, private :: ir_fabi_sun_paclftls - integer, private :: ir_fabd_sha_paclftls - integer, private :: ir_fabi_sha_paclftls + integer, private :: ir_watermem_siwm integer, private :: ir_prt_base ! Base index for all PRT variables @@ -192,6 +198,7 @@ module FatesRestartInterfaceMod procedure, public :: set_restart_vectors procedure, public :: create_patchcohort_structure procedure, public :: get_restart_vectors + procedure, public :: update_3dpatch_radiation ! private work functions procedure, private :: init_dim_kinds_maps procedure, private :: set_dim_indices @@ -609,6 +616,16 @@ subroutine define_restart_vars(this, initialize_variables) long_name='the number of cohorts per patch', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncohort_pa ) + call this%set_restart_var(vname='fates_solar_zenith_flag_pa', vtype=cohort_int, & + long_name='switch specifying if zenith is positive', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_solar_zenith_flag_pa ) + + call this%set_restart_var(vname='fates_solar_zenith_angle_pa', vtype=cohort_r8, & + long_name='the angle of the solar zenith for each patch', units='radians', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_solar_zenith_angle_pa ) + + + ! 1D cohort Variables ! ----------------------------------------------------------------------------------- @@ -753,6 +770,16 @@ subroutine define_restart_vars(this, initialize_variables) units='kgC/m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cwd_bg_pacw ) + call this%set_restart_var(vname='fates_gnd_alb_dif', vtype=cohort_r8, & + long_name='ground albedo of diffuse radiation vis and ir', & + units='fraction', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gnd_alb_dif_pasb ) + + call this%set_restart_var(vname='fates_gnd_alb_dir', vtype=cohort_r8, & + long_name='ground albedo of direct radiation vis and ir', & + units='fraction', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gnd_alb_dir_pasb ) + call this%set_restart_var(vname='fates_leaf_litter', vtype=cohort_r8, & long_name='leaf litter, by patch x pft (non-respiring)', & units='kgC/m2', flushval = flushzero, & @@ -796,31 +823,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='are of the ED patch', units='m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_area_pa ) - ! These dimensions are pa "patch" cl "canopy layer" ft "functional type" ls "layer sublevel" - call this%set_restart_var(vname='fates_f_sun', vtype=cohort_r8, & - long_name='fraction of sunlit leaves, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fsun_paclftls ) - - call this%set_restart_var(vname='fates_fabd_sun_z', vtype=cohort_r8, & - long_name='sun fraction of direct light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabd_sun_paclftls ) - - call this%set_restart_var(vname='fates_fabi_sun_z', vtype=cohort_r8, & - long_name='sun fraction of indirect light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabi_sun_paclftls ) - - call this%set_restart_var(vname='fates_fabd_sha_z', vtype=cohort_r8, & - long_name='shade fraction of direct light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabd_sha_paclftls ) - - call this%set_restart_var(vname='fates_fabi_sha_z', vtype=cohort_r8, & - long_name='shade fraction of indirect light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabi_sha_paclftls ) ! ! site x time level vars @@ -1031,14 +1033,13 @@ end subroutine set_restart_var subroutine set_restart_vectors(this,nc,nsites,sites) - use EDTypesMod, only : nclmax - use EDTypesMod, only : nlevleaf use FatesInterfaceMod, only : fates_maxElementsPerPatch use FatesInterfaceMod, only : numpft use EDTypesMod, only : ed_site_type use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd + use EDTypesMod, only : maxSWb use EDTypesMod, only : numWaterMem ! Arguments @@ -1063,7 +1064,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation + integer :: io_idx_pa_ib ! each SW band (vis/ir) per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site ! Some counters (for checking mostly) @@ -1107,6 +1108,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & + rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & + rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & @@ -1121,7 +1124,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_gpp_acc_hold_co => this%rvars(ir_gpp_acc_hold_co)%r81d, & rio_resp_acc_hold_co => this%rvars(ir_resp_acc_hold_co)%r81d, & rio_npp_acc_hold_co => this%rvars(ir_npp_acc_hold_co)%r81d, & - rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & @@ -1137,6 +1139,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_cwd_ag_pacw => this%rvars(ir_cwd_ag_pacw)%r81d, & rio_cwd_bg_pacw => this%rvars(ir_cwd_bg_pacw)%r81d, & + rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & + rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_leaf_litter_paft => this%rvars(ir_leaf_litter_paft)%r81d, & rio_root_litter_paft => this%rvars(ir_root_litter_paft)%r81d, & rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & @@ -1146,11 +1150,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_fsun_paclftls => this%rvars(ir_fsun_paclftls)%r81d, & - rio_fabd_sun_z_paclftls => this%rvars(ir_fabd_sun_paclftls)%r81d, & - rio_fabi_sun_z_paclftls => this%rvars(ir_fabi_sun_paclftls)%r81d, & - rio_fabd_sha_z_paclftls => this%rvars(ir_fabd_sha_paclftls)%r81d, & - rio_fabi_sha_z_paclftls => this%rvars(ir_fabi_sha_paclftls)%r81d, & rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) totalCohorts = 0 @@ -1174,8 +1173,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st ! write seed_bank info(site-level, but PFT-resolved) do i = 1,numpft @@ -1300,6 +1299,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! set cohorts per patch for IO rio_ncohort_pa( io_idx_co_1st ) = cohortsperpatch + ! Set zenith angle info + if ( cpatch%solar_zenith_flag ) then + rio_solar_zenith_flag_pa(io_idx_co_1st) = itrue + else + rio_solar_zenith_flag_pa(io_idx_co_1st) = ifalse + endif + rio_solar_zenith_angle_pa( io_idx_co_1st) = cpatch%solar_zenith_angle + if ( debug ) then write(fates_log(),*) 'offsetNumCohorts III ' & ,io_idx_co,cohortsperpatch @@ -1323,25 +1330,11 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do - if ( debug ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 1 ',io_idx_pa_sunz - - if ( debug ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft,nclmax - - do k = 1,nlevleaf ! nlevleaf currently 40 - do j = 1,numpft ! dependent on parameter file - do i = 1,nclmax ! nclmax currently 2 - rio_fsun_paclftls(io_idx_pa_sunz) = cpatch%f_sun(i,j,k) - rio_fabd_sun_z_paclftls(io_idx_pa_sunz) = cpatch%fabd_sun_z(i,j,k) - rio_fabi_sun_z_paclftls(io_idx_pa_sunz) = cpatch%fabi_sun_z(i,j,k) - rio_fabd_sha_z_paclftls(io_idx_pa_sunz) = cpatch%fabd_sha_z(i,j,k) - rio_fabi_sha_z_paclftls(io_idx_pa_sunz) = cpatch%fabi_sha_z(i,j,k) - io_idx_pa_sunz = io_idx_pa_sunz + 1 - end do - end do + do i = 1,maxSWb + rio_gnd_alb_dif_pasb(io_idx_pa_ib) = cpatch%gnd_alb_dif(i) + rio_gnd_alb_dir_pasb(io_idx_pa_ib) = cpatch%gnd_alb_dir(i) + io_idx_pa_ib = io_idx_pa_ib + 1 end do - - if ( debug ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 2 ',io_idx_pa_sunz - ! Set the first cohort index to the start of the next patch, increment ! by the maximum number of cohorts per patch @@ -1350,8 +1343,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! reset counters so that they are all advanced evenly. io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_co = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st if ( debug ) then write(fates_log(),*) 'CLTV io_idx_co_1st ', io_idx_co_1st @@ -1423,8 +1416,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd - use EDTypesMod, only : nlevleaf - use EDTypesMod, only : nclmax + use EDTypesMod, only : maxSWb use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : maxpft use EDTypesMod, only : area @@ -1620,8 +1612,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd - use EDTypesMod, only : nlevleaf - use EDTypesMod, only : nclmax + use EDTypesMod, only : maxSWb use FatesInterfaceMod, only : numpft use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : numWaterMem @@ -1656,7 +1647,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation + integer :: io_idx_pa_ib ! each SW radiation band per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site ! Some counters (for checking mostly) @@ -1693,6 +1684,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & + rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & + rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & @@ -1707,17 +1700,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_gpp_acc_hold_co => this%rvars(ir_gpp_acc_hold_co)%r81d, & rio_resp_acc_hold_co => this%rvars(ir_resp_acc_hold_co)%r81d, & rio_npp_acc_hold_co => this%rvars(ir_npp_acc_hold_co)%r81d, & - rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & rio_fmort_co => this%rvars(ir_fmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & - rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & rio_lmort_infra_co => this%rvars(ir_lmort_infra_co)%r81d, & - rio_ddbhdt_co => this%rvars(ir_ddbhdt_co)%r81d, & rio_resp_tstep_co => this%rvars(ir_resp_tstep_co)%r81d, & rio_pft_co => this%rvars(ir_pft_co)%int1d, & @@ -1725,6 +1715,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_cwd_ag_pacw => this%rvars(ir_cwd_ag_pacw)%r81d, & rio_cwd_bg_pacw => this%rvars(ir_cwd_bg_pacw)%r81d, & + rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & + rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_leaf_litter_paft => this%rvars(ir_leaf_litter_paft)%r81d, & rio_root_litter_paft => this%rvars(ir_root_litter_paft)%r81d, & rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & @@ -1734,11 +1726,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_fsun_paclftls => this%rvars(ir_fsun_paclftls)%r81d, & - rio_fabd_sun_z_paclftls => this%rvars(ir_fabd_sun_paclftls)%r81d, & - rio_fabi_sun_z_paclftls => this%rvars(ir_fabi_sun_paclftls)%r81d, & - rio_fabd_sha_z_paclftls => this%rvars(ir_fabd_sha_paclftls)%r81d, & - rio_fabi_sha_z_paclftls => this%rvars(ir_fabi_sha_paclftls)%r81d, & rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) totalcohorts = 0 @@ -1751,7 +1738,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st ! read seed_bank info(site-level, but PFT-resolved) @@ -1865,17 +1852,22 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! ! deal with patch level fields here ! - cpatch%livegrass = rio_livegrass_pa(io_idx_co_1st) - cpatch%age = rio_age_pa(io_idx_co_1st) - cpatch%area = rio_area_pa(io_idx_co_1st) - cpatch%age_class = get_age_class_index(cpatch%age) - + cpatch%livegrass = rio_livegrass_pa(io_idx_co_1st) + cpatch%age = rio_age_pa(io_idx_co_1st) + cpatch%area = rio_area_pa(io_idx_co_1st) + cpatch%age_class = get_age_class_index(cpatch%age) + + ! Set zenith angle info + cpatch%solar_zenith_flag = ( rio_solar_zenith_flag_pa(io_idx_co_1st) .eq. itrue ) + cpatch%solar_zenith_angle = rio_solar_zenith_angle_pa(io_idx_co_1st) + ! set cohorts per patch for IO if ( debug ) then write(fates_log(),*) 'CVTL III ' & ,io_idx_co,cohortsperpatch endif + ! ! deal with patch level fields of arrays here ! @@ -1896,23 +1888,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 enddo - if ( debug ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 1 ',io_idx_pa_sunz - - do k = 1,nlevleaf ! nlevleaf currently 40 - do j = 1,numpft - do i = 1,nclmax ! nclmax currently 2 - cpatch%f_sun(i,j,k) = rio_fsun_paclftls(io_idx_pa_sunz) - cpatch%fabd_sun_z(i,j,k) = rio_fabd_sun_z_paclftls(io_idx_pa_sunz) - cpatch%fabi_sun_z(i,j,k) = rio_fabi_sun_z_paclftls(io_idx_pa_sunz) - cpatch%fabd_sha_z(i,j,k) = rio_fabd_sha_z_paclftls(io_idx_pa_sunz) - cpatch%fabi_sha_z(i,j,k) = rio_fabi_sha_z_paclftls(io_idx_pa_sunz) - io_idx_pa_sunz = io_idx_pa_sunz + 1 - end do - end do + do i = 1,maxSWb + cpatch%gnd_alb_dif(i) = rio_gnd_alb_dif_pasb(io_idx_pa_ib) + cpatch%gnd_alb_dir(i) = rio_gnd_alb_dir_pasb(io_idx_pa_ib) + io_idx_pa_ib = io_idx_pa_ib + 1 end do - - if ( debug ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 2 ',io_idx_pa_sunz - + ! Now increment the position of the first cohort to that of the next ! patch @@ -1921,8 +1902,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! and max the number of allowed cohorts per patch io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_co = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st if ( debug ) then write(fates_log(),*) 'CVTL io_idx_co_1st ', io_idx_co_1st @@ -1979,4 +1960,111 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end associate end subroutine get_restart_vectors + + ! ==================================================================================== + + subroutine update_3dpatch_radiation(this, nc, nsites, sites, bc_out) + + ! ------------------------------------------------------------------------- + ! This subroutine populates output boundary conditions related to radiation + ! called upon restart reads. + ! ------------------------------------------------------------------------- + + use EDTypesMod, only : ed_site_type + use EDTypesMod, only : ed_patch_type + use EDSurfaceRadiationMod, only : PatchNormanRadiation + use FatesInterfaceMod, only : hlm_numSWb + + ! !ARGUMENTS: + class(fates_restart_interface_type) , intent(inout) :: this + integer , intent(in) :: nc + integer , intent(in) :: nsites + type(ed_site_type) , intent(inout), target :: sites(nsites) + type(bc_out_type) , intent(inout) :: bc_out(nsites) + + ! locals + ! ---------------------------------------------------------------------------------- + type(ed_patch_type),pointer :: currentPatch ! current patch + integer :: s ! site counter + integer :: ib ! radiation band counter + integer :: ifp ! patch counter + + do s = 1, nsites + + ifp = 0 + currentpatch => sites(s)%oldest_patch + do while (associated(currentpatch)) + ifp = ifp+1 + + currentPatch%f_sun (:,:,:) = 0._r8 + currentPatch%fabd_sun_z (:,:,:) = 0._r8 + currentPatch%fabd_sha_z (:,:,:) = 0._r8 + currentPatch%fabi_sun_z (:,:,:) = 0._r8 + currentPatch%fabi_sha_z (:,:,:) = 0._r8 + currentPatch%fabd (:) = 0._r8 + currentPatch%fabi (:) = 0._r8 + + ! zero diagnostic radiation profiles + currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_dir_z(:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_dif_z(:,:,:) = 0._r8 + + ! ----------------------------------------------------------- + ! When calling norman radiation from the short-timestep + ! we are passing in boundary conditions to set the following + ! variables: + ! currentPatch%solar_zenith_flag (is there daylight?) + ! currentPatch%solar_zenith_angle (what is the value?) + ! ----------------------------------------------------------- + + if(currentPatch%solar_zenith_flag)then + + bc_out(s)%albd_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%albi_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%fabi_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%fabd_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%ftdd_parb(ifp,:) = 1._r8 ! output HLM + bc_out(s)%ftid_parb(ifp,:) = 1._r8 ! output HLM + bc_out(s)%ftii_parb(ifp,:) = 1._r8 ! output HLM + + if (maxval(currentPatch%nrad(1,:))==0)then + !there are no leaf layers in this patch. it is effectively bare ground. + ! no radiation is absorbed + bc_out(s)%fabd_parb(ifp,:) = 0.0_r8 + bc_out(s)%fabi_parb(ifp,:) = 0.0_r8 + do ib = 1,hlm_numSWb + + ! REQUIRES A FIX HERE albd vs albi + + bc_out(s)%albd_parb(ifp,ib) = currentPatch%gnd_alb_dir(ib) + bc_out(s)%albd_parb(ifp,ib) = currentPatch%gnd_alb_dif(ib) + bc_out(s)%ftdd_parb(ifp,ib)= 1.0_r8 + bc_out(s)%ftid_parb(ifp,ib)= 1.0_r8 + bc_out(s)%ftii_parb(ifp,ib)= 1.0_r8 + enddo + else + + call PatchNormanRadiation (currentPatch, & + bc_out(s)%albd_parb(ifp,:), & + bc_out(s)%albi_parb(ifp,:), & + bc_out(s)%fabd_parb(ifp,:), & + bc_out(s)%fabi_parb(ifp,:), & + bc_out(s)%ftdd_parb(ifp,:), & + bc_out(s)%ftid_parb(ifp,:), & + bc_out(s)%ftii_parb(ifp,:)) + + endif ! is there vegetation? + + end if ! if the vegetation and zenith filter is active + + + currentPatch => currentPatch%younger + end do ! Loop linked-list patches + enddo ! Loop Sites + + return + end subroutine update_3dpatch_radiation + + end module FatesRestartInterfaceMod