diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 46db5845f1..163ee835d6 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1601,19 +1601,19 @@ subroutine leaf_area_profile( currentSite ) currentCohort => currentPatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft - min_chite = currentCohort%hite - currentCohort%hite * EDPftvarcon_inst%crown(ft) + min_chite = currentCohort%hite - currentCohort%hite * prt_params%crown(ft) max_chite = currentCohort%hite do iv = 1,N_HITE_BINS frac_canopy(iv) = 0.0_r8 ! this layer is in the middle of the canopy if(max_chite > maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*EDPftvarcon_inst%crown(ft))) + frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*prt_params%crown(ft))) ! this is the layer with the bottom of the canopy in it. elseif(min_chite < maxh(iv).and.min_chite > minh(iv).and.max_chite > maxh(iv))then - frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) + frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*prt_params%crown(ft)) ! this is the layer with the top of the canopy in it. elseif(max_chite > minh(iv).and.max_chite < maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) + frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*prt_params%crown(ft)) elseif(max_chite < maxh(iv).and.min_chite > minh(iv))then !the whole cohort is within this layer. frac_canopy(iv) = 1.0_r8 endif @@ -1709,11 +1709,11 @@ subroutine leaf_area_profile( currentSite ) layer_top_hite = currentCohort%hite - & ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & - EDPftvarcon_inst%crown(currentCohort%pft) ) + prt_params%crown(currentCohort%pft) ) layer_bottom_hite = currentCohort%hite - & ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & - EDPftvarcon_inst%crown(currentCohort%pft) ) + prt_params%crown(currentCohort%pft) ) fraction_exposed = 1.0_r8 if(currentSite%snow_depth > layer_top_hite)then diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 2fa98aa59f..ecdb731621 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -324,7 +324,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & call InitHydrCohort(currentSite,new_cohort) ! This calculates node heights - call UpdatePlantHydrNodes(new_cohort%co_hydr,new_cohort%pft, & + call UpdatePlantHydrNodes(new_cohort,new_cohort%pft, & new_cohort%hite,currentSite%si_hydr) ! This calculates volumes and lengths diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index e47934715d..6b315d4ef8 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -1983,7 +1983,7 @@ end subroutine h2d_martcano ! ===================================================================================== - subroutine CrownDepth(height,crown_depth) + subroutine CrownDepth(height,ft,crown_depth) ! ----------------------------------------------------------------------------------- ! This routine returns the depth of a plant's crown. Which is the length @@ -1993,14 +1993,20 @@ subroutine CrownDepth(height,crown_depth) ! optioned. ! ----------------------------------------------------------------------------------- - real(r8),intent(in) :: height ! The height of the plant [m] + real(r8),intent(in) :: height ! The height of the plant [m] + integer,intent(in) :: ft ! functional type index real(r8),intent(out) :: crown_depth ! The depth of the crown [m] - + ! Alternative Hypothesis: ! crown depth from Poorter, Bongers & Bongers ! crown_depth = exp(-1.169_r8)*cCohort%hite**1.098_r8 - - crown_depth = min(height,0.1_r8) + + ! Alternative Hypothesis: + ! Original FATES crown depth heigh used for hydraulics + ! crown_depth = min(height,0.1_r8) + + crown_depth = prt_params%crown(ft) * height + return end subroutine CrownDepth diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 98aaad6488..428f5f7d7a 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -85,8 +85,6 @@ module FatesPlantHydraulicsMod use FatesHydraulicsMemMod, only: aroot_p_media use FatesHydraulicsMemMod, only: rhiz_p_media use FatesHydraulicsMemMod, only: nlevsoi_hyd_max - use FatesHydraulicsMemMod, only: cohort_recruit_water_layer - use FatesHydraulicsMemMod, only: recruit_water_avail_layer use FatesHydraulicsMemMod, only: rwccap, rwcft use FatesHydraulicsMemMod, only: ignore_layer1 @@ -234,8 +232,8 @@ module FatesPlantHydraulicsMod ! The maximum allowable water balance error over a plant-soil continuum ! for a given step [kgs] (0.1 mg) - real(r8), parameter :: max_wb_step_err = 1.e-7_r8 - + real(r8), parameter :: max_wb_step_err = 2.e-7_r8 ! original is 1.e-7_r8, Junyan changed to 2.e-7_r8 + ! ! !PUBLIC MEMBER FUNCTIONS: public :: AccumulateMortalityWaterStorage @@ -352,7 +350,7 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) ccohort_hydr => ccohort%co_hydr ! This calculates node heights - call UpdatePlantHydrNodes(ccohort_hydr,ccohort%pft,ccohort%hite, & + call UpdatePlantHydrNodes(ccohort,ccohort%pft,ccohort%hite, & sites(s)%si_hydr) ! This calculates volumes and lengths @@ -508,14 +506,23 @@ subroutine InitPlantHydStates(site, cohort) do j=1, site_hydr%nlevrhiz - ! Match the potential of the absorbing root to the inner rhizosphere shell - cohort_hydr%psi_aroot(j) = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,1)) - - ! Calculate the mean total potential (include height) of absorbing roots - ! h_aroot_mean = h_aroot_mean + cohort_hydr%psi_aroot(j) + mpa_per_pa*denh2o*grav_earth*(-site_hydr%zi_rhiz(j)) - - cohort_hydr%th_aroot(j) = wrfa%p%th_from_psi(cohort_hydr%psi_aroot(j)) - cohort_hydr%ftc_aroot(j) = wkfa%p%ftc_from_psi(cohort_hydr%psi_aroot(j)) + ! Checking apperance of roots. Only proceed if there are roots in that layer + if(cohort_hydr%l_aroot_layer(j) > nearzero) then + + ! Match the potential of the absorbing root to the inner rhizosphere shell + cohort_hydr%psi_aroot(j) = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,1)) + + ! Calculate the mean total potential (include height) of absorbing roots + ! h_aroot_mean = h_aroot_mean + cohort_hydr%psi_aroot(j) + mpa_per_pa*denh2o*grav_earth*(-site_hydr%zi_rhiz(j)) + + cohort_hydr%th_aroot(j) = wrfa%p%th_from_psi(cohort_hydr%psi_aroot(j)) + cohort_hydr%ftc_aroot(j) = wkfa%p%ftc_from_psi(cohort_hydr%psi_aroot(j)) + else + cohort_hydr%psi_aroot(j) = psi_aroot_init + cohort_hydr%th_aroot(j) = 0 + + end if + end do else @@ -644,7 +651,7 @@ end subroutine UpdatePlantPsiFTCFromTheta ! ===================================================================================== - subroutine UpdatePlantHydrNodes(ccohort_hydr,ft,plant_height,csite_hydr) + subroutine UpdatePlantHydrNodes(ccohort,ft,plant_height,csite_hydr) ! -------------------------------------------------------------------------------- ! This subroutine calculates the nodal heights critical to hydraulics in the plant @@ -661,13 +668,14 @@ subroutine UpdatePlantHydrNodes(ccohort_hydr,ft,plant_height,csite_hydr) ! -------------------------------------------------------------------------------- ! Arguments - type(ed_cohort_hydr_type), intent(inout) :: ccohort_hydr - integer,intent(in) :: ft ! plant functional type index - real(r8), intent(in) :: plant_height ! [m] - type(ed_site_hydr_type), intent(in) :: csite_hydr + type(ed_cohort_type), intent(inout) :: ccohort + integer,intent(in) :: ft ! plant functional type index + real(r8), intent(in) :: plant_height ! [m] + type(ed_site_hydr_type), intent(in) :: csite_hydr ! Locals + type(ed_cohort_hydr_type), pointer :: ccohort_hydr integer :: nlevrhiz ! number of rhizosphere layers real(r8) :: roota ! root profile parameter a zeng2001_crootfr real(r8) :: rootb ! root profile parameter b zeng2001_crootfr @@ -679,15 +687,21 @@ subroutine UpdatePlantHydrNodes(ccohort_hydr,ft,plant_height,csite_hydr) real(r8) :: cumul_rf ! cumulative root distribution where depth is determined [-] real(r8) :: z_cumul_rf ! depth at which cumul_rf occurs [m] integer :: k ! Loop counter for compartments + real(r8) :: z_fr ! Maximum rooting depth of the plant [m] + ccohort_hydr => ccohort%co_hydr + + ! Crown Nodes ! in special case where n_hypool_leaf = 1, the node height of the canopy ! water pool is 1/2 the distance from the bottom of the canopy to the top of the tree roota = prt_params%fnrt_prof_a(ft) rootb = prt_params%fnrt_prof_b(ft) nlevrhiz = csite_hydr%nlevrhiz - call CrownDepth(plant_height,crown_depth) + !call CrownDepth(plant_height,ft,crown_depth) + crown_depth = min(plant_height,0.1_r8) + dz_canopy = crown_depth / real(n_hypool_leaf,r8) do k=1,n_hypool_leaf ccohort_hydr%z_lower_ag(k) = plant_height - dz_canopy*real(k,r8) @@ -707,10 +721,18 @@ subroutine UpdatePlantHydrNodes(ccohort_hydr,ft,plant_height,csite_hydr) ccohort_hydr%z_lower_ag(k) = ccohort_hydr%z_upper_ag(k) - dz_stem enddo + call MaximumRootingDepth(ccohort%dbh,ft,csite_hydr%zi_rhiz(nlevrhiz),z_fr) + ! Transporting Root Node depth [m] (negative from surface) - call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & + call bisect_rootfr(roota, rootb, z_fr, 0._r8, 1.E10_r8, & 0.001_r8, 0.001_r8, 0.5_r8, z_cumul_rf) + + if(z_cumul_rf > csite_hydr%zi_rhiz(nlevrhiz) ) then + print*,"z_cumul_rf > zi_rhiz(nlevrhiz)?",z_cumul_rf,csite_hydr%zi_rhiz(nlevrhiz) + stop + end if + z_cumul_rf = min(z_cumul_rf, abs(csite_hydr%zi_rhiz(nlevrhiz))) ccohort_hydr%z_node_troot = -z_cumul_rf @@ -767,7 +789,7 @@ subroutine UpdateSizeDepPlantHydProps(currentSite,ccohort,bc_in) call SavePreviousCompartmentVolumes(ccohort_hydr) ! This updates all of the z_node positions - call UpdatePlantHydrNodes(ccohort_hydr,ft,ccohort%hite,currentSite%si_hydr) + call UpdatePlantHydrNodes(ccohort,ft,ccohort%hite,currentSite%si_hydr) ! This updates plant compartment volumes, lengths and ! maximum conductances. Make sure for already @@ -826,8 +848,9 @@ subroutine UpdatePlantHydrLenVol(ccohort,site_hydr) real(r8) :: crown_depth ! Depth of the plant's crown [m] real(r8) :: norm ! total root fraction used <1 integer :: nlevrhiz ! number of rhizosphere levels - - + real(r8) :: dbh ! the dbh of current cohort [cm] + real(r8) :: z_fr ! rooting depth of a cohort [cm] + ! We allow the transporting root to donate a fraction of its volume to the absorbing ! roots to help mitigate numerical issues due to very small volumes. This is the ! fraction the transporting roots donate to those layers @@ -840,15 +863,16 @@ subroutine UpdatePlantHydrLenVol(ccohort,site_hydr) ccohort_hydr => ccohort%co_hydr ft = ccohort%pft - nlevrhiz = site_hydr%nlevrhiz - leaf_c = ccohort%prt%GetState(leaf_organ, carbon12_element) - sapw_c = ccohort%prt%GetState(sapw_organ, carbon12_element) - fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element) - struct_c = ccohort%prt%GetState(struct_organ, carbon12_element) + nlevrhiz = site_hydr%nlevrhiz + leaf_c = ccohort%prt%GetState(leaf_organ, carbon12_element) + sapw_c = ccohort%prt%GetState(sapw_organ, carbon12_element) + fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element) + struct_c = ccohort%prt%GetState(struct_organ, carbon12_element) + roota = prt_params%fnrt_prof_a(ft) + rootb = prt_params%fnrt_prof_b(ft) - roota = prt_params%fnrt_prof_a(ft) - rootb = prt_params%fnrt_prof_b(ft) + ! Leaf Volumes ! ----------------------------------------------------------------------------------- @@ -904,7 +928,8 @@ subroutine UpdatePlantHydrLenVol(ccohort,site_hydr) ! alternative cross section calculation ! a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * ccohort%hite ) * 1.e-4_r8 - call CrownDepth(ccohort%hite,crown_depth) + !call CrownDepth(ccohort%hite,ft,crown_depth) + crown_depth = min(ccohort%hite,0.1_r8) z_stem = ccohort%hite - crown_depth v_sapwood = a_sapwood * z_stem ! + 0.333_r8*a_sapwood*crown_depth ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem @@ -939,15 +964,30 @@ subroutine UpdatePlantHydrLenVol(ccohort,site_hydr) ! Partition the total absorbing root lengths and volumes into the active soil layers ! We have a condition, where we may ignore the first layer ! ------------------------------------------------------------------------------ + ! Further, incorporate maximum rooting depth parameterization into these + ! calculations. + + + call MaximumRootingDepth(ccohort%dbh,ft,site_hydr%zi_rhiz(nlevrhiz),z_fr) + norm = 1._r8 - & - zeng2001_crootfr(roota, rootb,site_hydr%zi_rhiz(1)-site_hydr%dz_rhiz(1), site_hydr%zi_rhiz(nlevrhiz)) + zeng2001_crootfr(roota, rootb,site_hydr%zi_rhiz(1)-site_hydr%dz_rhiz(1), z_fr ) do j=1,nlevrhiz - rootfr = norm*(zeng2001_crootfr(roota, rootb, site_hydr%zi_rhiz(j),site_hydr%zi_rhiz(nlevrhiz)) - & - zeng2001_crootfr(roota, rootb, site_hydr%zi_rhiz(j)-site_hydr%dz_rhiz(j),site_hydr%zi_rhiz(nlevrhiz))) + rootfr = norm*(zeng2001_crootfr(roota, rootb, site_hydr%zi_rhiz(j),z_fr) - & + zeng2001_crootfr(roota, rootb, site_hydr%zi_rhiz(j)-site_hydr%dz_rhiz(j),z_fr )) + if(debug)then + write(fates_log(),*) 'check rooting depth of cohort ' + write(fates_log(),*) 'dbh: ',ccohort%dbh,' sice class: ',ccohort%size_class + write(fates_log(),*) 'site_hydr%dz_rhiz(j) is: ', site_hydr%dz_rhiz(j) + write(fates_log(),*) 'z_max cohort: ',z_fr + write(fates_log(),*) 'layer: ',j,' depth (m): ',site_hydr%zi_rhiz(j),' rooting fraction:',rootfr + write(fates_log(),*) 'End of Junyan check' + end if + ccohort_hydr%l_aroot_layer(j) = rootfr*l_aroot_tot ! This is a hybrid absorbing root and transporting root volume @@ -1025,13 +1065,15 @@ subroutine UpdateSizeDepPlantHydStates(currentSite,ccohort) do j=1,currentSite%si_hydr%nlevrhiz - th_uncorr = ccohort_hydr%th_aroot(j) * & - ccohort_hydr%v_aroot_layer_init(j)/ccohort_hydr%v_aroot_layer(j) - ccohort_hydr%th_aroot(j) = constrain_water_contents(th_uncorr, small_theta_num, ft, aroot_p_media) - - csite_hydr%h2oveg_growturn_err = csite_hydr%h2oveg_growturn_err + & - denh2o*cCohort%n*AREA_INV*(ccohort_hydr%th_aroot(j)-th_uncorr)*ccohort_hydr%v_aroot_layer(j) + if (ccohort_hydr%v_aroot_layer(j) > nearzero) then + th_uncorr = ccohort_hydr%th_aroot(j) * & + ccohort_hydr%v_aroot_layer_init(j)/ccohort_hydr%v_aroot_layer(j) + ccohort_hydr%th_aroot(j) = constrain_water_contents(th_uncorr, small_theta_num, ft, aroot_p_media) + + csite_hydr%h2oveg_growturn_err = csite_hydr%h2oveg_growturn_err + & + denh2o*cCohort%n*AREA_INV*(ccohort_hydr%th_aroot(j)-th_uncorr)*ccohort_hydr%v_aroot_layer(j) + end if enddo @@ -1169,7 +1211,7 @@ subroutine FuseCohortHydraulics(currentSite,currentCohort, nextCohort, bc_in, ne call SavePreviousCompartmentVolumes(ccohort_hydr) ! This updates all of the z_node positions - call UpdatePlantHydrNodes(ccohort_hydr,ft,currentCohort%hite,site_hydr) + call UpdatePlantHydrNodes(currentCohort,ft,currentCohort%hite,site_hydr) ! This updates plant compartment volumes, lengths and ! maximum conductances. Make sure for already @@ -1331,7 +1373,8 @@ subroutine InitHydrSites(sites,bc_in) end subroutine InitHydrSites ! =================================================================================== -subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) + +subroutine HydrSiteColdStart(sites, bc_in ) ! Arguments @@ -1623,13 +1666,16 @@ subroutine RecruitWUptake(nsites,sites,bc_in,dtime,recruitflag) endif end do ! site loop - !write(fates_log(),*) 'Calculating recruit water' - !write(fates_log(),*) csite_hydr%recruit_w_uptake + if (debug) then + write(fates_log(),*) 'Calculating recruit uptake' + write(fates_log(),*) sum(csite_hydr%recruit_w_uptake(:)) + endif end subroutine RecruitWUptake !===================================================================================== + subroutine ConstrainRecruitNumber(csite,ccohort, bc_in) ! --------------------------------------------------------------------------- @@ -1645,82 +1691,83 @@ subroutine ConstrainRecruitNumber(csite,ccohort, bc_in) ! Locals type(ed_cohort_hydr_type), pointer :: ccohort_hydr type(ed_site_hydr_type), pointer :: csite_hydr - type(ed_patch_type), pointer :: cpatch + type(ed_patch_type), pointer :: cpatch real(r8) :: tmp1 - real(r8) :: watres_local ! minum water content [m3/m3] - real(r8) :: total_water ! total water in rhizosphere at a specific layer (m^3 ha-1) - real(r8) :: total_water_min ! total minimum water in rhizosphere at a specific layer (m^3) - real(r8) :: rootfr ! fraction of root in different soil layer - real(r8) :: recruitw ! water for newly recruited cohorts (kg water/m2/individual) - real(r8) :: n, nmin ! number of individuals in cohorts + real(r8) :: watres_local ! minum water content [m3/m3] + real(r8) :: total_water ! total water in rhizosphere at a specific layer (m^3 ha-1) + real(r8) :: total_water_min ! total minimum water in rhizosphere at a specific layer (m^3) + real(r8) :: rootfr ! fraction of root in different soil layer + real(r8) :: recruitw ! water for newly recruited cohorts (kg water/m2/individual) + real(r8) :: n, nmin ! number of individuals in cohorts real(r8) :: sum_l_aroot integer :: s, j, ft - integer :: el ! element loop index - integer :: element_id ! global element identifier index - real(r8) :: leaf_m, store_m, sapw_m ! Element mass in organ tissues - real(r8) :: fnrt_m, struct_m, repro_m ! Element mass in organ tissues - + integer :: el ! element loop index + integer :: element_id ! global element identifier index + real(r8) :: leaf_m, store_m, sapw_m ! Element mass in organ tissues + real(r8) :: fnrt_m, struct_m, repro_m ! Element mass in organ tissues - cpatch => ccohort%patchptr + cpatch => ccohort%patchptr csite_hydr => csite%si_hydr ccohort_hydr =>ccohort%co_hydr recruitw = (sum(ccohort_hydr%th_ag(:)*ccohort_hydr%v_ag(:)) + & ccohort_hydr%th_troot*ccohort_hydr%v_troot + & sum(ccohort_hydr%th_aroot(:)*ccohort_hydr%v_aroot_layer(:)))* & denh2o + sum_l_aroot = sum(ccohort_hydr%l_aroot_layer(:)) do j=1,csite_hydr%nlevrhiz - cohort_recruit_water_layer(j) = recruitw*ccohort_hydr%l_aroot_layer(j)/sum_l_aroot + csite_hydr%cohort_recruit_water_layer(j) = recruitw*ccohort_hydr%l_aroot_layer(j)/sum_l_aroot end do do j=1,csite_hydr%nlevrhiz - watres_local = csite_hydr%wrf_soil(j)%p%th_from_psi(bc_in%smpmin_si*denh2o*grav_earth*m_per_mm*mpa_per_pa) + watres_local = csite_hydr%wrf_soil(j)%p%th_from_psi(bc_in%smpmin_si*denh2o*grav_earth*m_per_mm*mpa_per_pa) + total_water = sum(csite_hydr%v_shell(j,:)*csite_hydr%h2osoi_liqvol_shell(j,:)) total_water_min = sum(csite_hydr%v_shell(j,:)*watres_local) - + !assumes that only 50% is available for recruit water.... - recruit_water_avail_layer(j)=0.5_r8*max(0.0_r8,total_water-total_water_min) - + csite_hydr%recruit_water_avail_layer(j)=0.5_r8*max(0.0_r8,total_water-total_water_min) + end do - + nmin = 1.0e+36 do j=1,csite_hydr%nlevrhiz - if(cohort_recruit_water_layer(j)>0.0_r8) then - n = recruit_water_avail_layer(j)/cohort_recruit_water_layer(j) + if(csite_hydr%cohort_recruit_water_layer(j)>nearzero) then + n = csite_hydr%recruit_water_avail_layer(j)/csite_hydr%cohort_recruit_water_layer(j) nmin = min(n, nmin) endif end do - ! If the minimum number of plants that are recruitable due to water - ! limitations, is less than what is currently recruitable (due to - ! carbon-nitrogen-phosphorus availability), then we apply a reduction. - ! We also have to add back in what had been taken, to the germination - ! seed pool - if(nmin < ccohort%n) then + ! If the minimum number of plants that are recruitable due to water + ! limitations, is less than what is currently recruitable (due to + ! carbon-nitrogen-phosphorus availability), then we apply a reduction. + ! We also have to add back in what had been taken, to the germination + ! seed pool + if(nmin < ccohort%n) then - do el = 1,num_elements + do el = 1,num_elements - element_id = element_list(el) + element_id = element_list(el) - leaf_m = ccohort%prt%GetState(leaf_organ, element_id) - store_m = ccohort%prt%GetState(store_organ, element_id) - sapw_m = ccohort%prt%GetState(sapw_organ, element_id) - fnrt_m = ccohort%prt%GetState(fnrt_organ, element_id) - struct_m = ccohort%prt%GetState(struct_organ, element_id) - repro_m = ccohort%prt%GetState(repro_organ, element_id) + leaf_m = ccohort%prt%GetState(leaf_organ, element_id) + store_m = ccohort%prt%GetState(store_organ, element_id) + sapw_m = ccohort%prt%GetState(sapw_organ, element_id) + fnrt_m = ccohort%prt%GetState(fnrt_organ, element_id) + struct_m = ccohort%prt%GetState(struct_organ, element_id) + repro_m = ccohort%prt%GetState(repro_organ, element_id) - cpatch%litter(el)%seed_germ(ccohort%pft) = cpatch%litter(el)%seed_germ(ccohort%pft) + & - (ccohort%n-nmin)/cpatch%area * & - (leaf_m+store_m+sapw_m+fnrt_m+struct_m+repro_m) + cpatch%litter(el)%seed_germ(ccohort%pft) = cpatch%litter(el)%seed_germ(ccohort%pft) + & + (ccohort%n-nmin)/cpatch%area * & + (leaf_m+store_m+sapw_m+fnrt_m+struct_m+repro_m) - end do - ccohort%n = nmin - end if + end do + ccohort%n = nmin + end if - return - end subroutine ConstrainRecruitNumber + return +end subroutine ConstrainRecruitNumber ! ===================================================================================== @@ -1744,7 +1791,7 @@ end subroutine SavePreviousRhizVolumes subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) ! - ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes. + ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes of the site. ! As fine root biomass (and thus absorbing root length) increases, this characteristic ! rhizosphere shrinks even though the total volume of soil tapped by fine roots remains ! the same. @@ -1777,7 +1824,9 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) csite_hydr => currentSite%si_hydr nlevrhiz = csite_hydr%nlevrhiz + ! Note, here is where the site level soil depth/layer is set ! update cohort-level root length density and accumulate it across cohorts and patches to the column level + csite_hydr%l_aroot_layer(:) = 0._r8 cPatch => currentSite%youngest_patch do while(associated(cPatch)) @@ -1791,12 +1840,12 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) enddo !patch ! update outer radii of column-level rhizosphere shells (same across patches and cohorts) + ! Provisions are made inside shellGeom() for layers with no roots do j = 1,nlevrhiz - ! proceed only if l_aroot_coh has changed - ! if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then + call shellGeom( csite_hydr%l_aroot_layer(j), csite_hydr%rs1(j), AREA, csite_hydr%dz_rhiz(j), & csite_hydr%r_out_shell(j,:), csite_hydr%r_node_shell(j,:),csite_hydr%v_shell(j,:)) - ! end if !has l_aroot_layer changed? + enddo @@ -1818,7 +1867,8 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) hksat_s = bc_in%hksat_sisl(j_bc) * m_per_mm * 1._r8/grav_earth * pa_per_mpa ! proceed only if the total absorbing root length (site-level) has changed in this layer - if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then + if( (csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j)) .and. & + csite_hydr%l_aroot_layer(j)>nearzero ) then ! Set the max conductance on the inner shell first. If the node radius ! on the shell is smaller than the root radius, just set the max conductance @@ -1842,9 +1892,6 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) log(csite_hydr%r_out_shell(j,k)/csite_hydr%r_node_shell(j,k ))*hksat_s enddo ! loop over rhizosphere shells - - - end if !has l_aroot_layer changed? enddo ! loop over soil layers @@ -1927,7 +1974,7 @@ subroutine UpdateSizeDepRhizHydStates(currentSite, bc_in) csite_hydr => currentSite%si_hydr - if(.false.) then + bypass_routine: if(.false.) then do j = 1, csite_hydr%nlevrhiz ! proceed only if l_aroot_coh has changed @@ -2064,8 +2111,9 @@ subroutine UpdateSizeDepRhizHydStates(currentSite, bc_in) end if enddo - end if !nshell > 1 + end if bypass_routine !nshell > 1 + end subroutine UpdateSizeDepRhizHydStates ! ==================================================================================== @@ -2181,9 +2229,14 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) csite_hydr => sites(s)%si_hydr + ! If there are just no plants in this site, don't bother shuffling water + if( sum(csite_hydr%l_aroot_layer) <= nearzero ) cycle + do j = 1,csite_hydr%nlevrhiz j_bc = j+csite_hydr%i_rhiz_t-1 + if (csite_hydr%l_aroot_layer(j) <= nearzero ) cycle + cumShellH2O=sum(csite_hydr%h2osoi_liqvol_shell(j,:) *csite_hydr%v_shell(j,:)) * denh2o*AREA_INV dwat_kgm2 = bc_in(s)%h2o_liq_sisl(j_bc) - cumShellH2O @@ -2481,7 +2534,10 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! |_____| | | 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(use_2d_hydrosolve) then call MatSolve2D(bc_in(s),site_hydr,ccohort,ccohort_hydr, & @@ -2585,27 +2641,51 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) root_flux = -sum(dth_layershell_col(1:site_hydr%nlevrhiz,:)*site_hydr%v_shell(:,:))*denh2o*AREA_INV + if(debug)then + write(fates_log(),*) 'root_flux: ', root_flux + end if + + ! Since not all layers have roots, we filter, therefore zero fluxes + bc_out(s)%qflx_soil2root_sisl(:) = 0._r8 + bc_out(s)%qflx_ro_sisl(:) = 0._r8 + do j=1,site_hydr%nlevrhiz j_bc = j+site_hydr%i_rhiz_t-1 - ! Update the site-level state variable - ! rhizosphere shell water content [m3/m3] - site_hydr%h2osoi_liqvol_shell(j,:) = site_hydr%h2osoi_liqvol_shell(j,:) + & - dth_layershell_col(j,:) + ! loginfo + if (debug) then + write(fates_log(),*) 'hydraulics_bc() position I' + write(fates_log(),*) 'layer: ', j + write(fates_log(),*) 'dth_layershell_col(j,:):', dth_layershell_col(j,:) + write(fates_log(),*) 'site_hydr%v_shell(j,:):', site_hydr%v_shell(j,:) + write(fates_log(),*) 'site_hydr%h2osoi_liqvol_shell: ', site_hydr%h2osoi_liqvol_shell(j,:) + write(fates_log(),*) 'dth_layershell_col(j,:) ', dth_layershell_col(j,:) + write(fates_log(),*) 'site_hydr%l_aroot_layer(j): ' , site_hydr%l_aroot_layer(j) + endif + + if (site_hydr%l_aroot_layer(j) > nearzero) then - bc_out(s)%qflx_soil2root_sisl(j_bc) = & - -(sum(dth_layershell_col(j,:)*site_hydr%v_shell(j,:))*denh2o*AREA_INV/dtime) + & - site_hydr%recruit_w_uptake(j) + ! Update the site-level state variable + ! rhizosphere shell water content [m3/m3] + site_hydr%h2osoi_liqvol_shell(j,:) = site_hydr%h2osoi_liqvol_shell(j,:) + & + dth_layershell_col(j,:) - ! Save the amount of liquid soil water known to the model after root uptake - ! This calculation also assumes that 1mm of water is 1kg - site_hydr%h2osoi_liq_prev(j) = bc_in(s)%h2o_liq_sisl(j_bc) - & - dtime*bc_out(s)%qflx_soil2root_sisl(j_bc) + bc_out(s)%qflx_soil2root_sisl(j_bc) = & + -(sum(dth_layershell_col(j,:)*site_hydr%v_shell(j,:))*denh2o*AREA_INV/dtime) + & + site_hydr%recruit_w_uptake(j) + ! Save the amount of liquid soil water known to the model after root uptake + ! This calculation also assumes that 1mm of water is 1kg + site_hydr%h2osoi_liq_prev(j) = bc_in(s)%h2o_liq_sisl(j_bc) - & + dtime*bc_out(s)%qflx_soil2root_sisl(j_bc) + + + end if + ! We accept that it is possible for gravity to push ! water into saturated soils, particularly at night when ! transpiration has stopped. In the real world, the water @@ -2614,24 +2694,25 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! this, if soils are pushed beyond saturation minus a small buffer ! then we remove that excess, send it to a runoff pool, and ! fix the node's water content to the saturation minus buffer value - + site_runoff = 0._r8 if(purge_supersaturation) then do i = 1,nshell if(site_hydr%h2osoi_liqvol_shell(j,i)>(bc_in(s)%watsat_sisl(j_bc)-thsat_buff)) then - + ! [m3/m3] * [kg/m3] * [m3/site] * [site/m2] => [kg/m2] site_runoff = site_runoff + & (site_hydr%h2osoi_liqvol_shell(j,i)-(bc_in(s)%watsat_sisl(j_bc)-thsat_buff)) * & site_hydr%v_shell(j,i)*AREA_INV*denh2o - + site_hydr%h2osoi_liqvol_shell(j,i) = bc_in(s)%watsat_sisl(j_bc)-thsat_buff end if end do - + bc_out(s)%qflx_ro_sisl(j_bc) = site_runoff/dtime end if + enddo @@ -2646,14 +2727,13 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) delta_soil_storage = sum(site_hydr%h2osoi_liqvol_shell(:,:) * & site_hydr%v_shell(:,:)) * denh2o * AREA_INV - prev_h2osoil - if(abs(delta_plant_storage - (root_flux - transp_flux)) > error_thresh ) then - + if(abs(delta_plant_storage - (root_flux - transp_flux)) > error_thresh ) then write(fates_log(),*) 'Site plant water balance does not close' - write(fates_log(),*) 'balance error: ',abs(delta_plant_storage - (root_flux - transp_flux)) write(fates_log(),*) 'delta plant storage: ',delta_plant_storage,' [kg/m2]' write(fates_log(),*) 'integrated root flux: ',root_flux,' [kg/m2]' write(fates_log(),*) 'transpiration flux: ',transp_flux,' [kg/m2]' write(fates_log(),*) 'end storage: ',site_hydr%h2oveg + write(fates_log(),*) 'pre_h2oveg', prev_h2oveg call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -2683,7 +2763,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! Now check on total error if( abs(wb_check_site) > 1.e-4_r8 ) then - write(fates_log(),*) 'FATES hydro water balance is not so great [kg/m2]' + write(fates_log(),*) 'FATES hydro water balance does not add up [kg/m2]' write(fates_log(),*) 'site_hydr%errh2o_hyd: ',wb_check_site write(fates_log(),*) 'delta_plant_storage: ',delta_plant_storage write(fates_log(),*) 'delta_soil_storage: ',delta_soil_storage @@ -2978,61 +3058,73 @@ subroutine OrderLayersForSolve1D(site_hydr,cohort,cohort_hydr,ordered, kbg_layer do j=1,site_hydr%nlevrhiz - ! Path is between the absorbing root - ! and the first rhizosphere shell nodes - ! Special case. Maximum conductance depends on the - ! potential gradient (same elevation, no geopotential - ! required. + if(cohort_hydr%l_aroot_layer(j)>nearzero)then - psi_inner_shell = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,1)) + ! Path is between the absorbing root + ! and the first rhizosphere shell nodes + ! Special case. Maximum conductance depends on the + ! potential gradient (same elevation, no geopotential + ! required. - ! Note, since their is no elevation difference between - ! the absorbing root and its layer, no need to calc - ! diff in total, just matric is fine [MPa] - if(cohort_hydr%psi_aroot(j) < psi_inner_shell) then - kmax_aroot = cohort_hydr%kmax_aroot_radial_in(j) - else - kmax_aroot = cohort_hydr%kmax_aroot_radial_out(j) - end if + psi_inner_shell = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,1)) - ! Get matric potential [Mpa] of the absorbing root - psi_aroot = wrf_plant(aroot_p_media,ft)%p%psi_from_th(cohort_hydr%th_aroot(j)) + ! Note, since their is no elevation difference between + ! the absorbing root and its layer, no need to calc + ! diff in total, just matric is fine [MPa] + if(cohort_hydr%psi_aroot(j) < psi_inner_shell) then + kmax_aroot = cohort_hydr%kmax_aroot_radial_in(j) + else + kmax_aroot = cohort_hydr%kmax_aroot_radial_out(j) + end if - ! Get Fraction of Total Conductivity [-] of the absorbing root - ftc_aroot = wkf_plant(aroot_p_media,ft)%p%ftc_from_psi(cohort_hydr%psi_aroot(j)) + ! Get matric potential [Mpa] of the absorbing root + psi_aroot = wrf_plant(aroot_p_media,ft)%p%psi_from_th(cohort_hydr%th_aroot(j)) - ! Calculate total effective conductance over path [kg s-1 MPa-1] - ! from absorbing root node to 1st rhizosphere shell - r_bg = 1._r8/(kmax_aroot*ftc_aroot) + ! Get Fraction of Total Conductivity [-] of the absorbing root + ftc_aroot = wkf_plant(aroot_p_media,ft)%p%ftc_from_psi(cohort_hydr%psi_aroot(j)) - ! Path is across the upper an lower rhizosphere comparment - ! on each side of the nodes. Since there is no flow across the outer - ! node to the edge, we ignore that last half compartment - aroot_frac_plant = cohort_hydr%l_aroot_layer(j)/site_hydr%l_aroot_layer(j) + ! Calculate total effective conductance over path [kg s-1 MPa-1] + ! from absorbing root node to 1st rhizosphere shell + r_bg = 1._r8/(kmax_aroot*ftc_aroot) - do k = 1,nshell + ! Path is across the upper an lower rhizosphere comparment + ! on each side of the nodes. Since there is no flow across the outer + ! node to the edge, we ignore that last half compartment + aroot_frac_plant = cohort_hydr%l_aroot_layer(j)/site_hydr%l_aroot_layer(j) - kmax_up = site_hydr%kmax_upper_shell(j,k)*aroot_frac_plant - kmax_lo = site_hydr%kmax_lower_shell(j,k)*aroot_frac_plant + do k = 1,nshell - psi_shell = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,k)) + kmax_up = site_hydr%kmax_upper_shell(j,k)*aroot_frac_plant + kmax_lo = site_hydr%kmax_lower_shell(j,k)*aroot_frac_plant - ftc_shell = site_hydr%wkf_soil(j)%p%ftc_from_psi(psi_shell) + psi_shell = site_hydr%wrf_soil(j)%p%psi_from_th(site_hydr%h2osoi_liqvol_shell(j,k)) - r_bg = r_bg + 1._r8/(kmax_up*ftc_shell) - if(k=5, rhizosphere z_node(i) = -site_hydr%zi_rhiz(ilayer) ! The volume of the Rhizosphere for a single plant v_node(i) = site_hydr%v_shell(ilayer,ishell)*aroot_frac_plant th_node_init(i) = site_hydr%h2osoi_liqvol_shell(ilayer,ishell) + if (th_node_init(i) < -nearzero) then + write(fates_log(),*) 'ImTaylorSolve1D(), print out shell theta' + write(fates_log(),*) 'layer: ',ilayer, 'shell:', ishell + write(fates_log(),*) 'th_node_init(i) is: ', th_node_init(i) + write(fates_log(),*) 'th_node_init(i) is: ', th_node_init(i) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if end do + ! Outer iteration loop ! This cuts timestep in half and resolve the solution with smaller substeps ! This loop is cleared when the model has found a solution solution_found = .false. iter = 0 - do while( .not.solution_found ) + solution_iteration: do while( .not.solution_found ) ! Gracefully quit if too many iterations have been used if(iter>max_iter)then @@ -3360,7 +3463,7 @@ subroutine ImTaylorSolve1D(site_hydr,cohort,cohort_hydr,dtime,q_top, & dftc_dpsi = site_hydr%wkf_soil(ilayer)%p%dftcdpsi_from_psi(psi_node(i)) dftc_dtheta_node(i) = dftc_dpsi * dpsi_dtheta_node(i) end do - + !-------------------------------------------------------------------------------- ! Part 2. Effective conductances over the path-length and Flux terms ! over the node-to-node paths @@ -3650,7 +3753,7 @@ subroutine ImTaylorSolve1D(site_hydr,cohort,cohort_hydr,dtime,q_top, & iter=iter+1 - end do + end do solution_iteration ! ----------------------------------------------------------- ! Do a final check on water balance error sumed over sub-steps @@ -3728,19 +3831,12 @@ subroutine ImTaylorSolve1D(site_hydr,cohort,cohort_hydr,dtime,q_top, & ! after all cohort-layers are complete. This allows each cohort ! to experience the same water conditions (for good or bad). - if(site_hydr%l_aroot_layer(ilayer) ilayer) + enddo loop_root_layers end associate return @@ -3812,6 +3908,9 @@ subroutine Report1DError(cohort, site_hydr, ilayer, z_node, v_node, & write(fates_log(),*) 'layer: ',ilayer write(fates_log(),*) 'wb_step_err = ',(q_top_eff*dt_step) - (w_tot_beg-w_tot_end) + write(fates_log(),*) 'q_top_eff*dt_step = ',q_top_eff*dt_step + write(fates_log(),*) 'w_tot_beg = ',w_tot_beg + write(fates_log(),*) 'w_tot_end = ',w_tot_end write(fates_log(),*) 'leaf water: ',leaf_water,' kg/plant' write(fates_log(),*) 'stem_water: ',stem_water,' kg/plant' write(fates_log(),*) 'troot_water: ',troot_water @@ -3865,13 +3964,8 @@ subroutine Report1DError(cohort, site_hydr, ilayer, z_node, v_node, & write(fates_log(),*) 'r4:',v_node(8) write(fates_log(),*) ' ' write(fates_log(),*) 'r5:',v_node(9) - write(fates_log(),*) 'inner shell kmaxs: ',site_hydr%kmax_lower_shell(:,1)*aroot_frac_plant - - - - deallocate(psi_node) deallocate(h_node) @@ -4145,7 +4239,43 @@ end subroutine RecruitWaterStorage ! Utility Functions ! ===================================================================================== -subroutine bisect_rootfr(a, b, lower_init, upper_init, xtol, ytol, crootfr, x_new) +subroutine MaximumRootingDepth(dbh,ft,z_max_soil,z_fr) + + ! --------------------------------------------------------------------------------- + ! Calculate the maximum rooting depth of the plant. + ! + ! This is an exponential which is constrained by the maximum soil depth: + ! site_hydr%zi_rhiz(nlevrhiz) + ! The dynamic root growth model by Junyan Ding, June 9, 2021 + ! --------------------------------------------------------------------------------- + + real(r8),intent(in) :: dbh ! Plant dbh + integer,intent(in) :: ft ! Funtional type index + real(r8),intent(in) :: z_max_soil ! Maximum depth of soil (pos convention) [m] + real(r8),intent(out) :: z_fr ! Maximum depth of plant's roots + ! (pos convention) [m] + + real(r8) :: dbh_rel ! Relative dbh of plant between the diameter at which we + ! define the shallowest rooting depth (dbh_0) and the diameter + ! at which we define the deepest rooting depth (dbh_max) + + associate( & + dbh_max => prt_params%allom_zroot_max_dbh(ft), & + dbh_0 => prt_params%allom_zroot_min_dbh(ft), & + z_fr_max => prt_params%allom_zroot_max_z(ft), & + z_fr_0 => prt_params%allom_zroot_min_z(ft), & + frk => prt_params%allom_zroot_k(ft)) + + dbh_rel = min(1._r8,(max(dbh,dbh_0) - dbh_0)/(dbh_max - dbh_0)) + + z_fr = min(z_max_soil, z_fr_max/(1._r8 + ((z_fr_max-z_fr_0)/z_fr_0)*exp(-frk*dbh_rel))) + + end associate + return +end subroutine MaximumRootingDepth + + +subroutine bisect_rootfr(a, b, z_max, lower_init, upper_init, xtol, ytol, crootfr, x_new) ! ! !DESCRIPTION: Bisection routine for getting the inverse of the cumulative root ! distribution. No analytical soln bc crootfr ~ exp(ax) + exp(bx). @@ -4153,7 +4283,8 @@ subroutine bisect_rootfr(a, b, lower_init, upper_init, xtol, ytol, crootfr, x_ne ! !USES: ! ! !ARGUMENTS - real(r8) , intent(in) :: a, b ! pft root distribution constants + real(r8) , intent(in) :: a, b ! pft root distribution constants + real(r8) , intent(in) :: z_max ! maximum rooting depth real(r8) , intent(in) :: lower_init ! lower bound of initial x estimate [m] real(r8) , intent(in) :: upper_init ! upper bound of initial x estimate [m] real(r8) , intent(in) :: xtol ! error tolerance for x_new [m] @@ -4175,12 +4306,12 @@ subroutine bisect_rootfr(a, b, lower_init, upper_init, xtol, ytol, crootfr, x_ne lower = lower_init upper = upper_init - f_lo = zeng2001_crootfr(a, b, lower) - crootfr - f_hi = zeng2001_crootfr(a, b, upper) - crootfr + f_lo = zeng2001_crootfr(a, b, lower, z_max) - crootfr + f_hi = zeng2001_crootfr(a, b, upper, z_max) - crootfr chg = upper - lower do while(abs(chg) .gt. xtol) x_new = 0.5_r8*(lower + upper) - f_new = zeng2001_crootfr(a, b, x_new) - crootfr + f_new = zeng2001_crootfr(a, b, x_new, z_max) - crootfr if(abs(f_new) .le. ytol) then EXIT end if @@ -4214,6 +4345,9 @@ function zeng2001_crootfr(a, b, z, z_max) result(crootfr) ! root fraction. if(present(z_max))then + ! If the soil depth is larger than the maximum rooting depth of the cohort, + ! then the cumulative root fraction of that layer equals that of the maximum rooting depth + crootfr = 1._r8 - .5_r8*(exp(-a*min(z,z_max)) + exp(-b*min(z,z_max))) crootfr_max = 1._r8 - .5_r8*(exp(-a*z_max) + exp(-b*z_max)) crootfr = crootfr/crootfr_max end if @@ -4235,7 +4369,7 @@ end function zeng2001_crootfr ! ===================================================================================== -subroutine shellGeom(l_aroot, rs1, area_site, dz, r_out_shell, r_node_shell, v_shell) +subroutine shellGeom(l_aroot_in, rs1_in, area_site, dz, r_out_shell, r_node_shell, v_shell) ! ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes. ! As fine root biomass (and thus absorbing root length) increases, this characteristic @@ -4246,9 +4380,9 @@ subroutine shellGeom(l_aroot, rs1, area_site, dz, r_out_shell, r_node_shell, v_s ! ! !ARGUMENTS: - real(r8) , intent(in) :: l_aroot ! Total length of absorbing roots - ! for the whole site, this layer (m) - real(r8) , intent(in) :: rs1 ! Fine root radius (m) + real(r8) , intent(in) :: l_aroot_in ! Total length of absorbing roots + ! for the whole site, this layer (m) + real(r8) , intent(in) :: rs1_in ! Fine root radius (m) real(r8) , intent(in) :: area_site ! Area of site (10,000 m2) real(r8) , intent(in) :: dz ! Width of current soil layer (m) real(r8) , intent(out) :: r_out_shell(:) ! Outer radius of each shell (m) @@ -4257,13 +4391,44 @@ subroutine shellGeom(l_aroot, rs1, area_site, dz, r_out_shell, r_node_shell, v_s ! for this layer ! ! !LOCAL VARIABLES: + real(r8) :: l_aroot ! effective length of absorbing root (m/layer) + real(r8) :: rs1 ! effective fine root ratius (m) integer :: k ! rhizosphere shell indicies integer :: nshells ! We don't use the global because of unit testing + + + ! When we have no roots, we may choose to use a nominal + ! value of 1cm per cubic meter to define the rhizosphere shells + ! this "should" help with the transition when roots grow into a layer + ! real(r8), parameter :: nominal_l_aroot = 0.01_r8 ! m/m3 + + !----------------------------------------------------------------------- nshells = size(r_out_shell,dim=1) + + if( l_aroot_in <= nearzero ) then + + ! Generate some nominal values for cases where we have no roots + ! The rational for this is to maintain shells and water contents in those + ! shells similar to what will be experienced when roots start to emerge + ! in these layers, so there will not be a shock to the system + ! Note! All root radii are currently the fine_root_radius const anyway (RGK 10-2021) + ! rs1 = fine_root_radius_const + ! l_aroot = nominal_l_aroot*dz + + r_out_shell(:) = 0._r8 + r_node_shell(:) = 0._r8 + v_shell(:) = 0._r8 + return + + else + rs1 = rs1_in + l_aroot = l_aroot_in + end if + ! update outer radii of column-level rhizosphere shells (same across patches and cohorts) r_out_shell(nshells) = (pi_const*l_aroot/(area_site*dz))**(-0.5_r8) ! eqn(8) S98 if(nshells > 1) then @@ -4297,8 +4462,8 @@ end subroutine shellGeom function xylemtaper(p, dz) result(chi_tapnotap) - ! !ARGUMENTS: - real(r8) , intent(in) :: p ! Taper exponent (see EDPftvar hydr_p_taper) [-] + ! !ARGUMENTS: + real(r8) , intent(in) :: p ! Savage et al. (2010) taper exponent real(r8) , intent(in) :: dz ! hydraulic distance from petiole to node of interest [m] ! ! !LOCAL VARIABLES: @@ -4677,7 +4842,11 @@ subroutine MatSolve2D(bc_in,site_hydr,cohort,cohort_hydr, & ! folume that this plant's rhizosphere accounts forPath is across the upper an lower rhizosphere comparment ! on each side of the nodes. Since there is no flow across the outer ! node to the edge, we ignore that last half compartment - aroot_frac_plant = cohort_hydr%l_aroot_layer(j)/site_hydr%l_aroot_layer(j) + if(cohort_hydr%l_aroot_layer(j)>nearzero)then + aroot_frac_plant = cohort_hydr%l_aroot_layer(j)/site_hydr%l_aroot_layer(j) + else + aroot_frac_plant = 0._r8 + end if do k = 1, n_hypool_aroot + nshell i = i + 1 diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 127dfa43f9..1d08ae2e51 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -46,6 +46,7 @@ module SFMainMod use PRTGenericMod, only : struct_organ use PRTGenericMod, only : SetState use FatesInterfaceTypesMod , only : numpft + use FatesAllometryMod, only : CrownDepth implicit none private @@ -908,7 +909,8 @@ subroutine crown_damage ( currentSite ) type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - + real(r8) :: crown_depth ! Depth of crown in meters + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) @@ -921,20 +923,21 @@ subroutine crown_damage ( currentSite ) if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only ! Flames lower than bottom of canopy. ! c%hite is height of cohort + + call CrownDepth(currentCohort%hite,currentCohort%pft,crown_depth) + if (currentPatch%Scorch_ht(currentCohort%pft) < & - (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft))) then + (currentCohort%hite-crown_depth)) then currentCohort%fraction_crown_burned = 0.0_r8 else ! Flames part of way up canopy. ! Equation 17 in Thonicke et al. 2010. ! flames over bottom of canopy but not over top. if ((currentCohort%hite > 0.0_r8).and.(currentPatch%Scorch_ht(currentCohort%pft) >= & - (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft)))) then + (currentCohort%hite-crown_depth))) then currentCohort%fraction_crown_burned = (currentPatch%Scorch_ht(currentCohort%pft) - & - currentCohort%hite*(1.0_r8 - & - EDPftvarcon_inst%crown(currentCohort%pft)))/(currentCohort%hite* & - EDPftvarcon_inst%crown(currentCohort%pft)) + (currentCohort%hite - crown_depth))/crown_depth else ! Flames over top of canopy. diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index c89e63df98..89c0cf2ea5 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -46,9 +46,7 @@ module EDPftvarcon real(r8), allocatable :: hgt_min(:) ! sapling height m real(r8), allocatable :: dleaf(:) ! leaf characteristic dimension length (m) real(r8), allocatable :: z0mr(:) ! ratio of roughness length of vegetation to height (-) - real(r8), allocatable :: displar(:) ! ratio of displacement height to canopy top height (-) - real(r8), allocatable :: crown(:) ! fraction of the height of the plant - ! that is occupied by crown. For fire model. + real(r8), allocatable :: displar(:) ! ratio of displacement height to canopy top height real(r8), allocatable :: bark_scaler(:) ! scaler from dbh to bark thickness. For fire model. real(r8), allocatable :: crown_kill(:) ! scaler on fire death. For fire model. real(r8), allocatable :: initd(:) ! initial seedling density @@ -331,10 +329,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_fire_crown_depth_frac' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_fire_bark_scaler' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -669,10 +663,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%hgt_min) - name = 'fates_fire_crown_depth_frac' - call fates_params%RetreiveParameterAllocate(name=name, & - data=this%crown) - name = 'fates_fire_bark_scaler' call fates_params%RetreiveParameterAllocate(name=name, & data=this%bark_scaler) @@ -1394,7 +1384,6 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'dleaf = ',EDPftvarcon_inst%dleaf write(fates_log(),fmt0) 'z0mr = ',EDPftvarcon_inst%z0mr write(fates_log(),fmt0) 'displar = ',EDPftvarcon_inst%displar - write(fates_log(),fmt0) 'crown = ',EDPftvarcon_inst%crown write(fates_log(),fmt0) 'bark_scaler = ',EDPftvarcon_inst%bark_scaler write(fates_log(),fmt0) 'crown_kill = ',EDPftvarcon_inst%crown_kill write(fates_log(),fmt0) 'initd = ',EDPftvarcon_inst%initd diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 55b69d78bc..862ddabfaf 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -46,8 +46,8 @@ module FatesHistoryInterfaceMod use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : hlm_model_day use FatesInterfaceTypesMod , only : nlevcoage - - ! FIXME(bja, 2016-10) need to remove CLM dependancy + use FatesAllometryMod , only : CrownDepth + use EDPftvarcon , only : EDPftvarcon_inst use PRTParametersMod , only : prt_params @@ -1775,6 +1775,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8) :: struct_m_net_alloc real(r8) :: repro_m_net_alloc real(r8) :: area_frac + real(r8) :: crown_depth type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort @@ -2219,8 +2220,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) + ccohort%c_area * AREA_INV ! calculate leaf height distribution, assuming leaf area is evenly distributed thru crown depth + call CrownDepth(ccohort%hite,ft,crown_depth) height_bin_max = get_height_index(ccohort%hite) - height_bin_min = get_height_index(ccohort%hite * (1._r8 - EDPftvarcon_inst%crown(ft))) + height_bin_min = get_height_index(ccohort%hite - crown_depth) do i_heightbin = height_bin_min, height_bin_max binbottom = ED_val_history_height_bin_edges(i_heightbin) if (i_heightbin .eq. nlevheight) then @@ -2230,8 +2232,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) endif ! what fraction of a cohort's crown is in this height bin? frac_canopy_in_bin = (min(bintop,ccohort%hite) - & - max(binbottom,ccohort%hite * (1._r8 - EDPftvarcon_inst%crown(ft)))) / & - (ccohort%hite * EDPftvarcon_inst%crown(ft)) + max(binbottom,ccohort%hite-crown_depth)) / & + (crown_depth) ! hio_leaf_height_dist_si_height(io_si,i_heightbin) = & hio_leaf_height_dist_si_height(io_si,i_heightbin) + & diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index 5de1165a16..f971b4f55b 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -70,9 +70,7 @@ module FatesHydraulicsMemMod ! ---------------------------------------------------------------------------------------------- !temporatory variables - real(r8), public :: cohort_recruit_water_layer(nlevsoi_hyd_max) ! the recruit water requirement for a - ! single individual at different layer (kg H2o/m2) - real(r8), public :: recruit_water_avail_layer(nlevsoi_hyd_max) ! the recruit water avaibility from soil (kg H2o/m2) + type, public :: ed_site_hydr_type @@ -185,10 +183,12 @@ module FatesHydraulicsMemMod real(r8), allocatable :: q_flux(:) real(r8), allocatable :: dftc_dpsi_node(:) real(r8), allocatable :: ftc_node(:) - - real(r8), allocatable :: kmax_up(:) real(r8), allocatable :: kmax_dn(:) + + ! Scratch arrays + real(r8) :: cohort_recruit_water_layer(nlevsoi_hyd_max) ! the recruit water requirement for a + real(r8) :: recruit_water_avail_layer(nlevsoi_hyd_max) ! the recruit water avaibility from soil (kg H2o/m2) contains @@ -496,6 +496,9 @@ end subroutine FlushSiteScratch subroutine SetConnections(this) + ! This routine should be updated + ! when new layers are added as plants grow into them? + class(ed_site_hydr_type),intent(inout) :: this integer :: k, j diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 39adab94f6..6d5fa2cb4e 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -137,19 +137,19 @@ variables: fates_allom_stmode:possible_values = "1: target storage proportional to trimmed maximum leaf biomass." ; double fates_allom_zroot_k(fates_pft) ; fates_allom_zroot_k:units = "unitless" ; - fates_allom_zroot_k:long_name = "scale coefficient of logistic rooting depth model (NOT USED)" ; + fates_allom_zroot_k:long_name = "scale coefficient of logistic rooting depth model" ; double fates_allom_zroot_max_dbh(fates_pft) ; fates_allom_zroot_max_dbh:units = "cm" ; - fates_allom_zroot_max_dbh:long_name = "dbh at which a plant reaches the maximum value for its maximum rooting depth (NOT USED)" ; + fates_allom_zroot_max_dbh:long_name = "dbh at which a plant reaches the maximum value for its maximum rooting depth" ; double fates_allom_zroot_max_z(fates_pft) ; fates_allom_zroot_max_z:units = "m" ; - fates_allom_zroot_max_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh (NOT USED). note: max_z=min_z=large, sets rooting depth to soil depth" ; + fates_allom_zroot_max_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh. note: max_z=min_z=large, sets rooting depth to soil depth" ; double fates_allom_zroot_min_dbh(fates_pft) ; fates_allom_zroot_min_dbh:units = "cm" ; - fates_allom_zroot_min_dbh:long_name = "dbh at which the maximum rooting depth for a recruit is defined (NOT USED)" ; + fates_allom_zroot_min_dbh:long_name = "dbh at which the maximum rooting depth for a recruit is defined" ; double fates_allom_zroot_min_z(fates_pft) ; fates_allom_zroot_min_z:units = "m" ; - fates_allom_zroot_min_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_min_dbh (NOT USED) note: max_z=min_z=large, sets rooting depth to soil depth" ; + fates_allom_zroot_min_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_min_dbh. note: max_z=min_z=large, sets rooting depth to soil depth" ; double fates_branch_turnover(fates_pft) ; fates_branch_turnover:units = "yr" ; fates_branch_turnover:long_name = "turnover time of branches" ; diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90 index dcf20dbd14..04a0f5dda0 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -101,7 +101,8 @@ module PRTParametersMod real(r8), allocatable :: c2b(:) ! Carbon to biomass multiplier [kg/kgC] real(r8), allocatable :: wood_density(:) ! wood density g cm^-3 ... real(r8), allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) - + real(r8), allocatable :: crown(:) ! fraction of the height of the plant + ! that is occupied by crown real(r8), allocatable :: slamax(:) ! Maximum specific leaf area of plant (at bottom) [m2/gC] real(r8), allocatable :: slatop(:) ! Specific leaf area at canopy top [m2/gC] real(r8), allocatable :: allom_sai_scaler(:) ! @@ -137,7 +138,6 @@ module PRTParametersMod real(r8), allocatable :: allom_agb3(:) ! Parameter 3 for agb allometry real(r8), allocatable :: allom_agb4(:) ! Parameter 3 for agb allometry - ! ------------------------ (NOT YET IMPLEMENTED) ------------------------- real(r8), allocatable :: allom_zroot_max_dbh(:) ! dbh at which maximum rooting depth saturates (largest possible) [cm] real(r8), allocatable :: allom_zroot_max_z(:) ! the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh [m] real(r8), allocatable :: allom_zroot_min_dbh(:) ! dbh at which the maximum rooting depth for a recruit is defined [cm] diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index dce172d47d..208ff848fb 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -178,6 +178,10 @@ subroutine PRTRegisterPFT(fates_params) name = 'fates_woody' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_fire_crown_depth_frac' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) name = 'fates_wood_density' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & @@ -430,6 +434,10 @@ subroutine PRTReceivePFT(fates_params) name = 'fates_fnrt_prof_mode' call fates_params%RetreiveParameterAllocate(name=name, & data=prt_params%fnrt_prof_mode) + + name = 'fates_fire_crown_depth_frac' + call fates_params%RetreiveParameterAllocate(name=name, & + data=prt_params%crown) name = 'fates_woody' call fates_params%RetreiveParameterAllocate(name=name, & @@ -903,6 +911,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'prt_phos_stoich_p2 = ',prt_params%phos_stoich_p2 write(fates_log(),fmt0) 'prt_alloc_priority = ',prt_params%alloc_priority write(fates_log(),fmt0) 'woody = ',prt_params%woody + write(fates_log(),fmt0) 'crown = ',prt_params%crown write(fates_log(),fmt0) 'roota_par = ',prt_params%fnrt_prof_a write(fates_log(),fmt0) 'rootb_par = ',prt_params%fnrt_prof_b write(fates_log(),fmt0) 'fnrt_prof_mode = ',prt_params%fnrt_prof_mode diff --git a/tools/BatchPatchParams.py b/tools/BatchPatchParams.py index 57edb7dfcb..ee78ebcbd0 100755 --- a/tools/BatchPatchParams.py +++ b/tools/BatchPatchParams.py @@ -51,8 +51,11 @@ def parse_syscall_str(fnamein,fnameout,param_name,param_val): sys_call_str = "../tools/modify_fates_paramfile.py"+" --fin " + fnamein + \ " --fout " + fnameout + " --var " + param_name + " --silent " +\ - " --val " + param_val + " --overwrite --all" + " --val " + "\" "+param_val+"\"" + " --overwrite --all" + + print(sys_call_str) + return(sys_call_str)