From c32431ff5b7f556f16beb9300c7ce5128ead0a36 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 16 May 2024 12:13:31 -0400 Subject: [PATCH 01/11] Removed index filtering of nocomp_bareground patches --- biogeochem/EDCanopyStructureMod.F90 | 195 ++- biogeochem/EDPatchDynamicsMod.F90 | 16 - biogeochem/EDPhysiologyMod.F90 | 57 +- biogeophys/EDAccumulateFluxesMod.F90 | 83 +- biogeophys/EDBtranMod.F90 | 235 ++-- biogeophys/FatesPlantHydraulicsMod.F90 | 302 +++-- biogeophys/FatesPlantRespPhotosynthMod.F90 | 1301 ++++++++++---------- fire/SFMainMod.F90 | 52 +- radiation/FatesNormanRadMod.F90 | 1 - radiation/FatesRadiationDriveMod.F90 | 461 +++---- 10 files changed, 1322 insertions(+), 1381 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 3fbd7b78bd..18b0968856 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -49,6 +49,7 @@ module EDCanopyStructureMod use PRTGenericMod, only : carbon12_element use FatesTwoStreamUtilsMod, only : FatesConstructRadElements use FatesRadiationMemMod , only : twostr_solver + use FatesInterfaceTypesMod, only : hlm_hio_ignore_val ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -1891,127 +1892,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 @@ -2033,11 +2027,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 @@ -2282,7 +2273,7 @@ function NumPotentialCanopyLayers(currentPatch,site_spread,include_substory) res real(r8) :: c_area real(r8) :: arealayer - z = 1 + z = 0 currentCohort => currentPatch%tallest do while (associated(currentCohort)) z = max(z,currentCohort%canopy_layer) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 0a0ac19e3a..8598371d14 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1834,22 +1834,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 0f080127c9..df12cb3b19 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 @@ -3156,44 +3155,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/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 6a1bb001c0..5d2665b12d 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -13,8 +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 private @@ -70,50 +68,49 @@ 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 - - if ( debug ) then - write(fates_log(),*) 'EDAccumFlux 64 ',ccohort%npp_tstep - write(fates_log(),*) 'EDAccumFlux 66 ',ccohort%gpp_tstep - write(fates_log(),*) 'EDAccumFlux 67 ',ccohort%resp_tstep + ifp = ifp+1 - endif - - ccohort%npp_acc = ccohort%npp_acc + ccohort%npp_tstep - ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep - ccohort%resp_acc = ccohort%resp_acc + ccohort%resp_tstep - - ccohort%sym_nfix_daily = ccohort%sym_nfix_daily + ccohort%sym_nfix_tstep + 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 + + if ( debug ) then + + write(fates_log(),*) 'EDAccumFlux 64 ',ccohort%npp_tstep + write(fates_log(),*) 'EDAccumFlux 66 ',ccohort%gpp_tstep + write(fates_log(),*) 'EDAccumFlux 67 ',ccohort%resp_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 + endif + + ccohort%npp_acc = ccohort%npp_acc + ccohort%npp_tstep + ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep + ccohort%resp_acc = ccohort%resp_acc + ccohort%resp_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 0aad6eb977..9837bf6b43 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 @@ -344,761 +343,759 @@ 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 - - ! Part I. Zero output boundary conditions - ! --------------------------------------------------------------------------- - bc_out(s)%rssun_pa(ifp) = 0._r8 - bc_out(s)%rssha_pa(ifp) = 0._r8 - - g_sb_leaves = 0._r8 - patch_la = 0._r8 - - ! Part II. Filter out patches - ! Patch level filter flag for photosynthesis calculations - ! has a short memory, flags: - ! 1 = patch has not been called - ! 2 = patch is currently marked for photosynthesis - ! 3 = patch has been called for photosynthesis already - ! --------------------------------------------------------------------------- - if_filter2: if(bc_in(s)%filter_photo_pa(ifp)==2)then - - - ! Part III. Calculate the number of sublayers for each pft and layer. - ! And then identify which layer/pft combinations have things in them. - ! Output: - ! currentPatch%ncan(:,:) - ! currentPatch%canopy_mask(:,:) - call UpdateCanopyNCanNRadPresent(currentPatch) - - - ! Part IV. Identify some environmentally derived parameters: - ! These quantities are biologically irrelevant - ! Michaelis-Menten constant for CO2 (Pa) - ! Michaelis-Menten constant for O2 (Pa) - ! CO2 compensation point (Pa) - ! leaf boundary layer conductance of h20 - ! constrained vapor pressure - - call GetCanopyGasParameters(bc_in(s)%forc_pbot, & ! in - bc_in(s)%oair_pa(ifp), & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%tgcm_pa(ifp), & ! in - bc_in(s)%eair_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - bc_in(s)%rb_pa(ifp), & ! in - mm_kco2, & ! out - mm_ko2, & ! out - co2_cpoint, & ! out - cf, & ! out - gb_mol, & ! out - ceair) ! out - - ! ------------------------------------------------------------------------ - ! Part VI: Loop over all leaf layers. - ! The concept of leaf layers is a result of the radiative transfer scheme. - ! A leaf layer has uniform radiation environment. Leaf layers are a group - ! of vegetation surfaces (stems and leaves) which inhabit the same - ! canopy-layer "CL", have the same functional type "ft" and within those - ! two partitions are further partitioned into vertical layers where - ! downwelling radiation attenuates in order. - ! In this phase we loop over the leaf layers and calculate the - ! photosynthesis and respiration of the layer (since all biophysical - ! properties are homogeneous). After this step, we can loop through - ! our cohort list, associate each cohort with its list of leaf-layers - ! and transfer these quantities to the cohort. - ! With plant hydraulics, we must realize that photosynthesis and - ! respiration will be different for leaves of each cohort in the leaf - ! layers, as they will have there own hydraulic limitations. - ! NOTE: Only need to flush mask on the number of used pfts, not the whole - ! scratch space. - ! ------------------------------------------------------------------------ - rate_mask_z(:,1:numpft,:) = .false. - - if_any_cohorts: if(currentPatch%countcohorts > 0.0)then - currentCohort => currentPatch%tallest - do_cohort_drive: do while (associated(currentCohort)) ! Cohort loop - - ! Identify the canopy layer (cl), functional type (ft) - ! and the leaf layer (IV) for this cohort - ft = currentCohort%pft - cl = currentCohort%canopy_layer - - ! Calculate the cohort specific elai profile - ! And the top and bottom edges of the veg area index - ! of each layer bin are. Note, if the layers - ! sink below the ground snow line, then the effective - ! LAI and SAI start to shrink to zero, as well as - ! the difference between vaitop and vaibot. - if(currentCohort%treesai>0._r8)then - do iv = 1,currentCohort%nv - call VegAreaLayer(currentCohort%treelai, & - currentCohort%treesai, & - currentCohort%height, & - iv, & - currentCohort%nv, & - currentCohort%pft, & - sites(s)%snow_depth, & - cohort_vaitop(iv), & - cohort_vaibot(iv), & - cohort_layer_elai(iv), & - cohort_layer_esai(iv)) - end do - - cohort_elai = sum(cohort_layer_elai(1:currentCohort%nv)) - cohort_esai = sum(cohort_layer_esai(1:currentCohort%nv)) + ifp = ifp+1 + NCL_p = currentPatch%NCL_p + + ! Part I. Zero output boundary conditions + ! --------------------------------------------------------------------------- + bc_out(s)%rssun_pa(ifp) = 0._r8 + bc_out(s)%rssha_pa(ifp) = 0._r8 + + g_sb_leaves = 0._r8 + patch_la = 0._r8 + + ! Part II. Filter out patches + ! Patch level filter flag for photosynthesis calculations + ! has a short memory, flags: + ! 1 = patch has not been called + ! 2 = patch is currently marked for photosynthesis + ! 3 = patch has been called for photosynthesis already + ! --------------------------------------------------------------------------- + if_filter2: if(bc_in(s)%filter_photo_pa(ifp)==2)then + + + ! Part III. Calculate the number of sublayers for each pft and layer. + ! And then identify which layer/pft combinations have things in them. + ! Output: + ! currentPatch%ncan(:,:) + ! currentPatch%canopy_mask(:,:) + call UpdateCanopyNCanNRadPresent(currentPatch) + + + ! Part IV. Identify some environmentally derived parameters: + ! These quantities are biologically irrelevant + ! Michaelis-Menten constant for CO2 (Pa) + ! Michaelis-Menten constant for O2 (Pa) + ! CO2 compensation point (Pa) + ! leaf boundary layer conductance of h20 + ! constrained vapor pressure + + call GetCanopyGasParameters(bc_in(s)%forc_pbot, & ! in + bc_in(s)%oair_pa(ifp), & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%tgcm_pa(ifp), & ! in + bc_in(s)%eair_pa(ifp), & ! in + bc_in(s)%esat_tv_pa(ifp), & ! in + bc_in(s)%rb_pa(ifp), & ! in + mm_kco2, & ! out + mm_ko2, & ! out + co2_cpoint, & ! out + cf, & ! out + gb_mol, & ! out + ceair) ! out + + ! ------------------------------------------------------------------------ + ! Part VI: Loop over all leaf layers. + ! The concept of leaf layers is a result of the radiative transfer scheme. + ! A leaf layer has uniform radiation environment. Leaf layers are a group + ! of vegetation surfaces (stems and leaves) which inhabit the same + ! canopy-layer "CL", have the same functional type "ft" and within those + ! two partitions are further partitioned into vertical layers where + ! downwelling radiation attenuates in order. + ! In this phase we loop over the leaf layers and calculate the + ! photosynthesis and respiration of the layer (since all biophysical + ! properties are homogeneous). After this step, we can loop through + ! our cohort list, associate each cohort with its list of leaf-layers + ! and transfer these quantities to the cohort. + ! With plant hydraulics, we must realize that photosynthesis and + ! respiration will be different for leaves of each cohort in the leaf + ! layers, as they will have there own hydraulic limitations. + ! NOTE: Only need to flush mask on the number of used pfts, not the whole + ! scratch space. + ! ------------------------------------------------------------------------ + rate_mask_z(:,1:numpft,:) = .false. + + if_any_cohorts: if(currentPatch%countcohorts > 0.0)then + currentCohort => currentPatch%tallest + do_cohort_drive: do while (associated(currentCohort)) ! Cohort loop + + ! Identify the canopy layer (cl), functional type (ft) + ! and the leaf layer (IV) for this cohort + ft = currentCohort%pft + cl = currentCohort%canopy_layer + + ! Calculate the cohort specific elai profile + ! And the top and bottom edges of the veg area index + ! of each layer bin are. Note, if the layers + ! sink below the ground snow line, then the effective + ! LAI and SAI start to shrink to zero, as well as + ! the difference between vaitop and vaibot. + if(currentCohort%treesai>0._r8)then + do iv = 1,currentCohort%nv + call VegAreaLayer(currentCohort%treelai, & + currentCohort%treesai, & + currentCohort%height, & + iv, & + currentCohort%nv, & + currentCohort%pft, & + sites(s)%snow_depth, & + cohort_vaitop(iv), & + cohort_vaibot(iv), & + cohort_layer_elai(iv), & + cohort_layer_esai(iv)) + end do + + cohort_elai = sum(cohort_layer_elai(1:currentCohort%nv)) + cohort_esai = sum(cohort_layer_esai(1:currentCohort%nv)) - else - cohort_layer_elai(:) = 0._r8 - cohort_layer_esai(:) = 0._r8 - cohort_vaitop(:) = 0._r8 - cohort_vaibot(:) = 0._r8 - cohort_elai = 0._r8 - cohort_esai = 0._r8 - end if - ! MLO. Assuming target to be related to leaf biomass when leaves are fully - ! flushed. But unsure whether this call is correct or not, shouldn't we get - ! the target value directly from the bstore_allom? - call bleaf(currentCohort%dbh,currentCohort%pft,& - currentCohort%crowndamage,currentCohort%canopy_trim,1.0_r8,store_c_target) - ! call bstore_allom(currentCohort%dbh,currentCohort%pft, & - ! currentCohort%canopy_trim,store_c_target) - - call storage_fraction_of_target(store_c_target, & - currentCohort%prt%GetState(store_organ, carbon12_element), & - frac) - call lowstorage_maintresp_reduction(frac,currentCohort%pft, & - maintresp_reduction_factor) - - ! are there any leaves of this pft in this layer? - canopy_mask_if: if(currentPatch%canopy_mask(cl,ft) == 1)then - - ! Loop over leaf-layers - leaf_layer_loop : do iv = 1,currentCohort%nv - - ! ------------------------------------------------------------ - ! If we are doing plant hydro-dynamics (or any run-type - ! where cohorts may generate different photosynthetic rates - ! of other cohorts in the same canopy-pft-layer combo), - ! we re-calculate the leaf biophysical rates for the - ! cohort-layer combo of interest. - ! but in the vanilla case, we only re-calculate if it has - ! not been done yet. - ! Other cases where we need to solve for every cohort - ! in every leaf layer: nutrient dynamic mode, multiple leaf - ! age classes - ! ------------------------------------------------------------ - - rate_mask_if: if ( .not.rate_mask_z(iv,ft,cl) .or. & - (hlm_use_planthydro.eq.itrue) .or. & - (radiation_model .eq. twostr_solver ) .or. & - (nleafage > 1) .or. & - (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then - - if (hlm_use_planthydro.eq.itrue ) then - - stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran ) - btran_eff = currentCohort%co_hydr%btran - - ! dinc_vai(:) is the total vegetation area index of each "leaf" layer - ! we convert to the leaf only portion of the increment - ! ------------------------------------------------------ - leaf_inc = dinc_vai(iv) * & - currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) - - ! Now calculate the cumulative top-down lai of the current layer's midpoint - lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1)) - - lai_layers_above = (dlower_vai(iv) - dinc_vai(iv)) * & - currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) - lai_current = min(leaf_inc, currentCohort%treelai - lai_layers_above) - cumulative_lai = lai_canopy_above + lai_layers_above + 0.5*lai_current - - leaf_psi = currentCohort%co_hydr%psi_ag(1) + else + cohort_layer_elai(:) = 0._r8 + cohort_layer_esai(:) = 0._r8 + cohort_vaitop(:) = 0._r8 + cohort_vaibot(:) = 0._r8 + cohort_elai = 0._r8 + cohort_esai = 0._r8 + end if + + ! MLO. Assuming target to be related to leaf biomass when leaves are fully + ! flushed. But unsure whether this call is correct or not, shouldn't we get + ! the target value directly from the bstore_allom? + call bleaf(currentCohort%dbh,currentCohort%pft,& + currentCohort%crowndamage,currentCohort%canopy_trim,1.0_r8,store_c_target) + ! call bstore_allom(currentCohort%dbh,currentCohort%pft, & + ! currentCohort%canopy_trim,store_c_target) - else + call storage_fraction_of_target(store_c_target, & + currentCohort%prt%GetState(store_organ, carbon12_element), & + frac) + call lowstorage_maintresp_reduction(frac,currentCohort%pft, & + maintresp_reduction_factor) + + ! are there any leaves of this pft in this layer? + canopy_mask_if: if(currentPatch%canopy_mask(cl,ft) == 1)then + + ! Loop over leaf-layers + leaf_layer_loop : do iv = 1,currentCohort%nv + + ! ------------------------------------------------------------ + ! If we are doing plant hydro-dynamics (or any run-type + ! where cohorts may generate different photosynthetic rates + ! of other cohorts in the same canopy-pft-layer combo), + ! we re-calculate the leaf biophysical rates for the + ! cohort-layer combo of interest. + ! but in the vanilla case, we only re-calculate if it has + ! not been done yet. + ! Other cases where we need to solve for every cohort + ! in every leaf layer: nutrient dynamic mode, multiple leaf + ! age classes + ! ------------------------------------------------------------ + + rate_mask_if: if ( .not.rate_mask_z(iv,ft,cl) .or. & + (hlm_use_planthydro.eq.itrue) .or. & + (radiation_model .eq. twostr_solver ) .or. & + (nleafage > 1) .or. & + (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then + + if (hlm_use_planthydro.eq.itrue ) then + + stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran ) + btran_eff = currentCohort%co_hydr%btran + + ! dinc_vai(:) is the total vegetation area index of each "leaf" layer + ! we convert to the leaf only portion of the increment + ! ------------------------------------------------------ + leaf_inc = dinc_vai(iv) * & + currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) + + ! Now calculate the cumulative top-down lai of the current layer's midpoint + lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + + lai_layers_above = (dlower_vai(iv) - dinc_vai(iv)) * & + currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) + lai_current = min(leaf_inc, currentCohort%treelai - lai_layers_above) + cumulative_lai = lai_canopy_above + lai_layers_above + 0.5*lai_current + + leaf_psi = currentCohort%co_hydr%psi_ag(1) + + else + + stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) ) + + btran_eff = currentPatch%btran_ft(ft) + ! For consistency sake, we use total LAI here, and not exposed + ! if the plant is under-snow, it will be effectively dormant for + ! the purposes of nscaler + + cumulative_lai = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + & + sum(currentPatch%tlai_profile(cl,ft,1:iv-1)) + & + 0.5*currentPatch%tlai_profile(cl,ft,iv) + + leaf_psi = fates_unset_r8 + + end if + + if(do_fates_salinity)then + btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) + endif + + + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used + ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al + ! (2010) Biogeosciences, 7, 1833-1859 + + kn = decay_coeff_vcmax(currentCohort%vcmax25top, & + prt_params%leafn_vert_scaler_coeff1(ft), & + prt_params%leafn_vert_scaler_coeff2(ft)) + + ! Scale for leaf nitrogen profile + nscaler = exp(-kn * cumulative_lai) + + ! Leaf maintenance respiration to match the base rate used in CN + ! but with the new temperature functions for C3 and C4 plants. - stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) ) + ! CN respiration has units: g C / g N [leaf] / s. This needs to be + ! converted from g C / g N [leaf] / s to umol CO2 / m**2 [leaf] / s - btran_eff = currentPatch%btran_ft(ft) - ! For consistency sake, we use total LAI here, and not exposed - ! if the plant is under-snow, it will be effectively dormant for - ! the purposes of nscaler + ! Then scale this value at the top of the canopy for canopy depth + ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) + select case(hlm_parteh_mode) + case (prt_carbon_allom_hyp) - cumulative_lai = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + & - sum(currentPatch%tlai_profile(cl,ft,1:iv-1)) + & - 0.5*currentPatch%tlai_profile(cl,ft,iv) + lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - leaf_psi = fates_unset_r8 + case (prt_cnp_flex_allom_hyp) + leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) + if( (leaf_c*slatop(ft)) > nearzero) then + leaf_n = currentCohort%prt%GetState(leaf_organ, nitrogen_element) + lnc_top = leaf_n / (slatop(ft) * leaf_c ) + else + lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) end if - if(do_fates_salinity)then - btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) - endif + ! If one wants to break coupling with dynamic N conentrations, + ! use the stoichiometry parameter + ! lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) + end select - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used - ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al - ! (2010) Biogeosciences, 7, 1833-1859 + ! Part VII: Calculate dark respiration (leaf maintenance) for this layer - kn = decay_coeff_vcmax(currentCohort%vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(ft), & - prt_params%leafn_vert_scaler_coeff2(ft)) + select case (maintresp_leaf_model) - ! Scale for leaf nitrogen profile - nscaler = exp(-kn * cumulative_lai) + case (lmrmodel_ryan_1991) - ! Leaf maintenance respiration to match the base rate used in CN - ! but with the new temperature functions for C3 and C4 plants. + call LeafLayerMaintenanceRespiration_Ryan_1991( lnc_top, & ! in + nscaler, & ! in + ft, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + lmr_z(iv,ft,cl)) ! out - ! CN respiration has units: g C / g N [leaf] / s. This needs to be - ! converted from g C / g N [leaf] / s to umol CO2 / m**2 [leaf] / s + case (lmrmodel_atkin_etal_2017) - ! Then scale this value at the top of the canopy for canopy depth - ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) - select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + call LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & ! in + cumulative_lai, & ! in + currentCohort%vcmax25top, & ! in + ft, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + currentPatch%tveg_lpa%GetMean(), & ! in + lmr_z(iv,ft,cl)) ! out + + case default + + write (fates_log(),*)'error, incorrect leaf respiration model specified' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end select + + ! Pre-process PAR absorbed per unit leaf area for different schemes + ! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves] + ! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves] + ! fsun = [m2 of sunlit leaves / m2 of total leaves] + ! laisun: m2 of exposed leaf, per m2 of crown. If this is the lowest layer + ! for the pft/canopy group, than the m2 per crown is probably not + ! as large as the layer above. + ! ------------------------------------------------------------------ + + if_radsolver: if(radiation_model.eq.norman_solver) then + + laisun = currentPatch%ed_laisun_z(cl,ft,iv) + laisha = currentPatch%ed_laisha_z(cl,ft,iv) + par_per_sunla = currentPatch%ed_parsun_z(cl,ft,iv) + par_per_shala = currentPatch%ed_parsha_z(cl,ft,iv) + canopy_area = currentPatch%canopy_area_profile(cl,ft,iv) + fsun = currentPatch%f_sun(cl,ft,iv) + + else ! Two-stream + + if(cohort_layer_elai(iv) > nearzero .and. currentPatch%solar_zenith_flag) then + + call FatesGetCohortAbsRad(currentPatch, currentCohort, ipar, & + cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, & + rb_abs, rd_abs, rb_abs_leaf, rd_abs_leaf, fsun) + + ! rd_abs_leaf: Watts of diffuse light absorbed by leaves over this + ! depth interval and ground footprint (m2) + ! rd_abs_leaf*fsun Watts of diffuse light absorbed by sunlit leaves + ! over this depth interval and ground footprint (m2) + ! rb_abs_leaf Watts of beam absorbed by sunlit leaves over this + ! depth interval and ground footprint (m2) + ! cohort_layer_elai*fsun Leaf area in sunlight within this interval and ground footprint + ! cohort_layer_elai*(1-fsun) Leaf area in shade within this interval and ground footprint + + laisun = (fsun*cohort_layer_elai(iv)) + laisha = ((1._r8 - fsun)*cohort_layer_elai(iv)) + if(fsun>nearzero) then + par_per_sunla = (rd_abs_leaf*fsun + rb_abs_leaf)! / laisun + else + par_per_sunla = 0._r8 + end if + par_per_shala = rd_abs_leaf*(1._r8-fsun) !/ laisha + canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) - lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) + else - case (prt_cnp_flex_allom_hyp) + par_per_sunla = 0._r8 + par_per_shala = 0._r8 + laisun = 0.5_r8*cohort_layer_elai(iv) + laisha = 0.5_r8*cohort_layer_elai(iv) + canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) + fsun = 0.5_r8 !avoid div0, should have no impact - leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) - if( (leaf_c*slatop(ft)) > nearzero) then - leaf_n = currentCohort%prt%GetState(leaf_organ, nitrogen_element) - lnc_top = leaf_n / (slatop(ft) * leaf_c ) - else - lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - end if + end if - ! If one wants to break coupling with dynamic N conentrations, - ! use the stoichiometry parameter - ! lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - - end select - - ! Part VII: Calculate dark respiration (leaf maintenance) for this layer - - select case (maintresp_leaf_model) - - case (lmrmodel_ryan_1991) - - call LeafLayerMaintenanceRespiration_Ryan_1991( lnc_top, & ! in - nscaler, & ! in - ft, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - lmr_z(iv,ft,cl)) ! out - - case (lmrmodel_atkin_etal_2017) - - call LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & ! in - cumulative_lai, & ! in - currentCohort%vcmax25top, & ! in - ft, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - currentPatch%tveg_lpa%GetMean(), & ! in - lmr_z(iv,ft,cl)) ! out - - case default - - write (fates_log(),*)'error, incorrect leaf respiration model specified' - call endrun(msg=errMsg(sourcefile, __LINE__)) - - end select - - ! Pre-process PAR absorbed per unit leaf area for different schemes - ! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves] - ! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves] - ! fsun = [m2 of sunlit leaves / m2 of total leaves] - ! laisun: m2 of exposed leaf, per m2 of crown. If this is the lowest layer - ! for the pft/canopy group, than the m2 per crown is probably not - ! as large as the layer above. - ! ------------------------------------------------------------------ - - if_radsolver: if(radiation_model.eq.norman_solver) then - - laisun = currentPatch%ed_laisun_z(cl,ft,iv) - laisha = currentPatch%ed_laisha_z(cl,ft,iv) - par_per_sunla = currentPatch%ed_parsun_z(cl,ft,iv) - par_per_shala = currentPatch%ed_parsha_z(cl,ft,iv) - canopy_area = currentPatch%canopy_area_profile(cl,ft,iv) - fsun = currentPatch%f_sun(cl,ft,iv) - - else ! Two-stream - - if(cohort_layer_elai(iv) > nearzero .and. currentPatch%solar_zenith_flag) then - - call FatesGetCohortAbsRad(currentPatch, currentCohort, ipar, & - cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, & - rb_abs, rd_abs, rb_abs_leaf, rd_abs_leaf, fsun) - - ! rd_abs_leaf: Watts of diffuse light absorbed by leaves over this - ! depth interval and ground footprint (m2) - ! rd_abs_leaf*fsun Watts of diffuse light absorbed by sunlit leaves - ! over this depth interval and ground footprint (m2) - ! rb_abs_leaf Watts of beam absorbed by sunlit leaves over this - ! depth interval and ground footprint (m2) - ! cohort_layer_elai*fsun Leaf area in sunlight within this interval and ground footprint - ! cohort_layer_elai*(1-fsun) Leaf area in shade within this interval and ground footprint - - laisun = (fsun*cohort_layer_elai(iv)) - laisha = ((1._r8 - fsun)*cohort_layer_elai(iv)) - if(fsun>nearzero) then - par_per_sunla = (rd_abs_leaf*fsun + rb_abs_leaf)! / laisun - else - par_per_sunla = 0._r8 - end if - par_per_shala = rd_abs_leaf*(1._r8-fsun) !/ laisha - canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) - - else + end if if_radsolver + + ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), + ! (2) maximum electron transport rate, (3) triose phosphate + ! utilization rate and (4) the initial slope of CO2 response curve + ! (C4 plants). Earlier we calculated their base rates as dictated + ! by their plant functional type and some simple scaling rules for + ! nitrogen limitation baesd on canopy position (not prognostic). + ! These rates are the specific rates used in the actual photosynthesis + ! calculations that take localized environmental effects (temperature) + ! into consideration. + + call LeafLayerBiophysicalRates(par_per_sunla, & ! in + ft, & ! in + currentCohort%vcmax25top, & ! in + currentCohort%jmax25top, & ! in + currentCohort%kp25top, & ! in + nscaler, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%dayl_factor_pa(ifp), & ! in + currentPatch%tveg_lpa%GetMean(), & ! in + currentPatch%tveg_longterm%GetMean(),& ! in + btran_eff, & ! in + vcmax_z, & ! out + jmax_z, & ! out + kp_z ) ! out + + ! Part IX: This call calculates the actual photosynthesis for the + ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. + + call LeafLayerPhotosynthesis(fsun, & ! in + par_per_sunla, & ! in + par_per_shala, & ! in + laisun, & ! in + laisha, & ! in + canopy_area, & ! in + ft, & ! in + vcmax_z, & ! in + jmax_z, & ! in + kp_z, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%esat_tv_pa(ifp), & ! in + bc_in(s)%forc_pbot, & ! in + bc_in(s)%cair_pa(ifp), & ! in + bc_in(s)%oair_pa(ifp), & ! in + btran_eff, & ! in + stomatal_intercept_btran, & ! in + cf, & ! in + gb_mol, & ! in + ceair, & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr_z(iv,ft,cl), & ! in + leaf_psi, & ! in + bc_in(s)%rb_pa(ifp), & ! in + currentPatch%psn_z(cl,ft,iv), & ! out + rs_z(iv,ft,cl), & ! out + anet_av_z(iv,ft,cl), & ! out + c13disc_z(cl,ft,iv)) ! out + + rate_mask_z(iv,ft,cl) = .true. + + end if rate_mask_if + end do leaf_layer_loop + + ! Zero cohort flux accumulators. + currentCohort%npp_tstep = 0.0_r8 + currentCohort%resp_tstep = 0.0_r8 + currentCohort%gpp_tstep = 0.0_r8 + currentCohort%rdark = 0.0_r8 + currentCohort%resp_m = 0.0_r8 + currentCohort%ts_net_uptake = 0.0_r8 + currentCohort%c13disc_clm = 0.0_r8 + + ! --------------------------------------------------------------- + ! Part VII: Transfer leaf flux rates (like maintenance respiration, + ! carbon assimilation and conductance) that are defined by the + ! leaf layer (which is area independent, ie /m2) onto each cohort + ! (where the rates become per cohort, ie /individual). Most likely + ! a sum over layers. + ! --------------------------------------------------------------- + nv = currentCohort%nv + + ! Temporary bypass to preserve B4B behavior + if(radiation_model.eq.norman_solver) then + + call ScaleLeafLayerFluxToCohort(nv, & !in + currentPatch%psn_z(cl,ft,1:nv), & !in + lmr_z(1:nv,ft,cl), & !in + rs_z(1:nv,ft,cl), & !in + currentPatch%elai_profile(cl,ft,1:nv), & !in + c13disc_z(cl, ft, 1:nv), & !in + currentCohort%c_area, & !in + currentCohort%n, & !in + bc_in(s)%rb_pa(ifp), & !in + maintresp_reduction_factor, & !in + currentCohort%g_sb_laweight, & !out + currentCohort%gpp_tstep, & !out + currentCohort%rdark, & !out + currentCohort%c13disc_clm, & !out + cohort_eleaf_area) !out - par_per_sunla = 0._r8 - par_per_shala = 0._r8 - laisun = 0.5_r8*cohort_layer_elai(iv) - laisha = 0.5_r8*cohort_layer_elai(iv) - canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) - fsun = 0.5_r8 !avoid div0, should have no impact - - end if + else - end if if_radsolver + call ScaleLeafLayerFluxToCohort(nv, & !in + currentPatch%psn_z(cl,ft,1:nv), & !in + lmr_z(1:nv,ft,cl), & !in + rs_z(1:nv,ft,cl), & !in + cohort_layer_elai(1:nv), & !in + c13disc_z(cl, ft, 1:nv), & !in + currentCohort%c_area, & !in + currentCohort%n, & !in + bc_in(s)%rb_pa(ifp), & !in + maintresp_reduction_factor, & !in + currentCohort%g_sb_laweight, & !out + currentCohort%gpp_tstep, & !out + currentCohort%rdark, & !out + currentCohort%c13disc_clm, & !out + cohort_eleaf_area) !out + end if - ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), - ! (2) maximum electron transport rate, (3) triose phosphate - ! utilization rate and (4) the initial slope of CO2 response curve - ! (C4 plants). Earlier we calculated their base rates as dictated - ! by their plant functional type and some simple scaling rules for - ! nitrogen limitation baesd on canopy position (not prognostic). - ! These rates are the specific rates used in the actual photosynthesis - ! calculations that take localized environmental effects (temperature) - ! into consideration. - call LeafLayerBiophysicalRates(par_per_sunla, & ! in - ft, & ! in - currentCohort%vcmax25top, & ! in - currentCohort%jmax25top, & ! in - currentCohort%kp25top, & ! in - nscaler, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%dayl_factor_pa(ifp), & ! in - currentPatch%tveg_lpa%GetMean(), & ! in - currentPatch%tveg_longterm%GetMean(),& ! in - btran_eff, & ! in - vcmax_z, & ! out - jmax_z, & ! out - kp_z ) ! out - - ! Part IX: This call calculates the actual photosynthesis for the - ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. - - call LeafLayerPhotosynthesis(fsun, & ! in - par_per_sunla, & ! in - par_per_shala, & ! in - laisun, & ! in - laisha, & ! in - canopy_area, & ! in - ft, & ! in - vcmax_z, & ! in - jmax_z, & ! in - kp_z, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - bc_in(s)%forc_pbot, & ! in - bc_in(s)%cair_pa(ifp), & ! in - bc_in(s)%oair_pa(ifp), & ! in - btran_eff, & ! in - stomatal_intercept_btran, & ! in - cf, & ! in - gb_mol, & ! in - ceair, & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr_z(iv,ft,cl), & ! in - leaf_psi, & ! in - bc_in(s)%rb_pa(ifp), & ! in - currentPatch%psn_z(cl,ft,iv), & ! out - rs_z(iv,ft,cl), & ! out - anet_av_z(iv,ft,cl), & ! out - c13disc_z(cl,ft,iv)) ! out - - rate_mask_z(iv,ft,cl) = .true. - - end if rate_mask_if - end do leaf_layer_loop - - ! Zero cohort flux accumulators. - currentCohort%npp_tstep = 0.0_r8 - currentCohort%resp_tstep = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%resp_m = 0.0_r8 - currentCohort%ts_net_uptake = 0.0_r8 - currentCohort%c13disc_clm = 0.0_r8 - - ! --------------------------------------------------------------- - ! Part VII: Transfer leaf flux rates (like maintenance respiration, - ! carbon assimilation and conductance) that are defined by the - ! leaf layer (which is area independent, ie /m2) onto each cohort - ! (where the rates become per cohort, ie /individual). Most likely - ! a sum over layers. - ! --------------------------------------------------------------- - nv = currentCohort%nv - - ! Temporary bypass to preserve B4B behavior - if(radiation_model.eq.norman_solver) then - - call ScaleLeafLayerFluxToCohort(nv, & !in - currentPatch%psn_z(cl,ft,1:nv), & !in - lmr_z(1:nv,ft,cl), & !in - rs_z(1:nv,ft,cl), & !in - currentPatch%elai_profile(cl,ft,1:nv), & !in - c13disc_z(cl, ft, 1:nv), & !in - currentCohort%c_area, & !in - currentCohort%n, & !in - bc_in(s)%rb_pa(ifp), & !in - maintresp_reduction_factor, & !in - currentCohort%g_sb_laweight, & !out - currentCohort%gpp_tstep, & !out - currentCohort%rdark, & !out - currentCohort%c13disc_clm, & !out - cohort_eleaf_area) !out - - else - - call ScaleLeafLayerFluxToCohort(nv, & !in - currentPatch%psn_z(cl,ft,1:nv), & !in - lmr_z(1:nv,ft,cl), & !in - rs_z(1:nv,ft,cl), & !in - cohort_layer_elai(1:nv), & !in - c13disc_z(cl, ft, 1:nv), & !in - currentCohort%c_area, & !in - currentCohort%n, & !in - bc_in(s)%rb_pa(ifp), & !in - maintresp_reduction_factor, & !in - currentCohort%g_sb_laweight, & !out - currentCohort%gpp_tstep, & !out - currentCohort%rdark, & !out - currentCohort%c13disc_clm, & !out - cohort_eleaf_area) !out - end if - - - ! Net Uptake does not need to be scaled, just transfer directly - currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC + ! Net Uptake does not need to be scaled, just transfer directly + currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC - else + else - ! In this case, the cohort had no leaves, - ! so no productivity,conductance, transpiration uptake - ! or dark respiration - cohort_eleaf_area = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%g_sb_laweight = 0.0_r8 - currentCohort%ts_net_uptake(:) = 0.0_r8 + ! In this case, the cohort had no leaves, + ! so no productivity,conductance, transpiration uptake + ! or dark respiration + cohort_eleaf_area = 0.0_r8 + currentCohort%gpp_tstep = 0.0_r8 + currentCohort%rdark = 0.0_r8 + currentCohort%g_sb_laweight = 0.0_r8 + currentCohort%ts_net_uptake(:) = 0.0_r8 - end if canopy_mask_if + end if canopy_mask_if - ! ------------------------------------------------------------------ - ! Part VIII: Calculate maintenance respiration in the sapwood and - ! fine root pools. - ! ------------------------------------------------------------------ + ! ------------------------------------------------------------------ + ! Part VIII: Calculate maintenance respiration in the sapwood and + ! fine root pools. + ! ------------------------------------------------------------------ - ! Calculate the amount of nitrogen in the above and below ground - ! stem and root pools, used for maint resp - ! We are using the fine-root C:N ratio as an approximation for - ! the sapwood pools. - ! Units are in (kgN/plant) - ! ------------------------------------------------------------------ + ! Calculate the amount of nitrogen in the above and below ground + ! stem and root pools, used for maint resp + ! We are using the fine-root C:N ratio as an approximation for + ! the sapwood pools. + ! Units are in (kgN/plant) + ! ------------------------------------------------------------------ - sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, carbon12_element) + sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) + fnrt_c = currentCohort%prt%GetState(fnrt_organ, carbon12_element) - if (hlm_use_tree_damage .eq. itrue) then + if (hlm_use_tree_damage .eq. itrue) then - ! Crown damage currenly only reduces the aboveground portion of - ! sapwood. Therefore we calculate the aboveground and the belowground portion - ! sapwood for use in stem respiration. - call GetCrownReduction(currentCohort%crowndamage, crown_reduction) + ! Crown damage currenly only reduces the aboveground portion of + ! sapwood. Therefore we calculate the aboveground and the belowground portion + ! sapwood for use in stem respiration. + call GetCrownReduction(currentCohort%crowndamage, crown_reduction) - else - crown_reduction = 0.0_r8 - end if + else + crown_reduction = 0.0_r8 + end if - ! If crown reduction is zero, undamaged sapwood target will equal sapwood carbon - agb_frac = prt_params%allom_agb_frac(currentCohort%pft) - branch_frac = param_derived%branch_frac(currentCohort%pft) - sapw_c_undamaged = sapw_c / (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) + ! If crown reduction is zero, undamaged sapwood target will equal sapwood carbon + agb_frac = prt_params%allom_agb_frac(currentCohort%pft) + branch_frac = param_derived%branch_frac(currentCohort%pft) + sapw_c_undamaged = sapw_c / (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) - ! Undamaged below ground portion - sapw_c_bgw = sapw_c_undamaged * (1.0_r8 - agb_frac) + ! Undamaged below ground portion + sapw_c_bgw = sapw_c_undamaged * (1.0_r8 - agb_frac) - ! Damaged aboveground portion - sapw_c_agw = sapw_c - sapw_c_bgw + ! Damaged aboveground portion + sapw_c_agw = sapw_c - sapw_c_bgw - select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + select case(hlm_parteh_mode) + case (prt_carbon_allom_hyp) - live_stem_n = sapw_c_agw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + live_stem_n = sapw_c_agw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - live_croot_n = sapw_c_bgw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + live_croot_n = sapw_c_bgw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) + fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) - case(prt_cnp_flex_allom_hyp) + case(prt_cnp_flex_allom_hyp) - live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & - currentCohort%prt%GetState(sapw_organ, nitrogen_element) + live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & + currentCohort%prt%GetState(sapw_organ, nitrogen_element) - live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & - currentCohort%prt%GetState(sapw_organ, nitrogen_element) + live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & + currentCohort%prt%GetState(sapw_organ, nitrogen_element) - fnrt_n = currentCohort%prt%GetState(fnrt_organ, nitrogen_element) + fnrt_n = currentCohort%prt%GetState(fnrt_organ, nitrogen_element) - if (hlm_use_tree_damage .eq. itrue) then + if (hlm_use_tree_damage .eq. itrue) then - sapw_n = currentCohort%prt%GetState(sapw_organ, nitrogen_element) + sapw_n = currentCohort%prt%GetState(sapw_organ, nitrogen_element) - sapw_n_undamaged = sapw_n / & - (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) + sapw_n_undamaged = sapw_n / & + (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) - sapw_n_bgw = sapw_n_undamaged * (1.0_r8 - agb_frac) - sapw_n_agw = sapw_n - sapw_n_bgw + sapw_n_bgw = sapw_n_undamaged * (1.0_r8 - agb_frac) + sapw_n_agw = sapw_n - sapw_n_bgw - live_croot_n = sapw_n_bgw + live_croot_n = sapw_n_bgw - live_stem_n = sapw_n_agw + live_stem_n = sapw_n_agw - end if + end if - ! If one wants to break coupling with dynamic N conentrations, - ! use the stoichiometry parameter - ! - ! live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & - ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - ! live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & - ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - ! fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) + ! If one wants to break coupling with dynamic N conentrations, + ! use the stoichiometry parameter + ! + ! live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & + ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + ! live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & + ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + ! fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) - case default + case default - end select + end select - !------------------------------------------------------------------------------ - ! Calculate Whole Plant Respiration - ! (this doesn't really need to be in this iteration at all, surely?) - ! Response: (RGK 12-2016): I think the positioning of these calls is - ! appropriate as of now. Maintenance calculations in sapwood and roots - ! vary by cohort and with changing temperature at the minimum, and there are - ! no sub-pools chopping up those pools any finer that need to be dealt with. - !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------ + ! Calculate Whole Plant Respiration + ! (this doesn't really need to be in this iteration at all, surely?) + ! Response: (RGK 12-2016): I think the positioning of these calls is + ! appropriate as of now. Maintenance calculations in sapwood and roots + ! vary by cohort and with changing temperature at the minimum, and there are + ! no sub-pools chopping up those pools any finer that need to be dealt with. + !------------------------------------------------------------------------------ - ! Live stem MR (kgC/plant/s) (above ground sapwood) - ! ------------------------------------------------------------------ - if ( int(woody(ft)) == itrue) then - tcwood = q10_mr**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8) - ! kgC/s = kgN * kgC/kgN/s - currentCohort%livestem_mr = live_stem_n * maintresp_nonleaf_baserate * tcwood * maintresp_reduction_factor - else - currentCohort%livestem_mr = 0._r8 - end if + ! Live stem MR (kgC/plant/s) (above ground sapwood) + ! ------------------------------------------------------------------ + if ( int(woody(ft)) == itrue) then + tcwood = q10_mr**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8) + ! kgC/s = kgN * kgC/kgN/s + currentCohort%livestem_mr = live_stem_n * maintresp_nonleaf_baserate * tcwood * maintresp_reduction_factor + else + currentCohort%livestem_mr = 0._r8 + end if - ! Fine Root MR (kgC/plant/s) - ! and calculate the N fixation rate as a function of the fixation-specific root respiration - ! for now use dev_arbitrary_pft as scaling term between 0 and 1 as additional increment of root respiration used for N fixation - ! ------------------------------------------------------------------ - currentCohort%froot_mr = 0._r8 - currentCohort%sym_nfix_tstep = 0._r8 + ! Fine Root MR (kgC/plant/s) + ! and calculate the N fixation rate as a function of the fixation-specific root respiration + ! for now use dev_arbitrary_pft as scaling term between 0 and 1 as additional increment of root respiration used for N fixation + ! ------------------------------------------------------------------ + currentCohort%froot_mr = 0._r8 + currentCohort%sym_nfix_tstep = 0._r8 - ! n_fixation is integrated over the course of the day - ! this variable is zeroed at the end of the FATES dynamics sequence + ! n_fixation is integrated over the course of the day + ! this variable is zeroed at the end of the FATES dynamics sequence - do j = 1,bc_in(s)%nlevsoil - tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) + do j = 1,bc_in(s)%nlevsoil + tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) + + fnrt_mr_layer = fnrt_n * maintresp_nonleaf_baserate * tcsoi * rootfr_ft(ft,j) * maintresp_reduction_factor - fnrt_mr_layer = fnrt_n * maintresp_nonleaf_baserate * tcsoi * rootfr_ft(ft,j) * maintresp_reduction_factor + ! calculate the cost of carbon for N fixation in each soil layer and calculate N fixation rate based on that [kgC / kgN] - ! calculate the cost of carbon for N fixation in each soil layer and calculate N fixation rate based on that [kgC / kgN] + call RootLayerNFixation(bc_in(s)%t_soisno_sl(j),ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,nfix_layer) - call RootLayerNFixation(bc_in(s)%t_soisno_sl(j),ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,nfix_layer) + currentCohort%froot_mr = currentCohort%froot_mr + fnrt_mr_nfix_layer + fnrt_mr_layer - currentCohort%froot_mr = currentCohort%froot_mr + fnrt_mr_nfix_layer + fnrt_mr_layer + currentCohort%sym_nfix_tstep = currentCohort%sym_nfix_tstep + nfix_layer - currentCohort%sym_nfix_tstep = currentCohort%sym_nfix_tstep + nfix_layer + enddo + ! Coarse Root MR (kgC/plant/s) (below ground sapwood) + ! ------------------------------------------------------------------ + if ( int(woody(ft)) == itrue) then + currentCohort%livecroot_mr = 0._r8 + do j = 1,bc_in(s)%nlevsoil + ! Soil temperature used to adjust base rate of MR + tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) + currentCohort%livecroot_mr = currentCohort%livecroot_mr + & + live_croot_n * maintresp_nonleaf_baserate * tcsoi * & + rootfr_ft(ft,j) * maintresp_reduction_factor enddo + else + currentCohort%livecroot_mr = 0._r8 + end if - ! Coarse Root MR (kgC/plant/s) (below ground sapwood) - ! ------------------------------------------------------------------ - if ( int(woody(ft)) == itrue) then - currentCohort%livecroot_mr = 0._r8 - do j = 1,bc_in(s)%nlevsoil - ! Soil temperature used to adjust base rate of MR - tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) - currentCohort%livecroot_mr = currentCohort%livecroot_mr + & - live_croot_n * maintresp_nonleaf_baserate * tcsoi * & - rootfr_ft(ft,j) * maintresp_reduction_factor - enddo - else - currentCohort%livecroot_mr = 0._r8 - end if + ! ------------------------------------------------------------------ + ! Part IX: Perform some unit conversions (rate to integrated) and + ! calcualate some fluxes that are sums and nets of the base fluxes + ! ------------------------------------------------------------------ + + if ( debug ) write(fates_log(),*) 'EDPhoto 904 ', currentCohort%resp_m + if ( debug ) write(fates_log(),*) 'EDPhoto 905 ', currentCohort%rdark + if ( debug ) write(fates_log(),*) 'EDPhoto 906 ', currentCohort%livestem_mr + if ( debug ) write(fates_log(),*) 'EDPhoto 907 ', currentCohort%livecroot_mr + if ( debug ) write(fates_log(),*) 'EDPhoto 908 ', currentCohort%froot_mr - ! ------------------------------------------------------------------ - ! Part IX: Perform some unit conversions (rate to integrated) and - ! calcualate some fluxes that are sums and nets of the base fluxes - ! ------------------------------------------------------------------ - if ( debug ) write(fates_log(),*) 'EDPhoto 904 ', currentCohort%resp_m - if ( debug ) write(fates_log(),*) 'EDPhoto 905 ', currentCohort%rdark - if ( debug ) write(fates_log(),*) 'EDPhoto 906 ', currentCohort%livestem_mr - if ( debug ) write(fates_log(),*) 'EDPhoto 907 ', currentCohort%livecroot_mr - if ( debug ) write(fates_log(),*) 'EDPhoto 908 ', currentCohort%froot_mr + ! add on whole plant respiration values in kgC/indiv/s-1 + currentCohort%resp_m = currentCohort%livestem_mr + & + currentCohort%livecroot_mr + & + currentCohort%froot_mr + ! no drought response right now.. something like: + ! resp_m = resp_m * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft) * & + ! EDPftvarcon_inst%resp_drought_response(ft)) - ! add on whole plant respiration values in kgC/indiv/s-1 - currentCohort%resp_m = currentCohort%livestem_mr + & - currentCohort%livecroot_mr + & - currentCohort%froot_mr + currentCohort%resp_m = currentCohort%resp_m + currentCohort%rdark - ! no drought response right now.. something like: - ! resp_m = resp_m * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft) * & - ! EDPftvarcon_inst%resp_drought_response(ft)) + ! save as a diagnostic the un-throttled maintenance respiration to be able to know how strong this is + currentCohort%resp_m_unreduced = currentCohort%resp_m / maintresp_reduction_factor - currentCohort%resp_m = currentCohort%resp_m + currentCohort%rdark + ! convert from kgC/indiv/s to kgC/indiv/timestep + currentCohort%resp_m = currentCohort%resp_m * dtime + currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime + currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime - ! save as a diagnostic the un-throttled maintenance respiration to be able to know how strong this is - currentCohort%resp_m_unreduced = currentCohort%resp_m / maintresp_reduction_factor + if ( debug ) write(fates_log(),*) 'EDPhoto 911 ', currentCohort%gpp_tstep + if ( debug ) write(fates_log(),*) 'EDPhoto 912 ', currentCohort%resp_tstep + if ( debug ) write(fates_log(),*) 'EDPhoto 913 ', currentCohort%resp_m - ! convert from kgC/indiv/s to kgC/indiv/timestep - currentCohort%resp_m = currentCohort%resp_m * dtime - currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime - currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime - if ( debug ) write(fates_log(),*) 'EDPhoto 911 ', currentCohort%gpp_tstep - if ( debug ) write(fates_log(),*) 'EDPhoto 912 ', currentCohort%resp_tstep - if ( debug ) write(fates_log(),*) 'EDPhoto 913 ', currentCohort%resp_m + currentCohort%resp_g_tstep = prt_params%grperc(ft) * & + (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) - currentCohort%resp_g_tstep = prt_params%grperc(ft) * & - (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) + currentCohort%resp_tstep = currentCohort%resp_m + & + currentCohort%resp_g_tstep ! kgC/indiv/ts + currentCohort%npp_tstep = currentCohort%gpp_tstep - & + currentCohort%resp_tstep ! kgC/indiv/ts + ! Accumulate the combined conductance (stomatal+leaf boundary layer) + ! Note that currentCohort%g_sb_laweight is weighted by the leaf area + ! of each cohort and has units of [m/s] * [m2 leaf] - currentCohort%resp_tstep = currentCohort%resp_m + & - currentCohort%resp_g_tstep ! kgC/indiv/ts - currentCohort%npp_tstep = currentCohort%gpp_tstep - & - currentCohort%resp_tstep ! kgC/indiv/ts + g_sb_leaves = g_sb_leaves + currentCohort%g_sb_laweight - ! Accumulate the combined conductance (stomatal+leaf boundary layer) - ! Note that currentCohort%g_sb_laweight is weighted by the leaf area - ! of each cohort and has units of [m/s] * [m2 leaf] + ! Accumulate the total effective leaf area from all cohorts + ! in this patch. Normalize by canopy area outside the loop + patch_la = patch_la + cohort_eleaf_area - g_sb_leaves = g_sb_leaves + currentCohort%g_sb_laweight + currentCohort => currentCohort%shorter + enddo do_cohort_drive - ! Accumulate the total effective leaf area from all cohorts - ! in this patch. Normalize by canopy area outside the loop - patch_la = patch_la + cohort_eleaf_area + end if if_any_cohorts - currentCohort => currentCohort%shorter - enddo do_cohort_drive + ! Normalize canopy total conductance by the effective LAI + ! The value here was integrated over each cohort x leaf layer + ! and was weighted by m2 of effective leaf area for each layer + ! 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(preserve_b4b) then + patch_la = patch_la/ currentPatch%total_canopy_area + end if - end if if_any_cohorts + ! Normalize canopy total conductance by the effective LAI + ! The value here was integrated over each cohort x leaf layer + ! and was weighted by m2 of effective leaf area for each layer - ! Normalize canopy total conductance by the effective LAI - ! The value here was integrated over each cohort x leaf layer - ! and was weighted by m2 of effective leaf area for each layer + if_any_lai: if(patch_la>tiny(patch_la)) then + + ! Normalize the leaf-area weighted canopy conductance + ! The denominator is the total effective leaf area in the canopy, + ! units of [m/s]*[m2] / [m2] = [m/s] ! 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(preserve_b4b) then - patch_la = patch_la/ currentPatch%total_canopy_area - end if - - ! Normalize canopy total conductance by the effective LAI - ! The value here was integrated over each cohort x leaf layer - ! and was weighted by m2 of effective leaf area for each layer - - if_any_lai: if(patch_la>tiny(patch_la)) then - - ! Normalize the leaf-area weighted canopy conductance - ! The denominator is the total effective leaf area in the canopy, - ! units of [m/s]*[m2] / [m2] = [m/s] - ! 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_preserve_b4b3: if(preserve_b4b) then - elai = calc_areaindex(currentPatch,'elai') - g_sb_leaves = g_sb_leaves / (elai*currentPatch%total_canopy_area) - else - g_sb_leaves = g_sb_leaves / max(0.1_r8*currentPatch%total_canopy_area,patch_la) - end if if_preserve_b4b3 - - - if_above_mincond: if( g_sb_leaves > (1._r8/rsmax0) ) then - - ! Combined mean leaf resistance is the inverse of mean leaf conductance - r_sb_leaves = 1.0_r8/g_sb_leaves - - if (r_sb_leaves (1._r8/rsmax0) ) then + + ! Combined mean leaf resistance is the inverse of mean leaf conductance + r_sb_leaves = 1.0_r8/g_sb_leaves + + if (r_sb_leaves currentPatch%younger end do diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index f1e7de14ea..081eb095a2 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -175,8 +175,6 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch => currentSite%oldest_patch; do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - litt_c => currentPatch%litter(element_pos(carbon12_element)) ! How much live grass is there? @@ -329,11 +327,10 @@ subroutine charecteristics_of_fuel ( currentSite ) ! FIX(SPM,032414) refactor... if(write_SF == itrue.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'problem with spitfire fuel averaging' - endif - endif !nocomp_pft_label check + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'problem with spitfire fuel averaging' + endif + currentPatch => currentPatch%younger - enddo !end patch loop end subroutine charecteristics_of_fuel @@ -366,6 +363,7 @@ subroutine wind_effect ( 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 currentPatch => currentPatch%younger endif @@ -383,15 +381,12 @@ subroutine wind_effect ( currentSite, bc_in) ! average_wspeed = 0.0_r8 tree_fraction = 0.0_r8 grass_fraction = 0.0_r8 - currentPatch=>currentSite%oldest_patch; + currentPatch=>currentSite%oldest_patch do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - currentPatch%total_tree_area = 0.0_r8 total_grass_area = 0.0_r8 currentCohort => currentPatch%tallest - do while(associated(currentCohort)) if (debug) write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area if( prt_params%woody(currentCohort%pft) == itrue)then @@ -411,8 +406,6 @@ subroutine wind_effect ( currentSite, bc_in) write(fates_log(),*) 'SF AREA ',AREA endif - endif !nocomp_pft_label check - currentPatch => currentPatch%younger enddo !currentPatch loop @@ -424,16 +417,13 @@ subroutine wind_effect ( currentSite, bc_in) grass_fraction, tree_fraction, bare_fraction endif - currentPatch=>currentSite%oldest_patch; - - do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then + currentPatch=>currentSite%oldest_patch + do while(associated(currentPatch)) currentPatch%total_tree_area = min(currentPatch%total_tree_area,currentPatch%area) + ! effect_wspeed in units m/min currentPatch%effect_wspeed = currentSite%wind * (tree_fraction*0.4_r8+(grass_fraction+bare_fraction)*0.6_r8) - - endif ! nocomp_pft_label check currentPatch => currentPatch%younger enddo !end patch loop @@ -473,8 +463,6 @@ subroutine rate_of_spread ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - ! remove mineral content from net fuel load per Thonicke 2010 for ir calculation currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals @@ -579,9 +567,7 @@ 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 currentPatch => currentPatch%younger - enddo !end patch loop end subroutine rate_of_spread @@ -609,8 +595,6 @@ subroutine ground_fuel_consumption ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - currentPatch%burnt_frac_litter(:) = 1.0_r8 ! Calculate fraction of litter is burnt for all classes. ! Equation B1 in Thonicke et al. 2010--- @@ -670,9 +654,7 @@ 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; + currentPatch=>currentPatch%younger enddo !end patch loop end subroutine ground_fuel_consumption @@ -765,8 +747,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 @@ -871,7 +851,6 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) endif endif ! NF ignitions check - endif ! nocomp_pft_label check currentPatch => currentPatch%younger @@ -904,8 +883,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; @@ -938,7 +915,6 @@ subroutine crown_scorching ( currentSite ) end do endif !fire - endif !nocomp_pft_label currentPatch => currentPatch%younger; enddo !end patch loop @@ -962,7 +938,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 @@ -1004,7 +979,6 @@ subroutine crown_damage ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check currentPatch => currentPatch%younger; @@ -1031,8 +1005,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)) @@ -1057,9 +1029,8 @@ subroutine cambial_damage_kill ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check - currentPatch=>currentPatch%younger; + currentPatch=>currentPatch%younger enddo !end patch loop @@ -1084,8 +1055,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)) @@ -1105,7 +1074,6 @@ subroutine post_fire_mortality ( currentSite ) enddo !end cohort loop endif !fire? - endif !nocomp_pft_label check currentPatch => currentPatch%younger diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index 24cf38fbb7..c4a3b21c1a 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 cb642c1289..a747c71b35 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 @@ -100,33 +99,53 @@ subroutine FatesNormalizedCanopyRadiation(nsites, sites, bc_in, bc_out ) ifp = 0 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) + - ifp = ifp+1 + if_noveg: if(currentPatch%total_canopy_area < nearzero) then - ! 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 + ! For a bare-ground patch, either imposed by a no-comp mode + ! or because we are in a desert and no plants can grown, + ! there is no light is absorbed by the canopy + ! (which isn't there anyway), all light it transmitted. The + ! albedo will be the soil albedo. Note: On the host-side + ! these numbers will be weighted by a zero anyway, so + ! any initialization will probably work fine because + ! the host has its own bare-ground methods, but we set these + ! values for constistency sake. - 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) + 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) + end do + bc_out(s)%fabi_parb(ifp,:) = 0._r8 + bc_out(s)%fabd_parb(ifp,:) = 0._r8 + bc_out(s)%ftdd_parb(ifp,:) = 1._r8 + bc_out(s)%ftid_parb(ifp,:) = 1._r8 + bc_out(s)%ftii_parb(ifp,:) = 1._r8 + + + else !if_noveg + if(radiation_model.eq.twostr_solver) then call currentPatch%twostr%CanopyPrep(bc_in(s)%fcansno_pa(ifp)) @@ -240,7 +259,7 @@ subroutine FatesNormalizedCanopyRadiation(nsites, sites, bc_in, bc_out ) end if if_nrad endif if_zenith_flag - end if if_notbareground + end if if_noveg currentPatch => currentPatch%younger end do ! Loop linked-list patches @@ -287,218 +306,212 @@ 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 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, nlevleaf - if(area_vlpfcl(iv,ft,cl) 1._r8)then + write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & + sunlai,shalai + endif + end if - end associate + elai = calc_areaindex(cpatch,'elai') - end if if_zenithflag - endif if_norm_twostr - - end if if_notbareground - - cpatch => cpatch%younger - enddo + 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 + + 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, nlevleaf + if(area_vlpfcl(iv,ft,cl) cpatch%younger enddo - return - end subroutine FatesSunShadeFracs + + enddo + return + +end subroutine FatesSunShadeFracs end module FatesRadiationDriveMod From 19297fe32002be75a731418c2bcb127a4171dedc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 22 Nov 2024 17:33:51 -0500 Subject: [PATCH 02/11] Removing patch index 0 filters --- main/EDMainMod.F90 | 11 ++--------- main/FatesHistoryInterfaceMod.F90 | 10 +++------- main/FatesInterfaceMod.F90 | 16 +++++++--------- 3 files changed, 12 insertions(+), 25 deletions(-) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index cd98993a6b..c7dc7f8066 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -209,15 +209,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/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index c2014b75e2..f9a43f0583 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -5002,8 +5002,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 @@ -5011,11 +5009,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 eea65eaa0a..81a4fb324f 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -792,7 +792,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 @@ -2113,8 +2113,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ifp=0 site_npp = 0._r8 cpatch => 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 +2127,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 +2161,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 +2175,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 From 968f858e432b7906f3708f7b167fbb39ac66d198 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 2 Dec 2024 12:31:10 -0500 Subject: [PATCH 03/11] removed another nocomp_bareground filter --- biogeochem/FatesSoilBGCFluxMod.F90 | 194 ++++++++++++++--------------- 1 file changed, 96 insertions(+), 98 deletions(-) diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index 2c5b7d9b18..cb9fbe7fee 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 @@ -288,107 +287,106 @@ 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_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_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 From 9bb650f5116642bd32cf37c1cc33ba4bd2ded6b0 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 2 Dec 2024 13:56:58 -0500 Subject: [PATCH 04/11] zeroing some ch4 bc arrays so the bg-patch will be pre-zeroed --- main/FatesInterfaceMod.F90 | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 81a4fb324f..eacb67b915 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -654,20 +654,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 From fcec87346b4ee2827b484696a85b99c4f129cf37 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 2 Dec 2024 13:58:22 -0500 Subject: [PATCH 05/11] updates to maxpatch counting --- main/EDParamsMod.F90 | 2 +- main/FatesInterfaceMod.F90 | 17 ++++++++++++----- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index c3d8a76273..e1b3dc5d91 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -845,7 +845,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/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index eea65eaa0a..82d2e714bd 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -783,6 +783,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 @@ -812,12 +820,11 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade if(hlm_use_nocomp==itrue) then maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) - maxpatch_total = sum(maxpatches_by_landuse(:)) + maxpatch_total = sum(maxpatches_by_landuse(:)+1) - !if(maxpatch_primary Date: Mon, 2 Dec 2024 14:15:23 -0500 Subject: [PATCH 06/11] updated comment --- main/FatesInterfaceMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b4cb68b561..64841b1ec1 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -824,7 +824,7 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade maxpatch_total = maxpatches_by_landuse(1) end if - ! maxpatch_total does not include the bare ground (so add 1) + ! maxpatch_total does not include HLM-side bare ground (so add 1) fates_maxPatchesPerSite = maxpatch_total+1 end if From 3759083eaefd455ea51487d40bbe5ac71a47ae0d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 2 Dec 2024 18:28:48 -0500 Subject: [PATCH 07/11] added use_luh condition to maxpatch_total --- main/FatesInterfaceMod.F90 | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 64841b1ec1..780430ada7 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -813,18 +813,27 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,param_reade ! 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 - maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) - maxpatch_total = sum(maxpatches_by_landuse(:)+1) + !! REVIEWERS: NEED THIS NEXT LINE? + !maxpatches_by_landuse(secondaryland) = max(maxpatches_by_landuse(secondaryland),fates_numpft) + ! We hold the bareground area in one of the patches + maxpatch_total = sum(maxpatches_by_landuse(:))+1 else - ! Check-in with team on this... (RGK) - maxpatch_total = maxpatches_by_landuse(1) + 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. + fates_maxPatchesPerSite = maxpatch_total+1 end if From 1b128b9c38adca83a8c41b5fea114bd01f4f0977 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 3 Dec 2024 10:43:04 -0700 Subject: [PATCH 08/11] Updated check on patch counts for non luh scenarios --- main/EDPftvarcon.F90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) 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 From 6798251f08baba1ab5eae055a73804f96431cdba Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 3 Dec 2024 14:23:29 -0700 Subject: [PATCH 09/11] reverting numpotentialcanopy layer minimum to 1 --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 14283d9025..46ee85d045 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -2290,7 +2290,7 @@ function NumPotentialCanopyLayers(currentPatch,site_spread,include_substory) res real(r8) :: c_area real(r8) :: arealayer - z = 0 + z = 1 currentCohort => currentPatch%tallest do while (associated(currentCohort)) z = max(z,currentCohort%canopy_layer) From 5d892a1dbbc1ead2eadedabaa901202a13fe4aea Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 10 Dec 2024 11:37:37 -0500 Subject: [PATCH 10/11] Updated patch allocation initialization to include discussion points with ckoven, and to have a b4b preservation bypass. --- main/EDInitMod.F90 | 69 +++++++++++++++++++++----------------- main/FatesInterfaceMod.F90 | 60 +++++++++++++++++++++------------ 2 files changed, 77 insertions(+), 52 deletions(-) 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/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3650e4a179..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 @@ -810,32 +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_luh==ifalse) then - maxpatches_by_landuse(secondaryland:n_landuse_cats) = 0 - end if - - if(hlm_use_nocomp==itrue) then - maxpatches_by_landuse(primaryland) = max(maxpatches_by_landuse(primaryland),fates_numpft) - !! REVIEWERS: NEED THIS NEXT LINE? - !maxpatches_by_landuse(secondaryland) = max(maxpatches_by_landuse(secondaryland),fates_numpft) - - ! We hold the bareground area in one of the patches - maxpatch_total = sum(maxpatches_by_landuse(:))+1 + 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 - maxpatch_total = sum(maxpatches_by_landuse(:)) - end if + + ! 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. + + end if if_preserve_b4b - ! 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. - fates_maxPatchesPerSite = maxpatch_total+1 + end if end if From 4bdf4c2d1a1cc9ce75114271b5d25e9cf02611bd Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 10 Dec 2024 09:41:30 -0700 Subject: [PATCH 11/11] added patch pointer protections in spitfire to accomodate instances where only the bareground nocomp patch exists --- fire/SFMainMod.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index e2590224fc..90fd1ceabd 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -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 @@ -501,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