diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index e859a7339f..46ee85d045 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -49,7 +49,7 @@ module EDCanopyStructureMod use PRTGenericMod, only : carbon12_element use FatesTwoStreamUtilsMod, only : FatesConstructRadElements use FatesRadiationMemMod , only : twostr_solver - use FatesRadiationMemMod , only : num_rad_stream_types + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -1909,127 +1909,120 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) c = fcolumn(s) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! ignore the bare-ground-PFT patch entirely for these BC outs + ifp = ifp+1 - ifp = ifp+1 + if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then + if(debug)then + write(fates_log(),*) 'ED: canopy area bigger than area', & + currentPatch%total_canopy_area ,currentPatch%area + end if + currentPatch%total_canopy_area = currentPatch%area + endif - if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then - if(debug)then - write(fates_log(),*) 'ED: canopy area bigger than area', & - currentPatch%total_canopy_area ,currentPatch%area - end if - currentPatch%total_canopy_area = currentPatch%area - endif + if (associated(currentPatch%tallest)) then + bc_out(s)%htop_pa(ifp) = currentPatch%tallest%height + else + ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height? + bc_out(s)%htop_pa(ifp) = 0.1_r8 + endif - if (associated(currentPatch%tallest)) then - bc_out(s)%htop_pa(ifp) = currentPatch%tallest%height - else - ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height? - bc_out(s)%htop_pa(ifp) = 0.1_r8 - endif + bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8)) - bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8)) + ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic + ! Roughness length and displacement height used by the HLM + ! use total LAI + SAI to weight the leaft characteristic dimension + ! Avoid this if running in satellite phenology mode + ! ---------------------------------------------------------------------------- - ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic - ! Roughness length and displacement height used by the HLM - ! use total LAI + SAI to weight the leaft characteristic dimension - ! Avoid this if running in satellite phenology mode - ! ---------------------------------------------------------------------------- + if (currentPatch%total_canopy_area > nearzero) then + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if (currentCohort%canopy_layer .eq. 1) then + weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area) + bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + & + EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%height * weight + bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + & + EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%height * weight + endif + currentCohort => currentCohort%taller + end do - if (currentPatch%total_canopy_area > nearzero) then - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if (currentCohort%canopy_layer .eq. 1) then - weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area) - bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + & - EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%height * weight - bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + & - EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%height * weight - endif - currentCohort => currentCohort%taller - end do + ! for lai, scale to total LAI + SAI in patch. first add up all the LAI and SAI in the patch + total_patch_leaf_stem_area = 0._r8 + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + total_patch_leaf_stem_area = total_patch_leaf_stem_area + & + (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area + currentCohort => currentCohort%taller + end do - ! for lai, scale to total LAI + SAI in patch. first add up all the LAI and SAI in the patch - total_patch_leaf_stem_area = 0._r8 + ! make sure there is some leaf and stem area + if (total_patch_leaf_stem_area > nearzero) then currentCohort => currentPatch%shortest do while(associated(currentCohort)) - total_patch_leaf_stem_area = total_patch_leaf_stem_area + & - (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area + ! weight dleaf by the relative totals of leaf and stem area + weight = (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area / total_patch_leaf_stem_area + bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + & + EDPftvarcon_inst%dleaf(currentCohort%pft) * weight currentCohort => currentCohort%taller end do - - ! make sure there is some leaf and stem area - if (total_patch_leaf_stem_area > nearzero) then - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - ! weight dleaf by the relative totals of leaf and stem area - weight = (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area / total_patch_leaf_stem_area - bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + & - EDPftvarcon_inst%dleaf(currentCohort%pft) * weight - currentCohort => currentCohort%taller - end do - else - ! dummy case - bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1) - endif else - ! if no canopy, then use dummy values (first PFT) of aerodynamic properties - bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp) - bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp) + ! dummy case bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1) endif - ! ----------------------------------------------------------------------------- - - ! We are assuming here that grass is all located underneath tree canopies. - ! The alternative is to assume it is all spatial distinct from tree canopies. - ! In which case, the bare area would have to be reduced by the grass area... - ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants - ! currentPatch%area/AREA is the fraction of the soil covered by this patch. - - if(currentPatch%area.gt.0.0_r8)then - bc_out(s)%canopy_fraction_pa(ifp) = & - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA) - else - bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8 - endif - - bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * & - (currentPatch%area/AREA) + else + ! if no canopy, then use dummy values (first PFT) of aerodynamic properties + bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp) + bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp) + bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1) + endif + ! ----------------------------------------------------------------------------- - total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area + ! We are assuming here that grass is all located underneath tree canopies. + ! The alternative is to assume it is all spatial distinct from tree canopies. + ! In which case, the bare area would have to be reduced by the grass area... + ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants + ! currentPatch%area/AREA is the fraction of the soil covered by this patch. - total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp) + if(currentPatch%area.gt.0.0_r8)then + bc_out(s)%canopy_fraction_pa(ifp) = & + min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA) + else + bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8 + endif - bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label + bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * & + (currentPatch%area/AREA) - ! Calculate area indices for output boundary to HLM - ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles - ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called) + total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area - bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai') - bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai') - bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai') - bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai') + total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp) - ! Fraction of vegetation free of snow. This is used to flag those - ! patches which shall under-go photosynthesis - ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let - ! FATES internal variables decide if photosynthesis is possible - ! we are essentially calculating it inside FATES to tell the - ! host to tell itself when to do things (circuitous). Just have - ! to determine where else it is used + bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label - if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then - bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8 - else - bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8 - end if + ! Calculate area indices for output boundary to HLM + ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles + ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called) - else ! nocomp or SP, and currentPatch%nocomp_pft_label .eq. 0 + bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai') + bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai') + bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai') + bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai') - total_patch_area = total_patch_area + currentPatch%area/AREA + ! Fraction of vegetation free of snow. This is used to flag those + ! patches which shall under-go photosynthesis + ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let + ! FATES internal variables decide if photosynthesis is possible + ! we are essentially calculating it inside FATES to tell the + ! host to tell itself when to do things (circuitous). Just have + ! to determine where else it is used + if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then + bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8 + else + bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8 end if + currentPatch => currentPatch%younger end do @@ -2051,11 +2044,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) currentPatch => sites(s)%oldest_patch ifp = 0 do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only - ifp = ifp+1 - bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area - endif ! veg patch - + ifp = ifp+1 + bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area currentPatch => currentPatch%younger end do diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 9b02c9e14a..55bbbbe557 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1899,22 +1899,6 @@ subroutine set_patchno( currentSite ) currentPatch => currentPatch%younger enddo - if(hlm_use_fixed_biogeog.eq.itrue .and. hlm_use_nocomp.eq.itrue)then - patchno = 1 - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label.eq.nocomp_bareground)then - ! for bareground patch, we make the patch number 0 - ! we also do not count this in the veg. patch numbering scheme. - currentPatch%patchno = 0 - else - currentPatch%patchno = patchno - patchno = patchno + 1 - endif - currentPatch => currentPatch%younger - enddo - endif - end subroutine set_patchno ! ============================================================================ diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 8d689a759a..9f0db5ff39 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -34,7 +34,6 @@ module EDPhysiologyMod use FatesConstantsMod, only : mpa_per_mm_suction use FatesConstantsMod, only : g_per_kg use FatesConstantsMod, only : ndays_per_year - use FatesConstantsMod, only : nocomp_bareground use FatesConstantsMod, only : nocomp_bareground_land use FatesConstantsMod, only : is_crop use FatesConstantsMod, only : area_error_2 @@ -3175,44 +3174,40 @@ subroutine fragmentation_scaler( currentPatch, bc_in) catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8 * ( t1 - 15.4_r8 )) catanf_30 = catanf(30._r8) - if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then + ! Use the hlm temp and moisture decomp fractions by default + if ( use_hlm_soil_scalar ) then - ! Use the hlm temp and moisture decomp fractions by default - if ( use_hlm_soil_scalar ) then + ! Calculate the fragmentation_scaler + currentPatch%fragmentation_scaler = min(1.0_r8,max(0.0_r8,bc_in%t_scalar_sisl * bc_in%w_scalar_sisl)) - ! Calculate the fragmentation_scaler - currentPatch%fragmentation_scaler = min(1.0_r8,max(0.0_r8,bc_in%t_scalar_sisl * bc_in%w_scalar_sisl)) + else + if ( .not. use_century_tfunc ) then + !calculate rate constant scalar for soil temperature,assuming that the base rate constants + !are assigned for non-moisture limiting conditions at 25C. + if (currentPatch%tveg24%GetMean() >= tfrz) then + t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8) + ! Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8) + else + t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8)) + ! Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8) + endif else + ! original century uses an arctangent function to calculate the + ! temperature dependence of decomposition + t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8) + endif - if ( .not. use_century_tfunc ) then - !calculate rate constant scalar for soil temperature,assuming that the base rate constants - !are assigned for non-moisture limiting conditions at 25C. - if (currentPatch%tveg24%GetMean() >= tfrz) then - t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8) - ! Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8) - else - t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8)) - ! Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8) - endif - else - ! original century uses an arctangent function to calculate the - ! temperature dependence of decomposition - t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8) - endif - - !Moisture Limitations - !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed - !soil moisture values, which is not realistic. - !litter decomp is proportional to water limitation on average... - w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8) + !Moisture Limitations + !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed + !soil moisture values, which is not realistic. + !litter decomp is proportional to water limitation on average... + w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8) ! Calculate the fragmentation_scaler - currentPatch%fragmentation_scaler(:) = min(1.0_r8,max(0.0_r8,t_scalar * w_scalar)) - - endif ! scalar + currentPatch%fragmentation_scaler(:) = min(1.0_r8,max(0.0_r8,t_scalar * w_scalar)) - endif ! not bare ground + endif ! scalar end subroutine fragmentation_scaler diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index 2dd3f816a4..3120f13e35 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -69,7 +69,6 @@ module FatesSoilBGCFluxMod use FatesConstantsMod, only : sec_per_day use FatesConstantsMod, only : years_per_day use FatesConstantsMod, only : itrue - use FatesConstantsMod, only : nocomp_bareground use FatesLitterMod, only : litter_type use FatesLitterMod , only : ncwd use FatesLitterMod , only : ndcmpy @@ -287,108 +286,107 @@ subroutine PrepCH4BCs(csite,bc_in,bc_out) fp = 0 cpatch => csite%oldest_patch do while (associated(cpatch)) - if_notbare: if(cpatch%nocomp_pft_label .ne. nocomp_bareground)then - ! Patch ordering when passing boundary conditions - ! always goes from oldest to youngest, following - ! the convention of EDPatchDynamics::set_patchno() - - fp = fp + 1 - - agnpp = 0._r8 - bgnpp = 0._r8 - woody_area = 0._r8 - plant_area = 0._r8 - - ccohort => cpatch%tallest - do while (associated(ccohort)) - - ! For consistency, only apply calculations to non-new - ! cohorts. New cohorts will not have respiration rates - ! at this point in the call sequence. - - if(.not.ccohort%isnew) then - - pft = ccohort%pft - - call set_root_fraction(csite%rootfrac_scr, pft, csite%zi_soil, & - bc_in%max_rooting_depth_index_col ) - - fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element) - - ! [kgC/day] - sapw_net_alloc = ccohort%prt%GetNetAlloc(sapw_organ, carbon12_element) * days_per_sec - store_net_alloc = ccohort%prt%GetNetAlloc(store_organ, carbon12_element) * days_per_sec - leaf_net_alloc = ccohort%prt%GetNetAlloc(leaf_organ, carbon12_element) * days_per_sec - fnrt_net_alloc = ccohort%prt%GetNetAlloc(fnrt_organ, carbon12_element) * days_per_sec - struct_net_alloc = ccohort%prt%GetNetAlloc(struct_organ, carbon12_element) * days_per_sec - repro_net_alloc = ccohort%prt%GetNetAlloc(repro_organ, carbon12_element) * days_per_sec - - ! [kgC/plant/day] -> [gC/m2/s] - agnpp = agnpp + ccohort%n/cpatch%area * (leaf_net_alloc + repro_net_alloc + & - prt_params%allom_agb_frac(pft)*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg - - ! [kgC/plant/day] -> [gC/m2/s] - bgnpp = bgnpp + ccohort%n/cpatch%area * (fnrt_net_alloc + & - (1._r8-prt_params%allom_agb_frac(pft))*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg - - if(hlm_use_ch4==itrue)then - - ! Fine root fraction over depth - bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = & - bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) + & - csite%rootfrac_scr(1:bc_in%nlevsoil) - - ! Fine root carbon, convert [kg/plant] -> [g/m2] - bc_out%frootc_pa(fp) = & - bc_out%frootc_pa(fp) + & - fnrt_c*ccohort%n/cpatch%area * g_per_kg - - ! (gC/m2/s) root respiration (fine root MR + total root GR) - ! RGK: We do not save root respiration and average over the day. Until we do - ! this is a best (bad) guess at fine root MR + total root GR - ! (kgC/indiv/yr) -> gC/m2/s - bc_out%root_resp(1:bc_in%nlevsoil) = bc_out%root_resp(1:bc_in%nlevsoil) + & - (ccohort%resp_m_acc_hold + ccohort%resp_g_acc_hold)*years_per_day*g_per_kg*days_per_sec* & - ccohort%n*area_inv*(1._r8-prt_params%allom_agb_frac(pft)) * csite%rootfrac_scr(1:bc_in%nlevsoil) - - - end if - - if( prt_params%woody(pft)==itrue ) then - woody_area = woody_area + ccohort%c_area - end if - plant_area = plant_area + ccohort%c_area - - - end if - - ccohort => ccohort%shorter - end do - - if(hlm_use_ch4==itrue)then - if( sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) > nearzero) then + + ! Patch ordering when passing boundary conditions + ! always goes from oldest to youngest, following + ! the convention of EDPatchDynamics::set_patchno() + + fp = fp + 1 + + agnpp = 0._r8 + bgnpp = 0._r8 + woody_area = 0._r8 + plant_area = 0._r8 + + ccohort => cpatch%tallest + do while (associated(ccohort)) + + ! For consistency, only apply calculations to non-new + ! cohorts. New cohorts will not have respiration rates + ! at this point in the call sequence. + + if(.not.ccohort%isnew) then + + pft = ccohort%pft + + call set_root_fraction(csite%rootfrac_scr, pft, csite%zi_soil, & + bc_in%max_rooting_depth_index_col ) + + fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element) + + ! [kgC/day] + sapw_net_alloc = ccohort%prt%GetNetAlloc(sapw_organ, carbon12_element) * days_per_sec + store_net_alloc = ccohort%prt%GetNetAlloc(store_organ, carbon12_element) * days_per_sec + leaf_net_alloc = ccohort%prt%GetNetAlloc(leaf_organ, carbon12_element) * days_per_sec + fnrt_net_alloc = ccohort%prt%GetNetAlloc(fnrt_organ, carbon12_element) * days_per_sec + struct_net_alloc = ccohort%prt%GetNetAlloc(struct_organ, carbon12_element) * days_per_sec + repro_net_alloc = ccohort%prt%GetNetAlloc(repro_organ, carbon12_element) * days_per_sec + + ! [kgC/plant/day] -> [gC/m2/s] + agnpp = agnpp + ccohort%n/cpatch%area * (leaf_net_alloc + repro_net_alloc + & + prt_params%allom_agb_frac(pft)*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg + + ! [kgC/plant/day] -> [gC/m2/s] + bgnpp = bgnpp + ccohort%n/cpatch%area * (fnrt_net_alloc + & + (1._r8-prt_params%allom_agb_frac(pft))*(sapw_net_alloc+store_net_alloc+struct_net_alloc)) * g_per_kg + + if(hlm_use_ch4==itrue)then + + ! Fine root fraction over depth bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = & - bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) / & - sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) + bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) + & + csite%rootfrac_scr(1:bc_in%nlevsoil) + + ! Fine root carbon, convert [kg/plant] -> [g/m2] + bc_out%frootc_pa(fp) = & + bc_out%frootc_pa(fp) + & + fnrt_c*ccohort%n/cpatch%area * g_per_kg + + ! (gC/m2/s) root respiration (fine root MR + total root GR) + ! RGK: We do not save root respiration and average over the day. Until we do + ! this is a best (bad) guess at fine root MR + total root GR + ! (kgC/indiv/yr) -> gC/m2/s + bc_out%root_resp(1:bc_in%nlevsoil) = bc_out%root_resp(1:bc_in%nlevsoil) + & + (ccohort%resp_m_acc_hold + ccohort%resp_g_acc_hold)*years_per_day*g_per_kg*days_per_sec* & + ccohort%n*area_inv*(1._r8-prt_params%allom_agb_frac(pft)) * csite%rootfrac_scr(1:bc_in%nlevsoil) + + end if - - ! RGK: These averages should switch to the new patch averaging methods - ! when available. Right now we are not doing any time averaging - ! because it would be mixing the memory of patches, which - ! would be arguably worse than just using the instantaneous value - - ! gC/m2/s - bc_out%annavg_agnpp_pa(fp) = agnpp - bc_out%annavg_bgnpp_pa(fp) = bgnpp - ! gc/m2/yr - bc_out%annsum_npp_pa(fp) = (bgnpp+agnpp)*days_per_year*sec_per_day - - if(plant_area>nearzero) then - bc_out%woody_frac_aere_pa(fp) = woody_area/plant_area + + if( prt_params%woody(pft)==itrue ) then + woody_area = woody_area + ccohort%c_area end if - + plant_area = plant_area + ccohort%c_area + + + end if + + ccohort => ccohort%shorter + end do + + if(hlm_use_ch4==itrue)then + if( sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) > nearzero) then + bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) = & + bc_out%rootfr_pa(fp,1:bc_in%nlevsoil) / & + sum(bc_out%rootfr_pa(fp,1:bc_in%nlevsoil)) end if - end if if_notbare + + ! RGK: These averages should switch to the new patch averaging methods + ! when available. Right now we are not doing any time averaging + ! because it would be mixing the memory of patches, which + ! would be arguably worse than just using the instantaneous value + + ! gC/m2/s + bc_out%annavg_agnpp_pa(fp) = agnpp + bc_out%annavg_bgnpp_pa(fp) = bgnpp + ! gc/m2/yr + bc_out%annsum_npp_pa(fp) = (bgnpp+agnpp)*days_per_year*sec_per_day + + if(plant_area>nearzero) then + bc_out%woody_frac_aere_pa(fp) = woody_area/plant_area + end if + end if + cpatch => cpatch%younger end do diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 5bab6b5191..32c5b38bcd 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -13,7 +13,6 @@ module EDAccumulateFluxesMod use FatesGlobals, only : fates_log use shr_log_mod , only : errMsg => shr_log_errMsg use FatesConstantsMod , only : r8 => fates_r8 - use FatesConstantsMod , only : nocomp_bareground implicit none @@ -70,41 +69,39 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) cpatch => sites(s)%oldest_patch do while (associated(cpatch)) - if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then - ifp = ifp+1 - - if( bc_in(s)%filter_photo_pa(ifp) == 3 ) then - ccohort => cpatch%shortest - do while(associated(ccohort)) - - ! Accumulate fluxes from hourly to daily values. - ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day - - ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep - ccohort%resp_m_acc = ccohort%resp_m_acc + ccohort%resp_m_tstep - - ccohort%sym_nfix_daily = ccohort%sym_nfix_daily + ccohort%sym_nfix_tstep - - ! weighted mean of D13C by gpp - if((ccohort%gpp_acc + ccohort%gpp_tstep) .eq. 0.0_r8) then - ccohort%c13disc_acc = 0.0_r8 - else - ccohort%c13disc_acc = ((ccohort%c13disc_acc * ccohort%gpp_acc) + & - (ccohort%c13disc_clm * ccohort%gpp_tstep)) / & - (ccohort%gpp_acc + ccohort%gpp_tstep) - endif - - do iv=1,ccohort%nv - if(ccohort%year_net_uptake(iv) == 999._r8)then ! note that there were leaves in this layer this year. - ccohort%year_net_uptake(iv) = 0._r8 - end if - ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) - enddo - - ccohort => ccohort%taller - enddo ! while(associated(ccohort)) - end if - end if ! not bare ground + ifp = ifp+1 + + if( bc_in(s)%filter_photo_pa(ifp) == 3 ) then + ccohort => cpatch%shortest + do while(associated(ccohort)) + + ! Accumulate fluxes from hourly to daily values. + ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day + + ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep + ccohort%resp_m_acc = ccohort%resp_m_acc + ccohort%resp_m_tstep + + ccohort%sym_nfix_daily = ccohort%sym_nfix_daily + ccohort%sym_nfix_tstep + + ! weighted mean of D13C by gpp + if((ccohort%gpp_acc + ccohort%gpp_tstep) .eq. 0.0_r8) then + ccohort%c13disc_acc = 0.0_r8 + else + ccohort%c13disc_acc = ((ccohort%c13disc_acc * ccohort%gpp_acc) + & + (ccohort%c13disc_clm * ccohort%gpp_tstep)) / & + (ccohort%gpp_acc + ccohort%gpp_tstep) + endif + + do iv=1,ccohort%nv + if(ccohort%year_net_uptake(iv) == 999._r8)then ! note that there were leaves in this layer this year. + ccohort%year_net_uptake(iv) = 0._r8 + end if + ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) + enddo + + ccohort => ccohort%taller + enddo ! while(associated(ccohort)) + end if cpatch => cpatch%younger end do ! while(associated(cpatch)) end do diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 3e2401a033..17f4356b1d 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -8,7 +8,6 @@ module EDBtranMod use EDPftvarcon , only : EDPftvarcon_inst use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod , only : itrue,ifalse,nearzero - use FatesConstantsMod , only : nocomp_bareground use EDTypesMod , only : ed_site_type use FatesPatchMod, only : fates_patch_type use EDParamsMod, only : maxpft @@ -131,128 +130,128 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) smpso => EDPftvarcon_inst%smpso & ! INTERF-TODO: THESE SHOULD BE FATES PARAMETERS ) - do s = 1,nsites + do s = 1,nsites - allocate(root_resis(numpft,bc_in(s)%nlevsoil)) + allocate(root_resis(numpft,bc_in(s)%nlevsoil)) - bc_out(s)%rootr_pasl(:,:) = 0._r8 + bc_out(s)%rootr_pasl(:,:) = 0._r8 - ifp = 0 - cpatch => sites(s)%oldest_patch - do while (associated(cpatch)) - if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then ! only for veg patches - ifp=ifp+1 + ifp = 0 + cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) - ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) - - do ft = 1,numpft - - call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, & - bc_in(s)%max_rooting_depth_index_col ) - - cpatch%btran_ft(ft) = 0.0_r8 - do j = 1,bc_in(s)%nlevsoil - - ! Calculations are only relevant where liquid water exists - ! see clm_fates%wrap_btran for calculation with CLM/ALM - - if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then - - smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) - - rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & - (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) - - root_resis(ft,j) = sites(s)%rootfrac_scr(j)*rresis - - ! root water uptake is not linearly proportional to root density, - ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) - - cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + root_resis(ft,j) - - else - root_resis(ft,j) = 0._r8 - end if - - end do !j - - ! Normalize root resistances to get layer contribution to ET - do j = 1,bc_in(s)%nlevsoil - if (cpatch%btran_ft(ft) > nearzero) then - root_resis(ft,j) = root_resis(ft,j)/cpatch%btran_ft(ft) - else - root_resis(ft,j) = 0._r8 - end if - end do - - end do !PFT - - ! PFT-averaged point level root fraction for extraction purposese. - ! The cohort's conductance g_sb_laweighted, contains a weighting factor - ! based on the cohort's leaf area. units: [m/s] * [m2] - - pftgs(1:maxpft) = 0._r8 - ccohort => cpatch%tallest - do while(associated(ccohort)) - pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%g_sb_laweight - ccohort => ccohort%shorter - enddo - - ! Process the boundary output, this is necessary for calculating the soil-moisture - ! sink term across the different layers in driver/host. Photosynthesis will - ! pass the host a total transpiration for the patch. This needs rootr to be - ! distributed over the soil layers. - sum_pftgs = sum(pftgs(1:numpft)) - - do j = 1, bc_in(s)%nlevsoil - bc_out(s)%rootr_pasl(ifp,j) = 0._r8 - do ft = 1,numpft - if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & - root_resis(ft,j) * pftgs(ft)/sum_pftgs - else - bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & - root_resis(ft,j) * 1._r8/real(numpft,r8) - end if - enddo - enddo - - ! Calculate the BTRAN that is passed back to the HLM - ! used only for diagnostics. If plant hydraulics is turned off - ! we are using the patchxpft level btran calculation - - if(hlm_use_planthydro.eq.ifalse) then - !weight patch level output BTRAN for the - bc_out(s)%btran_pa(ifp) = 0.0_r8 - do ft = 1,numpft - if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum_pftgs - else - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft - end if - enddo - end if - - temprootr = sum(bc_out(s)%rootr_pasl(ifp,1:bc_in(s)%nlevsoil)) - - if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then - - if(debug) write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs - - do j = 1,bc_in(s)%nlevsoil - bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j)/temprootr - enddo - - end if - endif ! not bare ground - cpatch => cpatch%younger - end do - - deallocate(root_resis) + ifp=ifp+1 - end do + ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) + + do ft = 1,numpft + + call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, & + bc_in(s)%max_rooting_depth_index_col ) + + cpatch%btran_ft(ft) = 0.0_r8 + do j = 1,bc_in(s)%nlevsoil + + ! Calculations are only relevant where liquid water exists + ! see clm_fates%wrap_btran for calculation with CLM/ALM + + if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then + + smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) + + rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & + (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) + + root_resis(ft,j) = sites(s)%rootfrac_scr(j)*rresis + + ! root water uptake is not linearly proportional to root density, + ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) + + cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + root_resis(ft,j) + + else + root_resis(ft,j) = 0._r8 + end if + + end do !j + + ! Normalize root resistances to get layer contribution to ET + do j = 1,bc_in(s)%nlevsoil + if (cpatch%btran_ft(ft) > nearzero) then + root_resis(ft,j) = root_resis(ft,j)/cpatch%btran_ft(ft) + else + root_resis(ft,j) = 0._r8 + end if + end do + + end do !PFT + + ! PFT-averaged point level root fraction for extraction purposese. + ! The cohort's conductance g_sb_laweighted, contains a weighting factor + ! based on the cohort's leaf area. units: [m/s] * [m2] + + pftgs(1:maxpft) = 0._r8 + ccohort => cpatch%tallest + do while(associated(ccohort)) + pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%g_sb_laweight + ccohort => ccohort%shorter + enddo + + ! Process the boundary output, this is necessary for calculating the soil-moisture + ! sink term across the different layers in driver/host. Photosynthesis will + ! pass the host a total transpiration for the patch. This needs rootr to be + ! distributed over the soil layers. + sum_pftgs = sum(pftgs(1:numpft)) + + do j = 1, bc_in(s)%nlevsoil + bc_out(s)%rootr_pasl(ifp,j) = 0._r8 + do ft = 1,numpft + if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail + !bit-retart test as a result? FIX(RF,032414) + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & + root_resis(ft,j) * pftgs(ft)/sum_pftgs + else + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & + root_resis(ft,j) * 1._r8/real(numpft,r8) + end if + enddo + enddo + + ! Calculate the BTRAN that is passed back to the HLM + ! used only for diagnostics. If plant hydraulics is turned off + ! we are using the patchxpft level btran calculation + + if(hlm_use_planthydro.eq.ifalse) then + !weight patch level output BTRAN for the + bc_out(s)%btran_pa(ifp) = 0.0_r8 + do ft = 1,numpft + if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail + !bit-retart test as a result? FIX(RF,032414) + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum_pftgs + else + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft + end if + enddo + end if + + temprootr = sum(bc_out(s)%rootr_pasl(ifp,1:bc_in(s)%nlevsoil)) + + if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then + + if(debug) write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs + + do j = 1,bc_in(s)%nlevsoil + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j)/temprootr + enddo + + end if + + cpatch => cpatch%younger + end do + + deallocate(root_resis) + + end do if(hlm_use_planthydro.eq.itrue) then call BTranForHLMDiagnosticsFromCohortHydr(nsites,sites,bc_out) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 454bf5d28d..f34d2250b5 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -42,7 +42,6 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : cm3_per_m3 use FatesConstantsMod, only : kg_per_g use FatesConstantsMod, only : fates_unset_r8 - use FatesConstantsMod, only : nocomp_bareground use EDParamsMod , only : hydr_kmax_rsurf1 use EDParamsMod , only : hydr_kmax_rsurf2 @@ -2499,179 +2498,178 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ifp = 0 cpatch => sites(s)%oldest_patch do while (associated(cpatch)) + + ifp = ifp + 1 - if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then - ifp = ifp + 1 - - ! ---------------------------------------------------------------------------- - ! Objective: Partition the transpiration flux - ! specfied by the land model to the cohorts. The weighting - ! factor we use to downscale is the cohort combo term: g_sb_laweight - ! This term is the stomatal conductance multiplied by total leaf - ! area. gscan_patch is the sum over all cohorts, used to normalize. - ! ---------------------------------------------------------------------------- - - gscan_patch = 0.0_r8 - ccohort=>cpatch%tallest - do while(associated(ccohort)) - ccohort_hydr => ccohort%co_hydr - ccohort_hydr%psi_ag(1) = wrf_plant(leaf_p_media,ccohort%pft)%p%psi_from_th(ccohort_hydr%th_ag(1)) - gscan_patch = gscan_patch + ccohort%g_sb_laweight - ccohort => ccohort%shorter - enddo !cohort - - ! The HLM predicted transpiration flux even though no leaves are present? - if(bc_in(s)%qflx_transp_pa(ifp) > 1.e-10_r8 .and. gscan_patchcpatch%tallest - co_loop1: do while(associated(ccohort)) + gscan_patch = 0.0_r8 + ccohort=>cpatch%tallest + do while(associated(ccohort)) + ccohort_hydr => ccohort%co_hydr + ccohort_hydr%psi_ag(1) = wrf_plant(leaf_p_media,ccohort%pft)%p%psi_from_th(ccohort_hydr%th_ag(1)) + gscan_patch = gscan_patch + ccohort%g_sb_laweight + ccohort => ccohort%shorter + enddo !cohort - ccohort_hydr => ccohort%co_hydr - ft = ccohort%pft + ! The HLM predicted transpiration flux even though no leaves are present? + if(bc_in(s)%qflx_transp_pa(ifp) > 1.e-10_r8 .and. gscan_patchcpatch%tallest + co_loop1: do while(associated(ccohort)) - ! This can cause large transpiration due to small g_sb_laweight - if(ccohort%g_sb_laweight>nearzero) then - qflx_tran_veg_indiv = bc_in(s)%qflx_transp_pa(ifp) * cpatch%total_canopy_area * & - (ccohort%g_sb_laweight/gscan_patch)/ccohort%n - else - qflx_tran_veg_indiv = 0._r8 - end if + ccohort_hydr => ccohort%co_hydr + ft = ccohort%pft - ! Save the transpiration flux for diagnostics (currently its a constant boundary condition) - ccohort_hydr%qtop = qflx_tran_veg_indiv*dtime - - transp_flux = transp_flux + (qflx_tran_veg_indiv*dtime)*ccohort%n*AREA_INV - - ! VERTICAL LAYER CONTRIBUTION TO TOTAL ROOT WATER UPTAKE OR LOSS - ! _____ - ! | | - ! |leaf | - ! |_____| - ! / - ! \ - ! / - ! __\__ - ! | | - ! |stem | - ! |_____| - !------/----------------_____--------------------------------- - ! \ | | | | | | | - ! / _/\/\|aroot| | |shell | shell | shell | layer j-1 - ! \ _/ |_____| | | k-1 | k | k+1 | - !------/------_/--------_____-------------------------------------- - ! \ _/ | | | | | | | - ! __/__ / _/\/\/\/\/|aroot| | | shell | shell | shell | layer j - ! | |_/ |_____| | | k-1 | k | k+1 | - !---|troot|-------------_____---------------------------------------------- - ! |_____|\_ | | | | | | | - ! \/\/\/\/\/|aroot| | | shell | shell | shell | layer j+1 - ! |_____| | | k-1 | k | k+1 | - !--------------------------------------------------------------------------- - - ! This routine will update the theta values for 1 cohort's flow-path - ! from leaf to the current soil layer. This does NOT - ! update cohort%th_* - - if(hydr_solver == hydr_solver_2DNewton) then - - call MatSolve2D(csite_hydr,ccohort,ccohort_hydr, & - dtime,qflx_tran_veg_indiv, & - sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & - dth_layershell_col) - - elseif(hydr_solver == hydr_solver_2DPicard) then - - call PicardSolve2D(csite_hydr,ccohort,ccohort_hydr, & - dtime,qflx_tran_veg_indiv, & - sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & - dth_layershell_col,csite_hydr%num_nodes) - - elseif(hydr_solver == hydr_solver_1DTaylor ) then - - ! --------------------------------------------------------------------------------- - ! Approach: do nlevsoi_hyd sequential solutions to Richards' equation, - ! each of which encompass all plant nodes and soil nodes for a given soil layer j, - ! with the timestep fraction for each layer-specific solution proportional to each - ! layer's contribution to the total root-soil conductance - ! Water potential in plant nodes is updated after each solution - ! As such, the order across soil layers in which the solution is conducted matters. - ! For now, the order proceeds across soil layers in order of decreasing root-soil conductance - ! NET EFFECT: total water removed from plant-soil system remains the same: it - ! sums up to total transpiration (qflx_tran_veg_indiv*dtime) - ! root water uptake in each layer is proportional to each layer's total - ! root length density and soil matric potential - ! root hydraulic redistribution emerges within this sequence when a - ! layers have transporting-to-absorbing root water potential gradients of opposite sign - ! ----------------------------------------------------------------------------------- - - call OrderLayersForSolve1D(csite_hydr, ccohort, ccohort_hydr, ordered, kbg_layer) - - call ImTaylorSolve1D(lat,lon,recruitflag,csite_hydr,ccohort,ccohort_hydr, & - dtime,qflx_tran_veg_indiv,ordered, kbg_layer, & - sapflow,rootuptake(1:nlevrhiz), & - wb_err_plant,dwat_plant, & - dth_layershell_col) + ! Relative transpiration of this cohort from the whole patch + ! Note that g_sb_laweight / gscan_patch is the weighting that gives cohort contribution per area + ! [mm H2O/plant/s] = [mm H2O/ m2 / s] * [m2 / patch] * [cohort/plant] * [patch/cohort] - end if + ! This can cause large transpiration due to small g_sb_laweight + if(ccohort%g_sb_laweight>nearzero) then + qflx_tran_veg_indiv = bc_in(s)%qflx_transp_pa(ifp) * cpatch%total_canopy_area * & + (ccohort%g_sb_laweight/gscan_patch)/ccohort%n + else + qflx_tran_veg_indiv = 0._r8 + end if - ! Remember the error for the cohort - ccohort_hydr%errh2o = ccohort_hydr%errh2o + wb_err_plant + ! Save the transpiration flux for diagnostics (currently its a constant boundary condition) + ccohort_hydr%qtop = qflx_tran_veg_indiv*dtime + + transp_flux = transp_flux + (qflx_tran_veg_indiv*dtime)*ccohort%n*AREA_INV + + ! VERTICAL LAYER CONTRIBUTION TO TOTAL ROOT WATER UPTAKE OR LOSS + ! _____ + ! | | + ! |leaf | + ! |_____| + ! / + ! \ + ! / + ! __\__ + ! | | + ! |stem | + ! |_____| + !------/----------------_____--------------------------------- + ! \ | | | | | | | + ! / _/\/\|aroot| | |shell | shell | shell | layer j-1 + ! \ _/ |_____| | | k-1 | k | k+1 | + !------/------_/--------_____-------------------------------------- + ! \ _/ | | | | | | | + ! __/__ / _/\/\/\/\/|aroot| | | shell | shell | shell | layer j + ! | |_/ |_____| | | k-1 | k | k+1 | + !---|troot|-------------_____---------------------------------------------- + ! |_____|\_ | | | | | | | + ! \/\/\/\/\/|aroot| | | shell | shell | shell | layer j+1 + ! |_____| | | k-1 | k | k+1 | + !--------------------------------------------------------------------------- + + ! This routine will update the theta values for 1 cohort's flow-path + ! from leaf to the current soil layer. This does NOT + ! update cohort%th_* + + if(hydr_solver == hydr_solver_2DNewton) then + + call MatSolve2D(csite_hydr,ccohort,ccohort_hydr, & + dtime,qflx_tran_veg_indiv, & + sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & + dth_layershell_col) + + elseif(hydr_solver == hydr_solver_2DPicard) then + + call PicardSolve2D(csite_hydr,ccohort,ccohort_hydr, & + dtime,qflx_tran_veg_indiv, & + sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & + dth_layershell_col,csite_hydr%num_nodes) + + elseif(hydr_solver == hydr_solver_1DTaylor ) then + + ! --------------------------------------------------------------------------------- + ! Approach: do nlevsoi_hyd sequential solutions to Richards' equation, + ! each of which encompass all plant nodes and soil nodes for a given soil layer j, + ! with the timestep fraction for each layer-specific solution proportional to each + ! layer's contribution to the total root-soil conductance + ! Water potential in plant nodes is updated after each solution + ! As such, the order across soil layers in which the solution is conducted matters. + ! For now, the order proceeds across soil layers in order of decreasing root-soil conductance + ! NET EFFECT: total water removed from plant-soil system remains the same: it + ! sums up to total transpiration (qflx_tran_veg_indiv*dtime) + ! root water uptake in each layer is proportional to each layer's total + ! root length density and soil matric potential + ! root hydraulic redistribution emerges within this sequence when a + ! layers have transporting-to-absorbing root water potential gradients of opposite sign + ! ----------------------------------------------------------------------------------- + + call OrderLayersForSolve1D(csite_hydr, ccohort, ccohort_hydr, ordered, kbg_layer) + + call ImTaylorSolve1D(lat,lon,recruitflag,csite_hydr,ccohort,ccohort_hydr, & + dtime,qflx_tran_veg_indiv,ordered, kbg_layer, & + sapflow,rootuptake(1:nlevrhiz), & + wb_err_plant,dwat_plant, & + dth_layershell_col) - ! Update total error in [kg/m2 ground] - csite_hydr%errh2o_hyd = csite_hydr%errh2o_hyd + wb_err_plant*ccohort%n*AREA_INV + end if + + ! Remember the error for the cohort + ccohort_hydr%errh2o = ccohort_hydr%errh2o + wb_err_plant + + ! Update total error in [kg/m2 ground] + csite_hydr%errh2o_hyd = csite_hydr%errh2o_hyd + wb_err_plant*ccohort%n*AREA_INV - ! Accumulate site level diagnostic of plant water change [kg/m2] - ! (this is zerod) - csite_hydr%dwat_veg = csite_hydr%dwat_veg + dwat_plant*ccohort%n*AREA_INV + ! Accumulate site level diagnostic of plant water change [kg/m2] + ! (this is zerod) + csite_hydr%dwat_veg = csite_hydr%dwat_veg + dwat_plant*ccohort%n*AREA_INV - ! Update total site-level stored plant water [kg/m2] - ! (this is not zerod, but incremented) - csite_hydr%h2oveg = csite_hydr%h2oveg + dwat_plant*ccohort%n*AREA_INV + ! Update total site-level stored plant water [kg/m2] + ! (this is not zerod, but incremented) + csite_hydr%h2oveg = csite_hydr%h2oveg + dwat_plant*ccohort%n*AREA_INV - sc = ccohort%size_class + sc = ccohort%size_class - ! Sapflow diagnostic [kg/ha/s] - csite_hydr%sapflow_scpf(sc,ft) = csite_hydr%sapflow_scpf(sc,ft) + sapflow*ccohort%n/dtime + ! Sapflow diagnostic [kg/ha/s] + csite_hydr%sapflow_scpf(sc,ft) = csite_hydr%sapflow_scpf(sc,ft) + sapflow*ccohort%n/dtime - ! Root uptake per pft x size class, over set layer depths [kg/ha/m/s] - ! These are normalized by depth (in case the desired horizon extends - ! beyond the actual rhizosphere) + ! Root uptake per pft x size class, over set layer depths [kg/ha/m/s] + ! These are normalized by depth (in case the desired horizon extends + ! beyond the actual rhizosphere) - csite_hydr%rootuptake0_scpf(sc,ft) = csite_hydr%rootuptake0_scpf(sc,ft) + & - SumBetweenDepths(csite_hydr,0._r8,0.1_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime + csite_hydr%rootuptake0_scpf(sc,ft) = csite_hydr%rootuptake0_scpf(sc,ft) + & + SumBetweenDepths(csite_hydr,0._r8,0.1_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime - csite_hydr%rootuptake10_scpf(sc,ft) = csite_hydr%rootuptake10_scpf(sc,ft) + & - SumBetweenDepths(csite_hydr,0.1_r8,0.5_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime + csite_hydr%rootuptake10_scpf(sc,ft) = csite_hydr%rootuptake10_scpf(sc,ft) + & + SumBetweenDepths(csite_hydr,0.1_r8,0.5_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime - csite_hydr%rootuptake50_scpf(sc,ft) = csite_hydr%rootuptake50_scpf(sc,ft) + & - SumBetweenDepths(csite_hydr,0.5_r8,1.0_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime + csite_hydr%rootuptake50_scpf(sc,ft) = csite_hydr%rootuptake50_scpf(sc,ft) + & + SumBetweenDepths(csite_hydr,0.5_r8,1.0_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime - csite_hydr%rootuptake100_scpf(sc,ft) = csite_hydr%rootuptake100_scpf(sc,ft) + & - SumBetweenDepths(csite_hydr,1.0_r8,1.e10_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime + csite_hydr%rootuptake100_scpf(sc,ft) = csite_hydr%rootuptake100_scpf(sc,ft) + & + SumBetweenDepths(csite_hydr,1.0_r8,1.e10_r8,rootuptake(1:nlevrhiz))*ccohort%n/dtime - ! --------------------------------------------------------- - ! Update water potential and frac total conductivity - ! of plant compartments - ! --------------------------------------------------------- + ! --------------------------------------------------------- + ! Update water potential and frac total conductivity + ! of plant compartments + ! --------------------------------------------------------- - call UpdatePlantPsiFTCFromTheta(ccohort,csite_hydr) + call UpdatePlantPsiFTCFromTheta(ccohort,csite_hydr) - ccohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(1)) + ccohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(1)) - ccohort => ccohort%shorter - enddo co_loop1 !cohort - endif ! not bareground patch + ccohort => ccohort%shorter + enddo co_loop1 !cohort + cpatch => cpatch%younger enddo !patch diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index d1413f1b15..5ef14eb1fd 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -31,7 +31,6 @@ module FATESPlantRespPhotosynthMod use FatesConstantsMod, only : rgas_J_K_mol use FatesConstantsMod, only : fates_unset_r8 use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : nocomp_bareground use FatesConstantsMod, only : photosynth_acclim_model_none use FatesConstantsMod, only : photosynth_acclim_model_kumarathunge_etal_2019 use FatesInterfaceTypesMod, only : hlm_use_planthydro @@ -347,7 +346,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ifp = 0 currentpatch => sites(s)%oldest_patch do while (associated(currentpatch)) - if_notbare: if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then + ifp = ifp+1 NCL_p = currentPatch%NCL_p @@ -1079,8 +1078,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end if if_filter2 - end if if_notbare - currentPatch => currentPatch%younger end do diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 97bbcb6d54..90fd1ceabd 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -49,7 +49,7 @@ module SFMainMod integer :: write_SF = ifalse ! for debugging logical :: debug = .false. ! for debugging - + ! ====================================================================================== contains @@ -126,8 +126,9 @@ subroutine UpdateFireWeather(currentSite, bc_in) currentPatch => currentSite%oldest_patch ! If the oldest patch is a bareground patch (i.e. nocomp mode is on) use the first vegetated patch - ! for the iofp index (i.e. the next younger patch) - if (currentPatch%nocomp_pft_label == nocomp_bareground) then + ! for the iofp index (i.e. the next younger patch). But, its possible that there is + ! 100% bareground fraction and no other patches. + if (currentPatch%nocomp_pft_label == nocomp_bareground .and. associated(currentPatch%younger) ) then currentPatch => currentPatch%younger endif @@ -173,33 +174,30 @@ subroutine UpdateFuelCharacteristics(currentSite) currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) - if (currentPatch%nocomp_pft_label /= nocomp_bareground) then - - ! calculate live grass [kgC/m2] - call currentPatch%UpdateLiveGrass() - - ! update fuel loading [kgC/m2] - litter => currentPatch%litter(element_pos(carbon12_element)) - call currentPatch%fuel%UpdateLoading(sum(litter%leaf_fines(:)), & - litter%ag_cwd(1), litter%ag_cwd(2), litter%ag_cwd(3), litter%ag_cwd(4), & - currentPatch%livegrass) - - ! sum up fuel classes and calculate fractional loading for each - call currentPatch%fuel%SumLoading() - call currentPatch%fuel%CalculateFractionalLoading() - - ! calculate fuel moisture [m3/m3] - call currentPatch%fuel%UpdateFuelMoisture(SF_val_SAV, SF_val_drying_ratio, & - currentSite%fireWeather) - - ! calculate geometric properties - call currentPatch%fuel%AverageBulkDensity_NoTrunks(SF_val_FBD) - call currentPatch%fuel%AverageSAV_NoTrunks(SF_val_SAV) + ! calculate live grass [kgC/m2] + call currentPatch%UpdateLiveGrass() + + ! update fuel loading [kgC/m2] + litter => currentPatch%litter(element_pos(carbon12_element)) + call currentPatch%fuel%UpdateLoading(sum(litter%leaf_fines(:)), & + litter%ag_cwd(1), litter%ag_cwd(2), litter%ag_cwd(3), litter%ag_cwd(4), & + currentPatch%livegrass) + + ! sum up fuel classes and calculate fractional loading for each + call currentPatch%fuel%SumLoading() + call currentPatch%fuel%CalculateFractionalLoading() + + ! calculate fuel moisture [m3/m3] + call currentPatch%fuel%UpdateFuelMoisture(SF_val_SAV, SF_val_drying_ratio, & + currentSite%fireWeather) + + ! calculate geometric properties + call currentPatch%fuel%AverageBulkDensity_NoTrunks(SF_val_FBD) + call currentPatch%fuel%AverageSAV_NoTrunks(SF_val_SAV) - end if - currentPatch => currentPatch%younger - end do - + currentPatch => currentPatch%younger + end do + end subroutine UpdateFuelCharacteristics !--------------------------------------------------------------------------------------- @@ -236,7 +234,7 @@ subroutine rate_of_spread (currentSite) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground .and. currentPatch%fuel%non_trunk_loading > nearzero) then + if_nontrunk_loading: if(currentPatch%fuel%non_trunk_loading > nearzero) then ! remove mineral content from net fuel load per Thonicke 2010 for ir calculation currentPatch%fuel%non_trunk_loading = currentPatch%fuel%non_trunk_loading * (1.0_r8 - SF_val_miner_total) !net of minerals @@ -346,9 +344,8 @@ subroutine rate_of_spread (currentSite) ! backward ROS wind not changed by vegetation currentPatch%ROS_back = currentPatch%ROS_front*exp(-0.012_r8*currentSite%wind) - end if ! nocomp_pft_label check + end if if_nontrunk_loading currentPatch => currentPatch%younger - enddo !end patch loop end subroutine rate_of_spread @@ -380,8 +377,6 @@ subroutine ground_fuel_consumption ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - currentPatch%fuel%frac_burnt(:) = 1.0_r8 ! Calculate fraction of litter is burnt for all classes. ! Equation B1 in Thonicke et al. 2010--- @@ -443,8 +438,6 @@ subroutine ground_fuel_consumption ( currentSite ) ! ignore 1000hr fuels. Just interested in fuels affecting ROS currentPatch%TFC_ROS = sum(FC_ground)-FC_ground(tr_sf) - end if ! nocomp_pft_label check - currentPatch=>currentPatch%younger; enddo !end patch loop @@ -509,7 +502,7 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) ! If the oldest patch is a bareground patch (i.e. nocomp mode is on) use the first vegetated patch ! for the iofp index (i.e. the next younger patch) - if(currentPatch%nocomp_pft_label .eq. nocomp_bareground)then + if(currentPatch%nocomp_pft_label .eq. nocomp_bareground .and. associated(currentPatch%younger))then currentPatch => currentPatch%younger endif @@ -538,8 +531,6 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) currentPatch => currentSite%oldest_patch; do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - ! ---initialize patch parameters to zero--- currentPatch%FI = 0._r8 currentPatch%fire = 0 @@ -644,10 +635,8 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) endif endif ! NF ignitions check - endif ! nocomp_pft_label check currentPatch => currentPatch%younger - enddo !end patch loop end subroutine area_burnt_intensity @@ -677,8 +666,6 @@ subroutine crown_scorching ( currentSite ) currentPatch => currentSite%oldest_patch; do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - tree_ag_biomass = 0.0_r8 if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; @@ -711,7 +698,6 @@ subroutine crown_scorching ( currentSite ) end do endif !fire - endif !nocomp_pft_label currentPatch => currentPatch%younger; enddo !end patch loop @@ -735,7 +721,6 @@ subroutine crown_damage ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then if (currentPatch%fire == 1) then currentCohort=>currentPatch%tallest @@ -777,7 +762,6 @@ subroutine crown_damage ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check currentPatch => currentPatch%younger; @@ -804,8 +788,6 @@ subroutine cambial_damage_kill ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; do while(associated(currentCohort)) @@ -830,7 +812,6 @@ subroutine cambial_damage_kill ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check currentPatch=>currentPatch%younger; @@ -857,8 +838,6 @@ subroutine post_fire_mortality ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest do while(associated(currentCohort)) @@ -878,8 +857,6 @@ subroutine post_fire_mortality ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check - currentPatch => currentPatch%younger enddo !end patch loop diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 9fc11491c1..81a53b8f52 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -481,42 +481,42 @@ subroutine set_site_properties( nsites, sites,bc_in ) use_fates_luh_if: if (hlm_use_luh .eq. itrue) then ! MAPPING OF FATES PFTs on to HLM_PFTs with land use ! add up the area associated with each FATES PFT - ! where pft_areafrac_lu is the area of land in each HLM PFT and land use type (from surface dataset) - ! hlm_pft_map is the area of that land in each FATES PFT (from param file) - - ! First check for NaNs in bc_in(s)%pft_areafrac_lu. If so, make everything bare ground. - if ( .not. (any( isnan( bc_in(s)%pft_areafrac_lu (:,:) )) .or. isnan( bc_in(s)%baregroundfrac))) then + ! where pft_areafrac_lu is the area of land in each HLM + ! PFT and land use type (from surface dataset) + ! hlm_pft_map is the area of that land in each FATES + ! PFT (from param file) + + ! First check for NaNs in bc_in(s)%pft_areafrac_lu. + ! If so, make everything bare ground. + if ( .not. (any( isnan( bc_in(s)%pft_areafrac_lu (:,:) )) .or. & + isnan( bc_in(s)%baregroundfrac))) then do i_landusetype = 1, n_landuse_cats if (.not. is_crop(i_landusetype)) then do hlm_pft = 1,size( EDPftvarcon_inst%hlm_pft_map,2) do fates_pft = 1,numpft ! loop round all fates pfts for all hlm pfts - sites(s)%area_pft(fates_pft,i_landusetype) = sites(s)%area_pft(fates_pft,i_landusetype) + & - EDPftvarcon_inst%hlm_pft_map(fates_pft,hlm_pft) * bc_in(s)%pft_areafrac_lu(hlm_pft,i_landusetype) + sites(s)%area_pft(fates_pft,i_landusetype) = & + sites(s)%area_pft(fates_pft,i_landusetype) + & + EDPftvarcon_inst%hlm_pft_map(fates_pft,hlm_pft) * & + bc_in(s)%pft_areafrac_lu(hlm_pft,i_landusetype) end do end do !hlm_pft else - ! for crops, we need to use different logic because the bc_in(s)%pft_areafrac_lu() information only exists for natural PFTs + ! for crops, we need to use different logic because + ! the bc_in(s)%pft_areafrac_lu() information only exists for natural PFTs sites(s)%area_pft(crop_lu_pft_vector(i_landusetype),i_landusetype) = 1._r8 endif end do sites(s)%area_bareground = bc_in(s)%baregroundfrac else - !if ( all( isnan( bc_in(s)%pft_areafrac_lu (:,:))) .and. isnan(bc_in(s)%baregroundfrac)) then - ! if given all NaNs, then make everything bare ground - sites(s)%area_bareground = 1._r8 - sites(s)%area_pft(:,:) = 0._r8 - write(fates_log(),*) 'Nan values for pftareafrac. dumping site info.' - call dump_site(sites(s)) - !else - ! ! if only some things are NaN but not all, then something terrible has probably happened. crash. - ! write(fates_log(),*) 'some but, not all, of the data in the PFT by LU matrix at this site is NaN.' - ! write(fates_log(),*) 'recommend checking the dataset to see what has happened.' - ! call endrun(msg=errMsg(sourcefile, __LINE__)) - !endif + + sites(s)%area_bareground = 1._r8 + sites(s)%area_pft(:,:) = 0._r8 + endif - else + else ! use_fates_luh_if + ! MAPPING OF FATES PFTs on to HLM_PFTs ! add up the area associated with each FATES PFT ! where pft_areafrac is the area of land in each HLM PFT and (from surface dataset) @@ -539,7 +539,8 @@ subroutine set_site_properties( nsites, sites,bc_in ) ! remove tiny patches to prevent numerical errors in terminate patches if (sites(s)%area_pft(ft, i_landusetype) .lt. min_nocomp_pftfrac_perlanduse & .and. sites(s)%area_pft(ft, i_landusetype) .gt. nearzero) then - if(debug) write(fates_log(),*) 'removing small numbers in site%area_pft',s,ft,i_landusetype,sites(s)%area_pft(ft, i_landusetype) + if(debug) write(fates_log(),*) 'removing small numbers in site%area_pft', & + s,ft,i_landusetype,sites(s)%area_pft(ft, i_landusetype) sites(s)%area_pft(ft, i_landusetype)=0.0_r8 endif @@ -551,17 +552,24 @@ subroutine set_site_properties( nsites, sites,bc_in ) end do end do - ! if in nocomp mode, and the number of nocomp PFTs of a given land use type is greater than the maximum number of patches - ! allowed to be allocated for that land use type, then only keep the number of PFTs correspondign to the number of patches - ! allowed on that land use type, starting with the PFTs with greatest area coverage and working down + ! if in nocomp mode, and the number of nocomp PFTs of a given + ! land use type is greater than the maximum number of patches + ! allowed to be allocated for that land use type, then only + ! keep the number of PFTs correspondign to the number of patches + ! allowed on that land use type, starting with the PFTs with + ! greatest area coverage and working down if (hlm_use_nocomp .eq. itrue) then do i_landusetype = 1, n_landuse_cats - ! count how many PFTs have areas greater than zero and compare to the number of patches allowed - if (COUNT(sites(s)%area_pft(:, i_landusetype) .gt. 0._r8) > max_nocomp_pfts_by_landuse(i_landusetype)) then + ! count how many PFTs have areas greater than zero and + ! compare to the number of patches allowed + if (COUNT(sites(s)%area_pft(:, i_landusetype) .gt. 0._r8) > & + max_nocomp_pfts_by_landuse(i_landusetype)) then ! write current vector to log file - if(debug) write(fates_log(),*) 'too many PFTs for LU type ', i_landusetype, sites(s)%area_pft(:, i_landusetype) + if(debug) write(fates_log(),*) 'too many PFTs for LU type ', & + i_landusetype, sites(s)%area_pft(:, i_landusetype) - ! start from largest area, put that PFT's area into a temp vector, and then work down to successively smaller-area PFTs, + ! start from largest area, put that PFT's area into a temp vector, + ! and then work down to successively smaller-area PFTs, ! at the end replace the original vector with the temp vector temp_vec(:) = 0._r8 do i_pftcount = 1, max_nocomp_pfts_by_landuse(i_landusetype) @@ -572,7 +580,8 @@ subroutine set_site_properties( nsites, sites,bc_in ) sites(s)%area_pft(:, i_landusetype) = temp_vec(:) ! write adjusted vector to log file - if(debug) write(fates_log(),*) 'new PFT vector for LU type', i_landusetype, sites(s)%area_pft(:, i_landusetype) + if(debug) write(fates_log(),*) 'new PFT vector for LU type', & + i_landusetype, sites(s)%area_pft(:, i_landusetype) endif end do end if diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 01111c0626..406e7f1dcc 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -205,15 +205,8 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) if (hlm_use_ed_st3.eq.ifalse.and.hlm_use_sp.eq.ifalse) then ! Bypass if ST3 - ! Check that the site doesn't consist solely of a single bareground patch. - ! If so, skip the fire model. Since the bareground patch should be the - ! oldest patch per set_patchno, we check that the youngest patch isn't zero. - ! If there are multiple patches on the site, the bareground patch is avoided - ! at the level of the fire_model subroutines. - if (currentSite%youngest_patch%patchno .ne. 0) then - call fire_model(currentSite, bc_in) - end if - + call fire_model(currentSite, bc_in) + ! Calculate disturbance and mortality based on previous timestep vegetation. ! disturbance_rates calls logging mortality and other mortalities, Yi Xu call disturbance_rates(currentSite, bc_in) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index cc906fecef..03cb8b7e7a 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -848,7 +848,7 @@ subroutine FatesReceiveParams(fates_params) data=tmp_vector_by_landuse2) maxpatches_by_landuse(:) = nint(tmp_vector_by_landuse2(:)) - maxpatch_total = sum(maxpatches_by_landuse(:)) + !maxpatch_total = sum(maxpatches_by_landuse(:)) deallocate(tmp_vector_by_landuse2) call fates_params%RetrieveParameterAllocate(name=ED_name_max_nocomp_pfts_by_landuse, & diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 3d91f2f1b9..26d8e8dbaa 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -12,6 +12,7 @@ module EDPftvarcon use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : itrue, ifalse + use FatesConstantsMod, only : primaryland use PRTParametersMod, only : prt_params use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun @@ -1830,7 +1831,7 @@ subroutine FatesCheckParams(is_master) use EDParamsMod , only : radiation_model, dayl_switch use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog,hlm_use_sp, hlm_name use FatesInterfaceTypesMod, only : hlm_use_inventory_init - use FatesInterfaceTypesMod, only : hlm_use_nocomp + use FatesInterfaceTypesMod, only : hlm_use_nocomp,hlm_use_luh use EDParamsMod , only : max_nocomp_pfts_by_landuse, maxpatches_by_landuse use FatesConstantsMod , only : n_landuse_cats @@ -2256,13 +2257,15 @@ subroutine FatesCheckParams(is_master) if ( hlm_use_nocomp .eq. itrue .and. hlm_use_sp.eq.ifalse) then do i_lu = 1, n_landuse_cats - if (max_nocomp_pfts_by_landuse(i_lu) .gt. maxpatches_by_landuse(i_lu)) then - write(fates_log(),*) 'The max number of nocomp PFTs must all be less than or equal to the number of patches, for a given land use type' - write(fates_log(),*) 'land use index:',i_lu - write(fates_log(),*) 'max_nocomp_pfts_by_landuse(i_lu):', max_nocomp_pfts_by_landuse(i_lu) - write(fates_log(),*) 'maxpatches_by_landuse(i_lu):', maxpatches_by_landuse(i_lu) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(i_lu.eq.primaryland .or. hlm_use_luh.eq.itrue)then + if (max_nocomp_pfts_by_landuse(i_lu) .gt. maxpatches_by_landuse(i_lu)) then + write(fates_log(),*) 'The max number of nocomp PFTs must all be less than or equal to the number of patches, for a given land use type' + write(fates_log(),*) 'land use index:',i_lu + write(fates_log(),*) 'max_nocomp_pfts_by_landuse(i_lu):', max_nocomp_pfts_by_landuse(i_lu) + write(fates_log(),*) 'maxpatches_by_landuse(i_lu):', maxpatches_by_landuse(i_lu) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if end do endif diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 55287c40f5..c7be667ef6 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -5072,8 +5072,6 @@ subroutine update_history_hifrq1(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) do while(associated(cpatch)) ipa = ipa + 1 - - hio_c_stomata_si(io_si) = hio_c_stomata_si(io_si) + & cpatch%c_stomata * cpatch%total_canopy_area * mol_per_umol * site_area_veg_inv @@ -5081,11 +5079,9 @@ subroutine update_history_hifrq1(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) cpatch%c_lblayer * cpatch%total_canopy_area * mol_per_umol * site_area_veg_inv ! Only accumulate the instantaneous vegetation temperature for vegetated patches - if (cpatch%patchno .ne. 0) then - hio_tveg(io_si) = hio_tveg(io_si) + & - (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm) * & - cpatch%total_canopy_area * site_area_veg_inv - end if + hio_tveg(io_si) = hio_tveg(io_si) + & + (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm) * & + cpatch%total_canopy_area * site_area_veg_inv ccohort => cpatch%shortest do while(associated(ccohort)) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index e77b3e34bc..793a7b1f5b 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -184,6 +184,8 @@ module FatesInterfaceMod public :: DetermineGridCellNeighbors logical :: debug = .false. ! for debugging this module + + logical, parameter :: preserve_b4b = .true. contains @@ -654,20 +656,17 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) ! Include the bare-ground patch for these patch-level boundary conditions ! (it will always be zero for all of these) - !if(hlm_use_ch4.eq.itrue) then - allocate(bc_out%annavg_agnpp_pa(0:maxpatch_total));bc_out%annavg_agnpp_pa(:)=nan - allocate(bc_out%annavg_bgnpp_pa(0:maxpatch_total));bc_out%annavg_bgnpp_pa(:)=nan - allocate(bc_out%annsum_npp_pa(0:maxpatch_total));bc_out%annsum_npp_pa(:)=nan - allocate(bc_out%frootc_pa(0:maxpatch_total));bc_out%frootc_pa(:)=nan - allocate(bc_out%root_resp(nlevsoil_in));bc_out%root_resp(:)=nan - allocate(bc_out%woody_frac_aere_pa(0:maxpatch_total));bc_out%woody_frac_aere_pa(:)=nan - allocate(bc_out%rootfr_pa(0:maxpatch_total,nlevsoil_in)) - bc_out%rootfr_pa(:,:)=nan - - ! Give the bare-ground root fractions a nominal fraction of unity over depth - bc_out%rootfr_pa(0,1:nlevsoil_in)=1._r8/real(nlevsoil_in,r8) - !end if - + allocate(bc_out%annavg_agnpp_pa(0:maxpatch_total));bc_out%annavg_agnpp_pa(:)=0._r8 + allocate(bc_out%annavg_bgnpp_pa(0:maxpatch_total));bc_out%annavg_bgnpp_pa(:)=0._r8 + allocate(bc_out%annsum_npp_pa(0:maxpatch_total));bc_out%annsum_npp_pa(:)=0._r8 + allocate(bc_out%frootc_pa(0:maxpatch_total));bc_out%frootc_pa(:)=0._r8 + allocate(bc_out%root_resp(nlevsoil_in));bc_out%root_resp(:)=0._r8 + allocate(bc_out%woody_frac_aere_pa(0:maxpatch_total));bc_out%woody_frac_aere_pa(:)=0._r8 + allocate(bc_out%rootfr_pa(0:maxpatch_total,nlevsoil_in)) + bc_out%rootfr_pa(:,:)=0._r8 + + ! Give the bare-ground root fractions a nominal fraction of unity over depth + bc_out%rootfr_pa(0,1:nlevsoil_in)=1._r8/real(nlevsoil_in,r8) bc_out%ema_npp = nan ! Fates -> BGC fragmentation mass fluxes @@ -783,6 +782,14 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade call FatesReadParameters(param_reader) fates_numpft = size(prt_params%wood_density,dim=1) + + ! maxpatch_total: the number of patches that FATES needs to allocate space for + ! + ! fates_maxPatchesPerSite: the number of patches that the HLM needs to allocate + ! space for. Why the fates_ prefix then? Because + ! it is for readability on the host side so it is clear + ! in the code that fates is dictating the allocation need + if(hlm_use_sp==itrue)then @@ -792,7 +799,7 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade maxpatches_by_landuse(primaryland) = fates_numpft maxpatches_by_landuse(secondaryland:n_landuse_cats) = 0 - maxpatch_total = fates_numpft + maxpatch_total = fates_numpft + 1 ! If this is an SP run, we actually need enough patches on the ! CLM/ELM side of the code to hold the LAI data. This @@ -805,24 +812,46 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade else - ! If we are using fixed biogeography or no-comp then we - ! can also apply those constraints to maxpatch_primaryland and secondary - ! and that value will match fates_maxPatchesPerSite - - if(hlm_use_nocomp==itrue) then - - maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) - maxpatch_total = sum(maxpatches_by_landuse(:)) + if_preserve_b4b: if(preserve_b4b)then + if(hlm_use_nocomp==itrue) then + maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) + maxpatch_total = sum(maxpatches_by_landuse(:))+1 + else + maxpatch_total = sum(maxpatches_by_landuse(:)) + end if + else + + ! If we are using fixed biogeography or no-comp then we + ! can also apply those constraints to maxpatch_primaryland and secondary + ! and that value will match fates_maxPatchesPerSite + + if(hlm_use_luh==ifalse) then + maxpatches_by_landuse(secondaryland:n_landuse_cats) = 0 + end if + + if(hlm_use_nocomp==itrue) then + + ! We do not modify maxpatches_by_landuse here, and/or ensure + ! that it is large enough to accomodate all of our pfts, because + ! we will compare it to the max_nocomp_pfts_by_landuse() array to + ! make sure all indices are larger. + + ! We hold the bareground area in one of the patches + maxpatch_total = sum(maxpatches_by_landuse(:))+1 + else + maxpatch_total = sum(maxpatches_by_landuse(:)) + end if + + ! maxpatch_total does not include HLM-side bare ground (so add 1) + ! Note: Yes, there are two bare-ground patch indices needed for no-comp, where + ! index 1 on the fates side will donate all its area to index 0 on the host + ! side, and index 1 on the host side will have zero area. - !if(maxpatch_primary sites(s)%oldest_patch - do while(associated(cpatch)) - if (cpatch%patchno .ne. 0) then + do_patch: do while(associated(cpatch)) ifp=ifp+1 call cpatch%tveg24%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) @@ -2128,14 +2156,14 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ! levels are the light on the exposed ground at the surface ! and the low levels are the intensity under the bottom-most ! vegetation. - + call SeedlingParPatch(cpatch, & bc_in(s)%solad_parb(ifp,ipar) + bc_in(s)%solai_parb(ifp,ipar), & seedling_par_high, par_high_frac, seedling_par_low,& & par_low_frac) - + new_seedling_layer_par = seedling_par_high*par_high_frac + seedling_par_low*par_low_frac - + call cpatch%seedling_layer_par24%UpdateRMean(new_seedling_layer_par) call cpatch%sdlng_mort_par%UpdateRMean(new_seedling_layer_par) call cpatch%sdlng2sap_par%UpdateRMean(new_seedling_layer_par) @@ -2162,7 +2190,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd) enddo !end pft loop - + end if ccohort => cpatch%tallest @@ -2176,10 +2204,9 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ccohort => ccohort%shorter end do - end if cpatch => cpatch%younger - enddo + enddo do_patch ! Smoothed [gc/m2/yr] if(sites(s)%ema_npp<-9000._r8)then diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index ae73748a42..448e84b28c 100644 --- a/radiation/FatesNormanRadMod.F90 +++ b/radiation/FatesNormanRadMod.F90 @@ -16,7 +16,6 @@ module FatesNormanRadMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue use FatesConstantsMod , only : pi_const - use FatesConstantsMod , only : nocomp_bareground use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : bc_out_type use FatesInterfaceTypesMod , only : numpft diff --git a/radiation/FatesRadiationDriveMod.F90 b/radiation/FatesRadiationDriveMod.F90 index b3e36c39b0..c3d2c9c041 100644 --- a/radiation/FatesRadiationDriveMod.F90 +++ b/radiation/FatesRadiationDriveMod.F90 @@ -17,7 +17,6 @@ module FatesRadiationDriveMod use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : itrue use FatesConstantsMod , only : pi_const - use FatesConstantsMod , only : nocomp_bareground use FatesConstantsMod , only : nearzero use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : bc_out_type @@ -101,147 +100,139 @@ subroutine FatesNormalizedCanopyRadiation(nsites, sites, bc_in, bc_out ) currentpatch => sites(s)%oldest_patch do while (associated(currentpatch)) - ! do not do albedo calculations for bare ground patch in SP mode - ! and (more impotantly) do not iterate ifp or it will mess up the indexing wherein - ! ifp=1 is the first vegetated patch. - - if_notbareground: if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then - - ifp = ifp+1 - - ! Zero diagnostics - 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 - currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8 - currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8 - - currentPatch%rad_error(:) = hlm_hio_ignore_val - - 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:num_swb) = bc_in(s)%albgr_dif_rb(1:num_swb) - currentPatch%gnd_alb_dir(1:num_swb) = bc_in(s)%albgr_dir_rb(1:num_swb) - currentPatch%fcansno = bc_in(s)%fcansno_pa(ifp) - - if(radiation_model.eq.twostr_solver) then - call currentPatch%twostr%CanopyPrep(bc_in(s)%fcansno_pa(ifp)) - call currentPatch%twostr%ZenithPrep(bc_in(s)%coszen_pa(ifp)) + ifp = ifp+1 + + ! Zero diagnostics + 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 + currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8 + + currentPatch%rad_error(:) = hlm_hio_ignore_val + + 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:num_swb) = bc_in(s)%albgr_dif_rb(1:num_swb) + currentPatch%gnd_alb_dir(1:num_swb) = bc_in(s)%albgr_dir_rb(1:num_swb) + currentPatch%fcansno = bc_in(s)%fcansno_pa(ifp) + + if(radiation_model.eq.twostr_solver) then + call currentPatch%twostr%CanopyPrep(bc_in(s)%fcansno_pa(ifp)) + call currentPatch%twostr%ZenithPrep(bc_in(s)%coszen_pa(ifp)) + end if + + if_zenith_flag: if(.not.currentPatch%solar_zenith_flag )then + + ! Sun below horizon, trivial solution + ! Note (RGK-MLO): Investigate twilight mechanics for + ! non-zero diffuse radiation when cosz<=0 + + ! Temporarily turn off to preserve b4b + ! preserve_b4b will be removed soon. This is kept here to prevent + ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) + if (.not.preserve_b4b) then + bc_out(s)%albd_parb(ifp,:) = 1._r8 + bc_out(s)%albi_parb(ifp,:) = 1._r8 + bc_out(s)%fabi_parb(ifp,:) = 0._r8 + bc_out(s)%fabd_parb(ifp,:) = 0._r8 + bc_out(s)%ftdd_parb(ifp,:) = 0._r8 + bc_out(s)%ftid_parb(ifp,:) = 0._r8 + bc_out(s)%ftii_parb(ifp,:) = 0._r8 end if + else + + 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_nrad: if (maxval(currentPatch%nrad(1,:))==0)then + ! there are no leaf layers in this patch. it is effectively bare ground. + bc_out(s)%fabd_parb(ifp,:) = 0.0_r8 + bc_out(s)%fabi_parb(ifp,:) = 0.0_r8 + currentPatch%rad_error(:) = 0.0_r8 + + do ib = 1,num_swb + bc_out(s)%albd_parb(ifp,ib) = bc_in(s)%albgr_dir_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) = 0.0_r8 + bc_out(s)%ftii_parb(ifp,ib) = 1.0_r8 + enddo - if_zenith_flag: if(.not.currentPatch%solar_zenith_flag )then - - ! Sun below horizon, trivial solution - ! Note (RGK-MLO): Investigate twilight mechanics for - ! non-zero diffuse radiation when cosz<=0 - - ! Temporarily turn off to preserve b4b - ! preserve_b4b will be removed soon. This is kept here to prevent - ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) - if (.not.preserve_b4b) then - bc_out(s)%albd_parb(ifp,:) = 1._r8 - bc_out(s)%albi_parb(ifp,:) = 1._r8 - bc_out(s)%fabi_parb(ifp,:) = 0._r8 - bc_out(s)%fabd_parb(ifp,:) = 0._r8 - bc_out(s)%ftdd_parb(ifp,:) = 0._r8 - bc_out(s)%ftid_parb(ifp,:) = 0._r8 - bc_out(s)%ftii_parb(ifp,:) = 0._r8 - end if else - 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_nrad: if (maxval(currentPatch%nrad(1,:))==0)then - ! there are no leaf layers in this patch. it is effectively bare ground. - bc_out(s)%fabd_parb(ifp,:) = 0.0_r8 - bc_out(s)%fabi_parb(ifp,:) = 0.0_r8 - currentPatch%rad_error(:) = 0.0_r8 - - do ib = 1,num_swb - bc_out(s)%albd_parb(ifp,ib) = bc_in(s)%albgr_dir_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) = 0.0_r8 - bc_out(s)%ftii_parb(ifp,ib) = 1.0_r8 - enddo - - else - - select case(radiation_model) - case(norman_solver) - - call PatchNormanRadiation (currentPatch, & - bc_out(s)%albd_parb(ifp,:), & ! Surface Albedo direct - bc_out(s)%albi_parb(ifp,:), & ! Surface Albedo (indirect) diffuse - bc_out(s)%fabd_parb(ifp,:), & ! Fraction direct absorbed by canopy per unit incident - bc_out(s)%fabi_parb(ifp,:), & ! Fraction diffuse absorbed by canopy per unit incident - bc_out(s)%ftdd_parb(ifp,:), & ! Down direct flux below canopy per unit direct at top - bc_out(s)%ftid_parb(ifp,:), & ! Down diffuse flux below canopy per unit direct at top - bc_out(s)%ftii_parb(ifp,:)) ! Down diffuse flux below canopy per unit diffuse at top - - case(twostr_solver) - - associate( twostr => currentPatch%twostr) - - !call twostr%CanopyPrep(bc_in(s)%fcansno_pa(ifp)) - !call twostr%ZenithPrep(bc_in(s)%coszen_pa(ifp)) - - do ib = 1,num_swb - - twostr%band(ib)%albedo_grnd_diff = bc_in(s)%albgr_dif_rb(ib) - twostr%band(ib)%albedo_grnd_beam = bc_in(s)%albgr_dir_rb(ib) - - call twostr%Solve(ib, & ! in - normalized_upper_boundary, & ! in - 1.0_r8,1.0_r8, & ! in - sites(s)%taulambda_2str, & ! inout (scratch) - sites(s)%omega_2str, & ! inout (scratch) - sites(s)%ipiv_2str, & ! inout (scratch) - bc_out(s)%albd_parb(ifp,ib), & ! out - bc_out(s)%albi_parb(ifp,ib), & ! out - currentPatch%rad_error(ib), & ! out - bc_out(s)%fabd_parb(ifp,ib), & ! out - bc_out(s)%fabi_parb(ifp,ib), & ! out - bc_out(s)%ftdd_parb(ifp,ib), & ! out - bc_out(s)%ftid_parb(ifp,ib), & ! out - bc_out(s)%ftii_parb(ifp,ib)) - - if(debug) then - currentPatch%twostr%band(ib)%Rbeam_atm = 1._r8 - currentPatch%twostr%band(ib)%Rdiff_atm = 1._r8 - call CheckPatchRadiationBalance(currentPatch, sites(s)%snow_depth, & - ib, bc_out(s)%fabd_parb(ifp,ib),bc_out(s)%fabi_parb(ifp,ib)) - currentPatch%twostr%band(ib)%Rbeam_atm = fates_unset_r8 - currentPatch%twostr%band(ib)%Rdiff_atm = fates_unset_r8 - - if(bc_out(s)%fabi_parb(ifp,ib)>1.0 .or. bc_out(s)%fabd_parb(ifp,ib)>1.0)then - write(fates_log(),*) 'absorbed fraction > 1.0?' - write(fates_log(),*) ifp,ib,bc_out(s)%fabi_parb(ifp,ib),bc_out(s)%fabd_parb(ifp,ib) - call twostr%Dump(ib,lat=sites(s)%lat,lon=sites(s)%lon) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + select case(radiation_model) + case(norman_solver) + + call PatchNormanRadiation (currentPatch, & + bc_out(s)%albd_parb(ifp,:), & ! Surface Albedo direct + bc_out(s)%albi_parb(ifp,:), & ! Surface Albedo (indirect) diffuse + bc_out(s)%fabd_parb(ifp,:), & ! Fraction direct absorbed by canopy per unit incident + bc_out(s)%fabi_parb(ifp,:), & ! Fraction diffuse absorbed by canopy per unit incident + bc_out(s)%ftdd_parb(ifp,:), & ! Down direct flux below canopy per unit direct at top + bc_out(s)%ftid_parb(ifp,:), & ! Down diffuse flux below canopy per unit direct at top + bc_out(s)%ftii_parb(ifp,:)) ! Down diffuse flux below canopy per unit diffuse at top + + case(twostr_solver) + + associate( twostr => currentPatch%twostr) + + !call twostr%CanopyPrep(bc_in(s)%fcansno_pa(ifp)) + !call twostr%ZenithPrep(bc_in(s)%coszen_pa(ifp)) + + do ib = 1,num_swb + + twostr%band(ib)%albedo_grnd_diff = bc_in(s)%albgr_dif_rb(ib) + twostr%band(ib)%albedo_grnd_beam = bc_in(s)%albgr_dir_rb(ib) + + call twostr%Solve(ib, & ! in + normalized_upper_boundary, & ! in + 1.0_r8,1.0_r8, & ! in + sites(s)%taulambda_2str, & ! inout (scratch) + sites(s)%omega_2str, & ! inout (scratch) + sites(s)%ipiv_2str, & ! inout (scratch) + bc_out(s)%albd_parb(ifp,ib), & ! out + bc_out(s)%albi_parb(ifp,ib), & ! out + currentPatch%rad_error(ib), & ! out + bc_out(s)%fabd_parb(ifp,ib), & ! out + bc_out(s)%fabi_parb(ifp,ib), & ! out + bc_out(s)%ftdd_parb(ifp,ib), & ! out + bc_out(s)%ftid_parb(ifp,ib), & ! out + bc_out(s)%ftii_parb(ifp,ib)) + + if(debug) then + currentPatch%twostr%band(ib)%Rbeam_atm = 1._r8 + currentPatch%twostr%band(ib)%Rdiff_atm = 1._r8 + call CheckPatchRadiationBalance(currentPatch, sites(s)%snow_depth, & + ib, bc_out(s)%fabd_parb(ifp,ib),bc_out(s)%fabi_parb(ifp,ib)) + currentPatch%twostr%band(ib)%Rbeam_atm = fates_unset_r8 + currentPatch%twostr%band(ib)%Rdiff_atm = fates_unset_r8 + + if(bc_out(s)%fabi_parb(ifp,ib)>1.0 .or. bc_out(s)%fabd_parb(ifp,ib)>1.0)then + write(fates_log(),*) 'absorbed fraction > 1.0?' + write(fates_log(),*) ifp,ib,bc_out(s)%fabi_parb(ifp,ib),bc_out(s)%fabd_parb(ifp,ib) + call twostr%Dump(ib,lat=sites(s)%lat,lon=sites(s)%lon) + call endrun(msg=errMsg(sourcefile, __LINE__)) end if + end if - end do - end associate + end do + end associate - end select - - end if if_nrad + end select - endif if_zenith_flag - end if if_notbareground + end if if_nrad + endif if_zenith_flag currentPatch => currentPatch%younger end do ! Loop linked-list patches enddo ! Loop Sites @@ -287,218 +278,209 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) do while (associated(cpatch)) - if_notbareground:if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then !only for veg patches - ! do not do albedo calculations for bare ground patch in SP mode - ! and (more impotantly) do not iterate ifp or it will mess up the indexing wherein - ! ifp=1 is the first vegetated patch. - ifp=ifp+1 - - ! Initialize diagnostics - cpatch%ed_parsun_z(:,:,:) = 0._r8 - cpatch%ed_parsha_z(:,:,:) = 0._r8 - cpatch%ed_laisun_z(:,:,:) = 0._r8 - cpatch%ed_laisha_z(:,:,:) = 0._r8 - cpatch%parprof_pft_dir_z(:,:,:) = 0._r8 - cpatch%parprof_pft_dif_z(:,:,:) = 0._r8 - - bc_out(s)%fsun_pa(ifp) = 0._r8 - + ifp=ifp+1 + + ! Initialize diagnostics + cpatch%ed_parsun_z(:,:,:) = 0._r8 + cpatch%ed_parsha_z(:,:,:) = 0._r8 + cpatch%ed_laisun_z(:,:,:) = 0._r8 + cpatch%ed_laisha_z(:,:,:) = 0._r8 + cpatch%parprof_pft_dir_z(:,:,:) = 0._r8 + cpatch%parprof_pft_dif_z(:,:,:) = 0._r8 + + bc_out(s)%fsun_pa(ifp) = 0._r8 + + ! preserve_b4b will be removed soon. This is kept here to prevent + ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) + if(.not.preserve_b4b)then + bc_out(s)%laisun_pa(ifp) = 0._r8 + bc_out(s)%laisha_pa(ifp) = calc_areaindex(cpatch,'elai') + end if + + sunlai = 0._r8 + shalai = 0._r8 + if_norm_twostr: if (radiation_model.eq.norman_solver) then + + ! Loop over patches to calculate laisun_z and laisha_z for each layer. + ! Derive canopy laisun, laisha, and fsun from layer sums. + ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from + ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value. + + ! cpatch%f_sun is calculated in the surface_albedo routine... + + do cl = 1, cpatch%ncl_p + do ft = 1,numpft + ! preserve_b4b will be removed soon. This is kept here to prevent + ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) + if(.not.preserve_b4b) then + sunlai = sunlai + sum(cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & + cpatch%f_sun(cl,ft,1:cpatch%nrad(cl,ft))) + shalai = shalai + sum(cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft))) + else + do iv = 1,cpatch%nrad(cl,ft) + cpatch%ed_laisun_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & + cpatch%f_sun(CL,ft,iv) + + cpatch%ed_laisha_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & + (1._r8 - cpatch%f_sun(CL,ft,iv)) + + end do + + !needed for the VOC emissions, etc. + sunlai = sunlai + sum(cpatch%ed_laisun_z(CL,ft,1:cpatch%nrad(CL,ft))) + shalai = shalai + sum(cpatch%ed_laisha_z(CL,ft,1:cpatch%nrad(CL,ft))) + + end if + end do + end do ! preserve_b4b will be removed soon. This is kept here to prevent ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) if(.not.preserve_b4b)then - bc_out(s)%laisun_pa(ifp) = 0._r8 - bc_out(s)%laisha_pa(ifp) = calc_areaindex(cpatch,'elai') + shalai = shalai-sunlai end if - sunlai = 0._r8 - shalai = 0._r8 - if_norm_twostr: if (radiation_model.eq.norman_solver) then - - ! Loop over patches to calculate laisun_z and laisha_z for each layer. - ! Derive canopy laisun, laisha, and fsun from layer sums. - ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from - ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value. - - ! cpatch%f_sun is calculated in the surface_albedo routine... - - do cl = 1, cpatch%ncl_p - do ft = 1,numpft - ! preserve_b4b will be removed soon. This is kept here to prevent - ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) - if(.not.preserve_b4b) then - sunlai = sunlai + sum(cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & - cpatch%f_sun(cl,ft,1:cpatch%nrad(cl,ft))) - shalai = shalai + sum(cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft))) - else - do iv = 1,cpatch%nrad(cl,ft) - cpatch%ed_laisun_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & - cpatch%f_sun(CL,ft,iv) - - cpatch%ed_laisha_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & - (1._r8 - cpatch%f_sun(CL,ft,iv)) - - end do - - !needed for the VOC emissions, etc. - sunlai = sunlai + sum(cpatch%ed_laisun_z(CL,ft,1:cpatch%nrad(CL,ft))) - shalai = shalai + sum(cpatch%ed_laisha_z(CL,ft,1:cpatch%nrad(CL,ft))) - - end if - end do - end do - ! preserve_b4b will be removed soon. This is kept here to prevent - ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) - if(.not.preserve_b4b)then - shalai = shalai-sunlai - end if - - if(sunlai+shalai > 0._r8)then - bc_out(s)%fsun_pa(ifp) = sunlai / (sunlai+shalai) - else - bc_out(s)%fsun_pa(ifp) = 0._r8 - endif - - if(debug)then - if(bc_out(s)%fsun_pa(ifp) > 1._r8)then - write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & - sunlai,shalai - endif - end if - - elai = calc_areaindex(cpatch,'elai') - - bc_out(s)%laisun_pa(ifp) = elai*bc_out(s)%fsun_pa(ifp) - bc_out(s)%laisha_pa(ifp) = elai*(1.0_r8-bc_out(s)%fsun_pa(ifp)) - - ! Absorbed PAR profile through canopy - ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo - ! are canopy integrated so that layer values equal big leaf values. - - do cl = 1, cpatch%ncl_p - do ft = 1,numpft - do iv = 1, cpatch%nrad(cl,ft) - - cpatch%ed_parsun_z(cl,ft,iv) = & - bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sun_z(cl,ft,iv) + & - bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sun_z(cl,ft,iv) - - cpatch%ed_parsha_z(cl,ft,iv) = & - bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sha_z(cl,ft,iv) + & - bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sha_z(cl,ft,iv) - - end do !iv - end do !ft - end do !cl - - ! Convert normalized radiation error units from fraction of radiation to W/m2 - do ib = 1,num_swb - cpatch%rad_error(ib) = cpatch%rad_error(ib) * & - (bc_in(s)%solad_parb(ifp,ib) + bc_in(s)%solai_parb(ifp,ib)) - end do - - ! output the actual PAR profiles through the canopy for diagnostic purposes - do cl = 1, cpatch%ncl_p - do ft = 1,numpft - do iv = 1, cpatch%nrad(cl,ft) - cpatch%parprof_pft_dir_z(cl,ft,iv) = (bc_in(s)%solad_parb(ifp,ipar) * & - cpatch%nrmlzd_parprof_pft_dir_z(idirect,cl,ft,iv)) + & - (bc_in(s)%solai_parb(ifp,ipar) * & - cpatch%nrmlzd_parprof_pft_dir_z(idiffuse,cl,ft,iv)) - - cpatch%parprof_pft_dif_z(cl,ft,iv) = (bc_in(s)%solad_parb(ifp,ipar) * & - cpatch%nrmlzd_parprof_pft_dif_z(idirect,cl,ft,iv)) + & - (bc_in(s)%solai_parb(ifp,ipar) * & - cpatch%nrmlzd_parprof_pft_dif_z(idiffuse,cl,ft,iv)) - - end do ! iv - end do ! ft - end do ! cl - + if(sunlai+shalai > 0._r8)then + bc_out(s)%fsun_pa(ifp) = sunlai / (sunlai+shalai) else + bc_out(s)%fsun_pa(ifp) = 0._r8 + endif + + if(debug)then + if(bc_out(s)%fsun_pa(ifp) > 1._r8)then + write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & + sunlai,shalai + endif + end if - ! If there is no sun out, we have a trivial solution - if_zenithflag: if(cpatch%solar_zenith_flag ) then - - ! Two-stream - ! ----------------------------------------------------------- - do ib = 1,num_swb - cpatch%twostr%band(ib)%Rbeam_atm = bc_in(s)%solad_parb(ifp,ib) - cpatch%twostr%band(ib)%Rdiff_atm = bc_in(s)%solai_parb(ifp,ib) - end do - - area_vlpfcl(:,:,:) = 0._r8 - cpatch%f_sun(:,:,:) = 0._r8 - - call FatesPatchFSun(cpatch, & - bc_out(s)%fsun_pa(ifp), & - bc_out(s)%laisun_pa(ifp), & - bc_out(s)%laisha_pa(ifp)) - - associate(twostr => cpatch%twostr) - - do_cl: do cl = 1,twostr%n_lyr - do_icol: do icol = 1,twostr%n_col(cl) - - ft = twostr%scelg(cl,icol)%pft - if_notair: if (ft>0) then - area_frac = twostr%scelg(cl,icol)%area - vai = twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai - nv = minloc(dlower_vai, DIM=1, MASK=(dlower_vai>vai)) - do iv = 1, nv - - vai_top = dlower_vai(iv)-dinc_vai(iv) - vai_bot = min(dlower_vai(iv),twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai) - - cpatch%parprof_pft_dir_z(cl,ft,iv) = cpatch%parprof_pft_dir_z(cl,ft,iv) + & - area_frac*twostr%GetRb(cl,icol,ivis,vai_top) - cpatch%parprof_pft_dif_z(cl,ft,iv) = cpatch%parprof_pft_dif_z(cl,ft,iv) + & - area_frac*twostr%GetRdDn(cl,icol,ivis,vai_top) + & - area_frac*twostr%GetRdUp(cl,icol,ivis,vai_top) - - call twostr%GetAbsRad(cl,icol,ipar,vai_top,vai_bot, & - Rb_abs,Rd_abs,Rd_abs_leaf,Rb_abs_leaf,R_abs_stem,R_abs_snow,leaf_sun_frac) - - cpatch%f_sun(cl,ft,iv) = cpatch%f_sun(cl,ft,iv) + & - area_frac*leaf_sun_frac - cpatch%ed_parsun_z(cl,ft,iv) = cpatch%ed_parsun_z(cl,ft,iv) + & - area_frac*(rd_abs_leaf*leaf_sun_frac + rb_abs_leaf) - cpatch%ed_parsha_z(cl,ft,iv) = cpatch%ed_parsha_z(cl,ft,iv) + & - area_frac*rd_abs_leaf*(1._r8-leaf_sun_frac) - - area_vlpfcl(iv,ft,cl) = area_vlpfcl(iv,ft,cl) + area_frac - end do - end if if_notair - end do do_icol - - do ft = 1,numpft - do_iv: do iv = 1,cpatch%nleaf(cl,ft) - if(area_vlpfcl(iv,ft,cl) cpatch%younger - enddo + ! Absorbed PAR profile through canopy + ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo + ! are canopy integrated so that layer values equal big leaf values. + + do cl = 1, cpatch%ncl_p + do ft = 1,numpft + do iv = 1, cpatch%nrad(cl,ft) + + cpatch%ed_parsun_z(cl,ft,iv) = & + bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sun_z(cl,ft,iv) + & + bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sun_z(cl,ft,iv) + + cpatch%ed_parsha_z(cl,ft,iv) = & + bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sha_z(cl,ft,iv) + & + bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sha_z(cl,ft,iv) + + end do !iv + end do !ft + end do !cl + + ! Convert normalized radiation error units from fraction of radiation to W/m2 + do ib = 1,num_swb + cpatch%rad_error(ib) = cpatch%rad_error(ib) * & + (bc_in(s)%solad_parb(ifp,ib) + bc_in(s)%solai_parb(ifp,ib)) + end do + ! output the actual PAR profiles through the canopy for diagnostic purposes + do cl = 1, cpatch%ncl_p + do ft = 1,numpft + do iv = 1, cpatch%nrad(cl,ft) + cpatch%parprof_pft_dir_z(cl,ft,iv) = (bc_in(s)%solad_parb(ifp,ipar) * & + cpatch%nrmlzd_parprof_pft_dir_z(idirect,cl,ft,iv)) + & + (bc_in(s)%solai_parb(ifp,ipar) * & + cpatch%nrmlzd_parprof_pft_dir_z(idiffuse,cl,ft,iv)) + cpatch%parprof_pft_dif_z(cl,ft,iv) = (bc_in(s)%solad_parb(ifp,ipar) * & + cpatch%nrmlzd_parprof_pft_dif_z(idirect,cl,ft,iv)) + & + (bc_in(s)%solai_parb(ifp,ipar) * & + cpatch%nrmlzd_parprof_pft_dif_z(idiffuse,cl,ft,iv)) + + end do ! iv + end do ! ft + end do ! cl + + else + + ! If there is no sun out, we have a trivial solution + if_zenithflag: if(cpatch%solar_zenith_flag ) then + + ! Two-stream + ! ----------------------------------------------------------- + do ib = 1,num_swb + cpatch%twostr%band(ib)%Rbeam_atm = bc_in(s)%solad_parb(ifp,ib) + cpatch%twostr%band(ib)%Rdiff_atm = bc_in(s)%solai_parb(ifp,ib) + end do + + area_vlpfcl(:,:,:) = 0._r8 + cpatch%f_sun(:,:,:) = 0._r8 + + call FatesPatchFSun(cpatch, & + bc_out(s)%fsun_pa(ifp), & + bc_out(s)%laisun_pa(ifp), & + bc_out(s)%laisha_pa(ifp)) + + associate(twostr => cpatch%twostr) + + do_cl: do cl = 1,twostr%n_lyr + do_icol: do icol = 1,twostr%n_col(cl) + + ft = twostr%scelg(cl,icol)%pft + if_notair: if (ft>0) then + area_frac = twostr%scelg(cl,icol)%area + vai = twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai + nv = minloc(dlower_vai, DIM=1, MASK=(dlower_vai>vai)) + do iv = 1, nv + + vai_top = dlower_vai(iv)-dinc_vai(iv) + vai_bot = min(dlower_vai(iv),twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai) + + cpatch%parprof_pft_dir_z(cl,ft,iv) = cpatch%parprof_pft_dir_z(cl,ft,iv) + & + area_frac*twostr%GetRb(cl,icol,ivis,vai_top) + cpatch%parprof_pft_dif_z(cl,ft,iv) = cpatch%parprof_pft_dif_z(cl,ft,iv) + & + area_frac*twostr%GetRdDn(cl,icol,ivis,vai_top) + & + area_frac*twostr%GetRdUp(cl,icol,ivis,vai_top) + + call twostr%GetAbsRad(cl,icol,ipar,vai_top,vai_bot, & + Rb_abs,Rd_abs,Rd_abs_leaf,Rb_abs_leaf,R_abs_stem,R_abs_snow,leaf_sun_frac) + + cpatch%f_sun(cl,ft,iv) = cpatch%f_sun(cl,ft,iv) + & + area_frac*leaf_sun_frac + cpatch%ed_parsun_z(cl,ft,iv) = cpatch%ed_parsun_z(cl,ft,iv) + & + area_frac*(rd_abs_leaf*leaf_sun_frac + rb_abs_leaf) + cpatch%ed_parsha_z(cl,ft,iv) = cpatch%ed_parsha_z(cl,ft,iv) + & + area_frac*rd_abs_leaf*(1._r8-leaf_sun_frac) + + area_vlpfcl(iv,ft,cl) = area_vlpfcl(iv,ft,cl) + area_frac + end do + end if if_notair + end do do_icol + + do ft = 1,numpft + do_iv: do iv = 1,cpatch%nleaf(cl,ft) + if(area_vlpfcl(iv,ft,cl) cpatch%younger + enddo enddo + return - end subroutine FatesSunShadeFracs end module FatesRadiationDriveMod