diff --git a/README b/README index 029a891482..67bc7e13c4 100644 --- a/README +++ b/README @@ -111,3 +111,4 @@ components/clm/src/biogeophys -- Biogeophysics (Hydrology) # (NOTE: ./xmlchange RESUBMIT=10 to set RESUBMIT to number # # of times to automatically resubmit -- 10 in this example) +END \ No newline at end of file diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 59ba94254f..d921e42533 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -12,6 +12,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : fusetol, nclmax use EDtypesMod , only : ncwd, numcohortsperpatch, udata + use EDtypesMod , only : sclass_ed,nlevsclass_ed,AREA ! implicit none private @@ -108,7 +109,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & endif ! Calculate live biomass allocation - call allocate_live_biomass(new_cohort) + call allocate_live_biomass(new_cohort,0) ! Assign canopy extent and depth new_cohort%c_area = c_area(new_cohort) @@ -134,6 +135,13 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & patchptr%shortest => new_cohort endif + ! Recuits do not have mortality rates, nor have they moved any + ! carbon when they are created. They will bias our statistics + ! until they have experienced a full day. We need a newly recruited flag. + ! This flag will be set to false after it has experienced + ! growth, disturbance and mortality. + new_cohort%isnew = .true. + call insert_cohort(new_cohort, patchptr%tallest, patchptr%shortest, tnull, snull, & storebigcohort, storesmallcohort) @@ -143,7 +151,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & end subroutine create_cohort !-------------------------------------------------------------------------------------! - subroutine allocate_live_biomass(cc_p) + subroutine allocate_live_biomass(cc_p,mode) ! ! !DESCRIPTION: ! Divide alive biomass between leaf, root and sapwood parts. @@ -153,6 +161,7 @@ subroutine allocate_live_biomass(cc_p) ! ! !ARGUMENTS type (ed_cohort_type), intent(inout), target :: cc_p ! current cohort pointer + integer , intent(in) :: mode ! ! !LOCAL VARIABLES: type (ed_cohort_type), pointer :: currentCohort @@ -170,8 +179,7 @@ subroutine allocate_live_biomass(cc_p) ft = currentcohort%pft leaf_frac = 1.0_r8/(1.0_r8 + EDecophyscon%sapwood_ratio(ft) * currentcohort%hite + pftcon%froot_leaf(ft)) - currentcohort%bl = currentcohort%balive*leaf_frac - ratio_balive = 1.0_r8 + !currentcohort%bl = currentcohort%balive*leaf_frac !for deciduous trees, there are no leaves if (pftcon%evergreen(ft) == 1) then @@ -179,11 +187,11 @@ subroutine allocate_live_biomass(cc_p) currentcohort%status_coh = 2 endif - !diagnore the root and stem biomass from the functional balance hypothesis. This is used when the leaves are + ! iagnore the root and stem biomass from the functional balance hypothesis. This is used when the leaves are !fully on. - currentcohort%br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac - currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & - currentcohort%laimemory)*leaf_frac + !currentcohort%br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac + !currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & + ! currentcohort%laimemory)*leaf_frac leaves_off_switch = 0 if (currentcohort%status_coh == 1.and.pftcon%stress_decid(ft) == 1.and.currentcohort%siteptr%dstatus==1) then !no leaves @@ -193,12 +201,46 @@ subroutine allocate_live_biomass(cc_p) leaves_off_switch = 1 !cold decid endif - if (leaves_off_switch==1) then - !the purpose of this section is to figure out the root and stem biomass when the leaves are off - !at this point, we know the former leaf mass (laimemory) and the current alive mass - !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted - !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass - !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF + ! Use different proportions if the leaves are on vs off + if(leaves_off_switch==0)then + + ! Tracking npp/gpp diagnostics only occur after growth derivatives is called + if(mode==1)then + ! it will not be able to put out as many leaves as it had previous timestep + currentcohort%npp_leaf = currentcohort%npp_leaf + & + max(0.0_r8,currentcohort%balive*leaf_frac - currentcohort%bl)/udata%deltat + end if + + currentcohort%bl = currentcohort%balive*leaf_frac + + !diagnose the root and stem biomass from the functional balance hypothesis. This is used when the leaves are + !fully on. + if(mode==1)then + + currentcohort%npp_froot = currentcohort%npp_froot + & + max(0._r8,pftcon%froot_leaf(ft)*(currentcohort%balive+currentcohort%laimemory)*leaf_frac - currentcohort%br)/udata%deltat + + currentcohort%npp_bsw = max(0._r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & + currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/udata%deltat + + currentcohort%npp_bdead = currentCohort%dbdeaddt + + end if + + currentcohort%br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac + currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & + currentcohort%laimemory)*leaf_frac + + + else ! Leaves are on (leaves_off_switch==1) + + !the purpose of this section is to figure out the root and stem biomass when the leaves are off + !at this point, we know the former leaf mass (laimemory) and the current alive mass + !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted + !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass + !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF + + currentcohort%bl = 0.0_r8 ideal_balive = currentcohort%laimemory * pftcon%froot_leaf(ft) + & currentcohort%laimemory* EDecophyscon%sapwood_ratio(ft) * currentcohort%hite @@ -208,7 +250,21 @@ subroutine allocate_live_biomass(cc_p) ratio_balive = currentcohort%balive / ideal_balive currentcohort%br = currentcohort%br * ratio_balive - currentcohort%bsw = currentcohort%bsw * ratio_balive + currentcohort%bsw = currentcohort%bsw * ratio_balive + + ! Diagnostics + if(mode==1)then + + currentcohort%npp_froot = currentcohort%npp_froot + & + max(0.0_r8,pftcon%froot_leaf(ft)*(ideal_balive + & + currentcohort%laimemory)*leaf_frac*ratio_balive-currentcohort%br)/udata%deltat + + currentcohort%npp_bsw = max(0.0_r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + & + currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/udata%deltat + + currentcohort%npp_bdead = currentCohort%dbdeaddt + + end if endif @@ -294,6 +350,14 @@ subroutine nan_cohort(cc_p) currentCohort%resp_clm = nan ! RESP: kgC/indiv/timestep currentCohort%resp_acc = nan ! RESP: kGC/cohort/day + currentCohort%npp_leaf = nan + currentCohort%npp_froot = nan + currentCohort%npp_bsw = nan + currentCohort%npp_bdead = nan + currentCohort%npp_bseed = nan + currentCohort%npp_store = nan + + !RESPIRATION currentCohort%rd = nan currentCohort%resp_m = nan ! Maintenance respiration. kGC/cohort/year @@ -387,6 +451,13 @@ subroutine zero_cohort(cc_p) currentcohort%gscan = 0._r8 currentcohort%treesai = 0._r8 + ! currentCohort%npp_leaf = 0._r8 + ! currentCohort%npp_froot = 0._r8 + ! currentCohort%npp_bsw = 0._r8 + ! currentCohort%npp_bdead = 0._r8 + ! currentCohort%npp_bseed = 0._r8 + ! currentCohort%npp_store = 0._r8 + end subroutine zero_cohort !-------------------------------------------------------------------------------------! @@ -533,7 +604,9 @@ subroutine fuse_cohorts(patchptr) iterate = 1 fusion_took_place = 0 currentPatch => patchptr - maxcohorts = currentPatch%NCL_p * numCohortsPerPatch + ! maxcohorts = currentPatch%NCL_p * numCohortsPerPatch + maxcohorts = numCohortsPerPatch + !---------------------------------------------------------------------! ! Keep doing this until nocohorts <= maxcohorts ! !---------------------------------------------------------------------! @@ -559,7 +632,17 @@ subroutine fuse_cohorts(patchptr) if (currentCohort%pft == nextc%pft) then ! check cohorts in same c. layer. before fusing - if (currentCohort%canopy_layer == nextc%canopy_layer) then + + if (currentCohort%canopy_layer == nextc%canopy_layer) then + + ! Note: because newly recruited cohorts that have not experienced + ! a day yet will have un-known flux quantities or change rates + ! we don't want them fusing with non-new cohorts. We allow them + ! to fuse with other new cohorts to keep the total number of cohorts + ! down. + + if( (.not.(currentCohort%isnew) .and. .not.(nextc%isnew) ) .or. & + ( currentCohort%isnew .and. nextc%isnew ) ) then fusion_took_place = 1 newn = currentCohort%n + nextc%n ! sum individuals in both cohorts. @@ -613,10 +696,25 @@ subroutine fuse_cohorts(patchptr) currentCohort%npp = (currentCohort%n*currentCohort%npp + nextc%n*nextc%npp)/newn currentCohort%gpp = (currentCohort%n*currentCohort%gpp + nextc%n*nextc%gpp)/newn currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim + nextc%n*nextc%canopy_trim)/newn - currentCohort%dmort = (currentCohort%n*currentCohort%dmort + nextc%n*nextc%dmort)/newn + currentCohort%dmort = (currentCohort%n*currentCohort%dmort + nextc%n*nextc%dmort)/newn currentCohort%fire_mort = (currentCohort%n*currentCohort%fire_mort + nextc%n*nextc%fire_mort)/newn currentCohort%leaf_litter = (currentCohort%n*currentCohort%leaf_litter + nextc%n*nextc%leaf_litter)/newn + ! mortality diagnostics + currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn + currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn + currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn + currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn + currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn + + ! npp diagnostics + currentCohort%npp_leaf = (currentCohort%n*currentCohort%npp_leaf + nextc%n*nextc%npp_leaf)/newn + currentCohort%npp_froot = (currentCohort%n*currentCohort%npp_froot + nextc%n*nextc%npp_froot)/newn + currentCohort%npp_bsw = (currentCohort%n*currentCohort%npp_bsw + nextc%n*nextc%npp_bsw)/newn + currentCohort%npp_bdead = (currentCohort%n*currentCohort%npp_bdead + nextc%n*nextc%npp_bdead)/newn + currentCohort%npp_bseed = (currentCohort%n*currentCohort%npp_bseed + nextc%n*nextc%npp_bseed)/newn + currentCohort%npp_store = (currentCohort%n*currentCohort%npp_store + nextc%n*nextc%npp_store)/newn + do i=1, nlevcan_ed if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then currentCohort%year_net_uptake(i) = min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i)) @@ -639,6 +737,8 @@ subroutine fuse_cohorts(patchptr) deallocate(nextc) endif + endif ! Not a recruit + endif !canopy layer endif !pft endif !index no. @@ -917,6 +1017,13 @@ subroutine copy_cohort( currentCohort,copyc ) n%year_net_uptake = o%year_net_uptake n%ts_net_uptake = o%ts_net_uptake + n%npp_leaf = o%npp_leaf + n%npp_froot = o%npp_froot + n%npp_bsw = o%npp_bsw + n%npp_bdead = o%npp_bdead + n%npp_bseed = o%npp_bseed + n%npp_store = o%npp_store + !RESPIRATION n%rd = o%rd n%resp_m = o%resp_m @@ -943,6 +1050,16 @@ subroutine copy_cohort( currentCohort,copyc ) n%c_area = o%c_area n%woody_turnover = o%woody_turnover + ! Mortality diagnostics + n%cmort = o%cmort + n%bmort = o%bmort + n%imort = o%imort + n%fmort = o%fmort + n%hmort = o%hmort + + ! Flags + n%isnew = o%isnew + ! VARIABLES NEEDED FOR INTEGRATION n%dndt = o%dndt n%dhdt = o%dhdt diff --git a/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 b/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 index aa23223d3e..98925f1684 100755 --- a/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 @@ -322,7 +322,11 @@ real(r8) function dDbhdBl( cohort_in ) dblddbh = 1.56_r8*0.0419_r8*(cohort_in%dbh**0.56_r8)*(EDecophyscon%wood_density(cohort_in%pft)**0.55_r8) dblddbh = dblddbh*cohort_in%canopy_trim - dDbhdBl = 1.0_r8/dblddbh + if( cohort_in%dbh 0._r8 ) then - if(Bleaf(cohort_in) > 0._r8.and.cohort_in%bstore <= Bleaf(cohort_in))then + if(Bleaf(cohort_in) > 0._r8 .and. cohort_in%bstore <= Bleaf(cohort_in))then frac = cohort_in%bstore/(Bleaf(cohort_in)) - smort = smort + max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac)) + cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac)) + else + cmort = 0.0_r8 endif + else write(iulog,*) 'dbh problem in mortality_rates', & cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer,cohort_in%indexnumber endif - mortality_rates = smort + bmort + !mortality_rates = bmort + hmort + cmort - end function mortality_rates + end subroutine mortality_rates ! ============================================================================ diff --git a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 index da8f6f8cca..780beee01f 100755 --- a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 @@ -5,10 +5,11 @@ module EDPatchDynamicsMod ! ============================================================================ use shr_kind_mod , only : r8 => shr_kind_r8; + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use clm_varctl , only : iulog use pftconMod , only : pftcon use EDCohortDynamicsMod , only : fuse_cohorts, sort_cohorts, insert_cohort - use EDtypesMod , only : ncwd, n_dbh_bins, ntol, numpft_ed, area, dbhmax + use EDtypesMod , only : ncwd, n_dbh_bins, ntol, numpft_ed, area, dbhmax, numPatchesPerGridCell use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, udata ! implicit none @@ -50,6 +51,9 @@ subroutine disturbance_rates( site_in) ! !LOCAL VARIABLES: type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type), pointer :: currentCohort + real(r8) :: cmort + real(r8) :: bmort + real(r8) :: hmort !--------------------------------------------------------------------- !MORTALITY @@ -65,9 +69,17 @@ subroutine disturbance_rates( site_in) ! Mortality for trees in the understorey. currentCohort%patchptr => currentPatch - currentCohort%dmort = mortality_rates(currentCohort) + call mortality_rates(currentCohort,cmort,hmort,bmort) + currentCohort%dmort = cmort+hmort+bmort currentCohort%c_area = c_area(currentCohort) + ! Initialize diagnostic mortality rates + currentCohort%cmort = cmort + currentCohort%bmort = bmort + currentCohort%hmort = hmort + currentCohort%imort = 0.0_r8 ! Impact mortality is always zero except in new patches + currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed + if(currentCohort%canopy_layer == 1)then currentPatch%disturbance_rates(1) = currentPatch%disturbance_rates(1) + & @@ -91,6 +103,27 @@ subroutine disturbance_rates( site_in) !Only use larger of two natural disturbance modes WHY? if(currentPatch%disturbance_rates(2) > currentPatch%disturbance_rates(1))then ! DISTURBANCE IS FIRE currentPatch%disturbance_rate = currentPatch%disturbance_rates(2) + + ! RGK 02-18-2014 + ! Since treefall mortality is not actually being applied + ! Go through and zero the diagnostic rates + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer == 1)then + currentCohort%cmort=0.0_r8 + currentCohort%hmort=0.0_r8 + currentCohort%bmort=0.0_r8 + end if + + ! This may be counter-intuitive, but the diagnostic fire-mortality rate + ! will stay zero in the patch that undergoes fire, this is because + ! the actual cohorts who experience the fire are only those in the + ! newly created patch so currentCohort%fmort = 0.0_r8 + ! Don't worry, the cohorts in the newly created patch will reflect burn + + currentCohort => currentCohort%taller + enddo !currentCohort + else currentPatch%disturbance_rate = currentPatch%disturbance_rates(1) ! DISTURBANCE IS MORTALITY endif @@ -202,7 +235,7 @@ subroutine spawn_patches( currentSite ) do while(associated(currentPatch)) patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate ! how much land is disturbed in this donor patch? - call average_patch_properties(currentPatch, new_patch, patch_site_areadis) + call average_patch_properties(currentPatch, new_patch, patch_site_areadis) ! MAY BE REDUNDANT CALL if (currentSite%disturbance_mortality > currentSite%disturbance_fire) then !mortality is dominant disturbance call mortality_litter_fluxes(currentPatch, new_patch, patch_site_areadis) else @@ -226,35 +259,83 @@ subroutine spawn_patches( currentSite ) !mortality is dominant disturbance if(currentPatch%disturbance_rates(1) > currentPatch%disturbance_rates(2))then - if(currentCohort%canopy_layer == 1)then - ! keep the trees that didn't die + if(currentCohort%canopy_layer == 1)then + ! In the donor patch we are left with fewer trees because the area has decreased + ! the plant density for large trees does not actually decrease in the donor patch + ! because this is the part of the original patch where no trees have actually fallen + ! The diagnostic cmort,bmort and hmort rates have already been saved + currentCohort%n = currentCohort%n * (1.0_r8 - min(1.0_r8,currentCohort%dmort * udata%deltat)) - nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. - else + + nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + nc%cmort = nan ! The mortality diagnostics are set to nan because the cohort should dissappear + nc%hmort = nan + nc%bmort = nan + nc%fmort = nan + nc%imort = nan + else + ! small trees if(pftcon%woody(currentCohort%pft) == 1)then + ! Number of trees in the understory of new patch, before we impose impact mortality and survivorship + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + ! remaining of understory plants of those that are knocked over by the overstorey trees dying... nc%n = (1.0_r8 - ED_val_understorey_death) * currentCohort%n * patch_site_areadis/currentPatch%area + + ! since the donor patch split and sent a fraction of its members + ! to the new patch and a fraction to be preserved in itself, + ! when reporting diagnostic rates, we must carry over the mortality rates from + ! the donor that were applied before the patch split. Remember this is only + ! for diagnostics. But think of it this way, the rates are weighted by + ! number density in EDCLMLink, and the number density of this new patch is donated + ! so with the number density must come the effective mortality rates. + + nc%fmort = 0.0_r8 ! Should had also been zero in the donor + nc%imort = ED_val_understorey_death/udata%deltat ! This was zero in the donor + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + ! understory trees that might potentially be knocked over in the disturbance. - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! The existing (donor) patch should not have any impact mortality, it should + ! only lose cohorts due to the decrease in area. This is not mortality. + ! Besides, the current and newly created patch sum to unity + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) else + ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - ! remaining of understory plants of those that are knocked over by the overstorey trees dying... - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - ! understory trees that might potentially be knocked over in the disturbance. + ! Those remaining in the existing currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + nc%fmort = nan ! These should not make it to the diagnostics + nc%imort = nan ! If they do.. they should invalidate it + nc%cmort = nan ! + nc%hmort = nan ! + nc%bmort = nan ! + endif endif else !fire - ! loss of individual from fire in new patch. - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area * (1.0_r8 - currentCohort%fire_mort) - ! loss of individuals from source patch + ! Number of members in the new patch, before we impose fire survivorship + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! loss of individuals from source patch due to area shrinking currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + ! loss of individual from fire in new patch. + nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) + + nc%fmort = currentCohort%fire_mort/udata%deltat + nc%imort = 0.0_r8 + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + endif if (nc%n > 0.0_r8) then @@ -919,7 +1000,8 @@ subroutine fuse_patches( csite ) integer :: fuse_flag !do patches get fused (1) or not (0). !--------------------------------------------------------------------- - maxpatch = 4 + !maxpatch = 4 + maxpatch = numPatchesPerGridCell currentSite => csite diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index c042e5e3d3..b949c5fd76 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -732,6 +732,9 @@ subroutine Growth_Derivatives( currentCohort) real(r8) :: f_store !fraction of NPP allocated to storage in this timestep (functionf of stored pool) real(r8) :: gr_fract !fraction of carbon balance that is allocated to growth (not reproduction) real(r8) :: target_balive !target leaf biomass under allometric optimum. + real(r8) :: cmort ! starvation mortality rate (fraction per year) + real(r8) :: bmort ! background mortality rate (fraction per year) + real(r8) :: hmort ! hydraulic failure mortality rate (fraction per year) real(r8) :: balive_loss !---------------------------------------------------------------------- @@ -740,7 +743,8 @@ subroutine Growth_Derivatives( currentCohort) ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology if (currentCohort%canopy_layer > 1)then - currentCohort%dndt = -1.0_r8 * mortality_rates(currentCohort) * currentCohort%n + call mortality_rates(currentCohort,cmort,hmort,bmort) + currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort) * currentCohort%n else currentCohort%dndt = 0._r8 endif @@ -749,7 +753,7 @@ subroutine Growth_Derivatives( currentCohort) currentCohort%hite = Hite(currentCohort) h = currentCohort%hite - call allocate_live_biomass(currentCohort) + call allocate_live_biomass(currentCohort,0) ! calculate target size of living biomass compartment for a given dbh. target_balive = Bleaf(currentCohort) * (1.0_r8 + pftcon%froot_leaf(currentCohort%pft) + & @@ -812,6 +816,14 @@ subroutine Growth_Derivatives( currentCohort) currentCohort%carbon_balance = currentCohort%npp - currentCohort%md * EDecophyscon%leaf_stor_priority(currentCohort%pft) + ! Allowing only carbon from NPP pool to account for npp flux into the maintenance pools + ! ie this does not include any use of storage carbon or balive to make up for missing carbon balance in the transfer + currentCohort%npp_leaf = min(currentCohort%npp*currentCohort%leaf_md/currentCohort%md, & + currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)) + currentCohort%npp_froot = min(currentCohort%npp*currentCohort%root_md/currentCohort%md, & + currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)) + + if (Bleaf(currentCohort) > 0._r8)then if ( DEBUG ) write(iulog,*) 'EDphys A ',currentCohort%carbon_balance @@ -849,7 +861,19 @@ subroutine Growth_Derivatives( currentCohort) if (currentCohort%carbon_balance > currentCohort%md*(1.0_r8- EDecophyscon%leaf_stor_priority(currentCohort%pft)))then ! Yes... currentCohort%carbon_balance = currentCohort%carbon_balance - currentCohort%md * (1.0_r8 - & EDecophyscon%leaf_stor_priority(currentCohort%pft)) + + currentCohort%npp_leaf = currentCohort%npp_leaf + & + currentCohort%leaf_md * (1.0_r8-EDecophyscon%leaf_stor_priority(currentCohort%pft)) + currentCohort%npp_froot = currentCohort%npp_froot + & + currentCohort%root_md * (1.0_r8-EDecophyscon%leaf_stor_priority(currentCohort%pft)) + else ! we can't maintain constant leaf area and root area. Balive is reduced + + currentCohort%npp_leaf = currentCohort%npp_leaf + & + max(0.0_r8,currentCohort%carbon_balance*(currentCohort%leaf_md/currentCohort%md)) + currentCohort%npp_froot = currentCohort%npp_froot + & + max(0.0_r8,currentCohort%carbon_balance*(currentCohort%root_md/currentCohort%md)) + balive_loss = currentCohort%md *(1.0_r8- EDecophyscon%leaf_stor_priority(currentCohort%pft))- currentCohort%carbon_balance currentCohort%carbon_balance = 0._r8 endif @@ -921,10 +945,16 @@ subroutine Growth_Derivatives( currentCohort) currentCohort%dbalivedt = 0._r8 endif + currentCohort%npp_bseed = currentCohort%seed_prod + currentCohort%npp_store = max(0.0_r8,currentCohort%storage_flux) + ! calculate change in diameter and height currentCohort%ddbhdt = currentCohort%dbdeaddt * dDbhdBd(currentCohort) currentCohort%dhdt = currentCohort%dbdeaddt * dHdBd(currentCohort) + ! If the cohort has grown, it is not new + currentCohort%isnew=.false. + end subroutine Growth_Derivatives ! ============================================================================ diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index fd369fd1ef..6bf500cc52 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -438,6 +438,7 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, & !RF - copied this from the CLM trunk code, but where did it come from, and how can we make these consistant? !jmax25top(FT) = (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrz),11._r8),35._r8)) * vcmax25top(FT) + jmax25top(FT) = 1.67_r8 * vcmax25top(FT) tpu25top(FT) = 0.167_r8 * vcmax25top(FT) kp25top(FT) = 20000._r8 * vcmax25top(FT) diff --git a/components/clm/src/ED/main/EDCLMLinkMod.F90 b/components/clm/src/ED/main/EDCLMLinkMod.F90 index 054a34dbb9..5257751e2d 100755 --- a/components/clm/src/ED/main/EDCLMLinkMod.F90 +++ b/components/clm/src/ED/main/EDCLMLinkMod.F90 @@ -5,11 +5,15 @@ module EDCLMLinkMod ! diagnostics, or as input to the land surface components. ! ============================================================================ - use shr_kind_mod , only : r8 => shr_kind_r8; + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_infnan_mod, only : isnan => shr_infnan_isnan use decompMod , only : bounds_type - use clm_varpar , only : nclmax, nlevcan_ed, numpft, numcft + use clm_varpar , only : nclmax, nlevcan_ed, numpft, numcft, mxpft use clm_varctl , only : iulog - use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type + use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, nlevsclass_ed + use EDtypesMod , only : sclass_ed,nlevsclass_ed,AREA + use EDParamsMod , only : ED_val_ag_biomass + ! implicit none private @@ -69,7 +73,38 @@ module EDCLMLinkMod real(r8), pointer, private :: deadstemc_patch (:) ! (gC/m2) dead stem C real(r8), pointer, private :: livestemn_patch (:) ! (gN/m2) live stem N real(r8), pointer, private :: npp_patch (:) ! (gC/m2/s) patch net primary production - real(r8), pointer, public :: gpp_patch (:) ! (gC/m2/s) patch gross primary production + + real(r8), pointer, private :: gpp_patch (:) ! (gC/m2/s) patch gross primary production + + real(r8), pointer :: ed_gpp_gd_scpf (:,:) ! [kg/m2/yr] gross primary production + real(r8), pointer :: ed_npp_totl_gd_scpf (:,:) ! [kg/m2/yr] net primary production (npp) + real(r8), pointer :: ed_npp_leaf_gd_scpf (:,:) ! [kg/m2/yr] npp flux into leaf pool + real(r8), pointer :: ed_npp_seed_gd_scpf (:,:) ! [kg/m2/yr] npp flux into flower,fruit,nut,seed + real(r8), pointer :: ed_npp_fnrt_gd_scpf (:,:) ! [kg/m2/yr] npp flux into fine roots + real(r8), pointer :: ed_npp_bgsw_gd_scpf (:,:) ! [kg/m2/yr] npp flux into below ground sapwood + real(r8), pointer :: ed_npp_bgdw_gd_scpf (:,:) ! [kg/m2/yr] npp flux into below ground structural (dead) wood + real(r8), pointer :: ed_npp_agsw_gd_scpf (:,:) ! [kg/m2/yr] npp flux into above ground sapwood + real(r8), pointer :: ed_npp_agdw_gd_scpf (:,:) ! [kg/m2/yr] npp flux into below ground structural (dead) wood + real(r8), pointer :: ed_npp_stor_gd_scpf (:,:) ! [kg/m2/yr] npp flux through the storage pool + real(r8), pointer :: ed_litt_leaf_gd_scpf (:,:) ! [kg/m2/yr] carbon flux of live leaves to litter + real(r8), pointer :: ed_litt_fnrt_gd_scpf (:,:) ! [kg/m2/yr] carbon flux of fine roots to litter + real(r8), pointer :: ed_litt_sawd_gd_scpf (:,:) ! [kg/m2/yr] carbon flux of sapwood to litter (above+below) + real(r8), pointer :: ed_litt_ddwd_gd_scpf (:,:) ! [kg/m2/yr] carbon flux of dead wood (above+below) to litter + real(r8), pointer :: ed_r_leaf_gd_scpf (:,:) ! [kg/m2/yr] total leaf respiration + real(r8), pointer :: ed_r_stem_gd_scpf (:,:) ! [kg/m2/yr] total above ground live wood (stem) respiration + real(r8), pointer :: ed_r_root_gd_scpf (:,:) ! [kg/m2/yr] total below ground live wood (root) respiration + real(r8), pointer :: ed_r_stor_gd_scpf (:,:) ! [kg/m2/yr] total storage respiration + + ! Carbon State Variables for direct comparison to inventory - dimensions: (disturbance patch, pft x size) + + real(r8), pointer :: ed_ddbh_gd_scpf (:,:) ! [cm/yr] diameter increment + real(r8), pointer :: ed_ba_gd_scpf (:,:) ! [m2/ha] basal area + real(r8), pointer :: ed_np_gd_scpf (:,:) ! [/m2] number of plants + real(r8), pointer :: ed_m1_gd_scpf (:,:) ! [Stems/ha/yr] Mean Background Mortality + real(r8), pointer :: ed_m2_gd_scpf (:,:) ! [Stems/ha/yr] Mean Hydraulic Mortaliry + real(r8), pointer :: ed_m3_gd_scpf (:,:) ! [Stems/ha/yr] Mean Carbon Starvation Mortality + real(r8), pointer :: ed_m4_gd_scpf (:,:) ! [Stems/ha/yr] Mean Impact Mortality + real(r8), pointer :: ed_m5_gd_scpf (:,:) ! [Stems/ha/yr] Mean Fire Mortality contains @@ -123,9 +158,13 @@ subroutine InitAllocate(this, bounds) ! ! !LOCAL VARIABLES: integer :: begp,endp + integer :: begc,endc !bounds + integer :: begg,endg !------------------------------------------------------------------------ begp = bounds%begp; endp = bounds%endp + begc = bounds%begc; endc = bounds%endc + allocate(this%trimming_patch (begp:endp)) ; this%trimming_patch (:) = 0.0_r8 allocate(this%canopy_spread_patch (begp:endp)) ; this%canopy_spread_patch (:) = 0.0_r8 @@ -171,6 +210,36 @@ subroutine InitAllocate(this, bounds) allocate(this%gpp_patch (begp:endp)) ; this%gpp_patch (:) = nan allocate(this%npp_patch (begp:endp)) ; this%npp_patch (:) = nan + begg = bounds%begg; endg = bounds%endg + allocate(this%ed_gpp_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_gpp_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_totl_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_totl_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_leaf_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_leaf_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_seed_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_seed_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_fnrt_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_fnrt_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_bgsw_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_bgsw_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_bgdw_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_bgdw_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_agsw_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_agsw_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_agdw_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_agdw_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_npp_stor_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_stor_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_litt_leaf_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_litt_leaf_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_litt_fnrt_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_litt_fnrt_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_litt_sawd_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_litt_sawd_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_litt_ddwd_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_litt_ddwd_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_r_leaf_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_r_leaf_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_r_stem_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_r_stem_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_r_root_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_r_root_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_r_stor_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_r_stor_gd_scpf (:,:) = 0.0_r8 + + ! Carbon State Variables for direct comparison to inventory - dimensions: (disturbance patch, pft x size) + allocate(this%ed_ddbh_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_ddbh_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_ba_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_ba_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_np_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_np_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_m1_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_m1_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_m2_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_m2_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_m3_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_m3_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_m4_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_m4_gd_scpf (:,:) = 0.0_r8 + allocate(this%ed_m5_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)) ; this%ed_m5_gd_scpf (:,:) = 0.0_r8 + end subroutine InitAllocate !------------------------------------------------------------------------ @@ -379,6 +448,78 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='net primary production', & ptr_patch=this%npp_patch) + + ! Carbon Flux (grid dimension x scpf) + ! ============================================================== + + call hist_addfld2d (fname='ED_GPP_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='gross primary production', & + ptr_gcell=this%ed_gpp_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_LEAF_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into leaves', & + ptr_gcell=this%ed_npp_leaf_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_SEED_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into seeds', & + ptr_gcell=this%ed_npp_seed_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_FNRT_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into fine roots', & + ptr_gcell=this%ed_npp_fnrt_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_BGSW_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into below-ground sapwood', & + ptr_gcell=this%ed_npp_bgsw_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_BGDW_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into below-ground deadwood', & + ptr_gcell=this%ed_npp_bgdw_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_NPP_AGSW_GD_SCPF',units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into above-ground sapwood', & + ptr_gcell=this%ed_npp_agsw_gd_scpf,default='inactive') + + call hist_addfld2d ( fname = 'ED_NPP_AGDW_GD_SCPF', units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into above-ground deadwood', & + ptr_gcell=this%ed_npp_agdw_gd_scpf,default='inactive') + + call hist_addfld2d ( fname = 'ED_NPP_STOR_GD_SCPF', units='kgC/m2/yr',type2d='levscpf',& + avgflag='A', long_name='NPP flux into storage', & + ptr_gcell=this%ed_npp_stor_gd_scpf,default='inactive') + + call hist_addfld2d (fname='ED_DDBH_GD_SCPF', units = 'cm/yr/ha', type2d = 'levscpf', & + avgflag='A', long_name='diameter growth increment and pft/size', & + ptr_gcell=this%ed_ddbh_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_BA_GD_SCPF',units = 'm2/ha', type2d = 'levscpf', & + avgflag='A', long_name='basal area by patch and pft/size', & + ptr_gcell=this%ed_ba_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_NPLANT_GD_SCPF',units = 'N/ha', type2d = 'levscpf', & + avgflag='A', long_name='stem number density by patch and pft/size', & + ptr_gcell=this%ed_np_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_M1_GD_SCPF',units = 'N/ha/yr', type2d = 'levscpf', & + avgflag='A', long_name='background mortality count by patch and pft/size', & + ptr_gcell=this%ed_m1_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_M2_GD_SCPF',units = 'N/ha/yr', type2d = 'levscpf', & + avgflag='A', long_name='hydraulic mortality count by patch and pft/size', & + ptr_gcell=this%ed_m2_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_M3_GD_SCPF',units = 'N/ha/yr', type2d = 'levscpf', & + avgflag='A', long_name='carbon starvation mortality count by patch and pft/size', & + ptr_gcell=this%ed_m3_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_M4_GD_SCPF',units = 'N/ha/yr', type2d = 'levscpf', & + avgflag='A', long_name='impact mortality count by patch and pft/size', & + ptr_gcell=this%ed_m4_gd_scpf, default='inactive') + + call hist_addfld2d (fname='ED_M5_GD_SCPF',units = 'N/ha/yr', type2d = 'levscpf', & + avgflag='A', long_name='fire mortality count by patch and pft/size', & + ptr_gcell=this%ed_m5_gd_scpf, default='inactive') + end subroutine InitHistory !----------------------------------------------------------------------- @@ -744,6 +885,7 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & use EDPhenologyType , only : ed_phenology_type use CanopyStateType , only : canopystate_type use PatchType , only : clmpatch => patch + use pftconMod , only : pftcon ! ! !ARGUMENTS: class(ed_clm_type) :: this @@ -758,6 +900,10 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & integer :: G,p,ft integer :: firstsoilpatch(bounds%begg:bounds%endg) real(r8) :: n_density ! individual of cohort per m2. + real(r8) :: n_perm2 ! individuals per m2 for the whole grid cell + real(r8) :: dbh ! actual dbh used to identify relevant size class + integer :: scpf ! size class x pft index + integer :: sc !----------------------------------------------------------------------- associate( & @@ -798,7 +944,27 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & gpp => this%gpp_patch , & ! Output: npp => this%npp_patch , & ! Output: - + + ed_gpp_scpf => this%ed_gpp_gd_scpf , & + ed_npp_totl_scpf => this%ed_npp_totl_gd_scpf , & + ed_npp_leaf_scpf => this%ed_npp_leaf_gd_scpf , & + ed_npp_seed_scpf => this%ed_npp_seed_gd_scpf , & + ed_npp_fnrt_scpf => this%ed_npp_fnrt_gd_scpf , & + ed_npp_bgsw_scpf => this%ed_npp_bgsw_gd_scpf , & + ed_npp_bgdw_scpf => this%ed_npp_bgdw_gd_scpf , & + ed_npp_agsw_scpf => this%ed_npp_agsw_gd_scpf , & + ed_npp_agdw_scpf => this%ed_npp_agdw_gd_scpf , & + ed_npp_stor_scpf => this%ed_npp_stor_gd_scpf , & + + ed_ddbh_gd_scpf => this%ed_ddbh_gd_scpf , & + ed_ba_gd_scpf => this%ed_ba_gd_scpf , & + ed_np_gd_scpf => this%ed_np_gd_scpf , & + ed_m1_gd_scpf => this%ed_m1_gd_scpf , & + ed_m2_gd_scpf => this%ed_m2_gd_scpf , & + ed_m3_gd_scpf => this%ed_m3_gd_scpf , & + ed_m4_gd_scpf => this%ed_m4_gd_scpf , & + ed_m5_gd_scpf => this%ed_m5_gd_scpf , & + tlai => canopystate_inst%tlai_patch , & ! InOut: elai => canopystate_inst%elai_patch , & ! InOut: tsai => canopystate_inst%tsai_patch , & ! InOut: @@ -846,6 +1012,26 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & ED_balive(:) = 0.0_r8 phen_cd_status(:) = 2 + ed_gpp_scpf(:,:) = 0.0_r8 + ed_npp_totl_scpf(:,:) = 0.0_r8 + ed_npp_leaf_scpf(:,:) = 0.0_r8 + ed_npp_seed_scpf(:,:) = 0.0_r8 + ed_npp_fnrt_scpf(:,:) = 0.0_r8 + ed_npp_bgsw_scpf(:,:) = 0.0_r8 + ed_npp_bgdw_scpf(:,:) = 0.0_r8 + ed_npp_agsw_scpf(:,:) = 0.0_r8 + ed_npp_agdw_scpf(:,:) = 0.0_r8 + ed_npp_stor_scpf(:,:) = 0.0_r8 + + ed_ddbh_gd_scpf(:,:) = 0.0_r8 + ed_ba_gd_scpf(:,:) = 0.0_r8 + ed_np_gd_scpf(:,:) = 0.0_r8 + ed_m1_gd_scpf(:,:) = 0.0_r8 + ed_m2_gd_scpf(:,:) = 0.0_r8 + ed_m3_gd_scpf(:,:) = 0.0_r8 + ed_m4_gd_scpf(:,:) = 0.0_r8 + ed_m5_gd_scpf(:,:) = 0.0_r8 + do g = bounds%begg,bounds%endg if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then @@ -907,8 +1093,10 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & ft = currentCohort%pft if(currentPatch%area>0._r8)then n_density = currentCohort%n/currentPatch%area + n_perm2 = currentCohort%n/AREA ! plant density using whole area (for grid cell averages) else n_density = 0.0_r8 + n_perm2 = 0.0_r8 endif if ( DEBUG ) then @@ -927,7 +1115,63 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & PFTleafbiomass(p,ft) = PFTleafbiomass(p,ft) + n_density * currentCohort%bl PFTstorebiomass(p,ft) = PFTstorebiomass(p,ft) + n_density * currentCohort%bstore PFTnindivs(p,ft) = PFTnindivs(p,ft) + currentCohort%n - currentCohort => currentCohort%taller + + dbh = currentCohort%dbh !-0.5*(1./365.25)*currentCohort%ddbhdt + sc = count(dbh-sclass_ed.ge.0.0) + scpf = (ft-1)*nlevsclass_ed+sc + + ! Flux Variables (must pass a NaN check on growth increment and not be recruits) + if( .not.(isnan(currentCohort%ddbhdt)) .and. .not.(currentCohort%isnew)) then + ed_gpp_scpf(g,scpf) = ed_gpp_scpf(g,scpf) + n_perm2*currentCohort%gpp ! [kgC/m2/yr] + ed_npp_totl_scpf(g,scpf) = ed_npp_totl_scpf(g,scpf) + currentcohort%npp*n_perm2 + ed_npp_leaf_scpf(g,scpf) = ed_npp_leaf_scpf(g,scpf) + currentcohort%npp_leaf*n_perm2 + ed_npp_fnrt_scpf(g,scpf) = ed_npp_fnrt_scpf(g,scpf) + currentcohort%npp_froot*n_perm2 + ed_npp_bgsw_scpf(g,scpf) = ed_npp_bgsw_scpf(g,scpf) + currentcohort%npp_bsw*(1._r8-ED_val_ag_biomass)*n_perm2 + ed_npp_agsw_scpf(g,scpf) = ed_npp_agsw_scpf(g,scpf) + currentcohort%npp_bsw*ED_val_ag_biomass*n_perm2 + ed_npp_bgdw_scpf(g,scpf) = ed_npp_bgdw_scpf(g,scpf) + currentcohort%npp_bdead*(1._r8-ED_val_ag_biomass)*n_perm2 + ed_npp_agdw_scpf(g,scpf) = ed_npp_agdw_scpf(g,scpf) + currentcohort%npp_bdead*ED_val_ag_biomass*n_perm2 + ed_npp_seed_scpf(g,scpf) = ed_npp_seed_scpf(g,scpf) + currentcohort%npp_bseed*n_perm2 + ed_npp_stor_scpf(g,scpf) = ed_npp_stor_scpf(g,scpf) + currentcohort%npp_store*n_perm2 + if( abs(currentcohort%npp-(currentcohort%npp_leaf+currentcohort%npp_froot+ & + currentcohort%npp_bsw+currentcohort%npp_bdead+ & + currentcohort%npp_bseed+currentcohort%npp_store))>1.e-9) then + write(iulog,*) 'NPP Partitions are not balancing' + write(iulog,*) 'Fractional Error: ',abs(currentcohort%npp-(currentcohort%npp_leaf+currentcohort%npp_froot+ & + currentcohort%npp_bsw+currentcohort%npp_bdead+ & + currentcohort%npp_bseed+currentcohort%npp_store))/currentcohort%npp + write(iulog,*) 'Terms: ',currentcohort%npp,currentcohort%npp_leaf,currentcohort%npp_froot, & + currentcohort%npp_bsw,currentcohort%npp_bdead, & + currentcohort%npp_bseed,currentcohort%npp_store + stop + end if + ! Woody State Variables (basal area and number density and mortality) + if (pftcon%woody(ft) == 1) then + + ed_m1_gd_scpf(g,scpf) = ed_m1_gd_scpf(g,scpf) + currentcohort%bmort*n_perm2*AREA + ed_m2_gd_scpf(g,scpf) = ed_m2_gd_scpf(g,scpf) + currentcohort%hmort*n_perm2*AREA + ed_m3_gd_scpf(g,scpf) = ed_m3_gd_scpf(g,scpf) + currentcohort%cmort*n_perm2*AREA + ed_m4_gd_scpf(g,scpf) = ed_m4_gd_scpf(g,scpf) + currentcohort%imort*n_perm2*AREA + ed_m5_gd_scpf(g,scpf) = ed_m5_gd_scpf(g,scpf) + currentcohort%fmort*n_perm2*AREA + + ! basal area [m2/ha] + ed_ba_gd_scpf(g,scpf) = ed_ba_gd_scpf(g,scpf) + & + 0.25*3.14159*((dbh/100.0)**2.0)*n_perm2*AREA + + ! number density [/ha] + ed_np_gd_scpf(g,scpf) = ed_np_gd_scpf(g,scpf) + AREA*n_perm2 + + ! Growth Incrments must have NaN check and woody check + if(currentCohort%ddbhdt == currentCohort%ddbhdt) then + ed_ddbh_gd_scpf(g,scpf) = ed_ddbh_gd_scpf(g,scpf) + & + currentCohort%ddbhdt*n_perm2*AREA + else + ed_ddbh_gd_scpf(g,scpf) = -999.9 + end if + end if + + end if + + currentCohort => currentCohort%taller enddo ! cohort loop !Patch specific variables that are already calculated diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index d96f5178fe..17d132a847 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -285,7 +285,7 @@ subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature currentCohort%gpp_acc = 0.0_r8 currentCohort%resp_acc = 0.0_r8 - call allocate_live_biomass(currentCohort) + call allocate_live_biomass(currentCohort,1) currentCohort => currentCohort%taller diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 515086f86b..ab816a9d7c 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -2,7 +2,7 @@ module EDTypesMod use shr_kind_mod , only : r8 => shr_kind_r8; use decompMod , only : bounds_type - use clm_varpar , only : nlevcan_ed, nclmax, numrad, nlevgrnd + use clm_varpar , only : nlevcan_ed, nclmax, numrad, nlevgrnd, mxpft use domainMod , only : domain_type use shr_sys_mod , only : shr_sys_flush @@ -23,9 +23,11 @@ module EDTypesMod ! for setting number of patches per gridcell and number of cohorts per patch ! for I/O and converting to a vector + integer, parameter :: numPatchesPerGridCell = 10 ! integer, parameter :: numCohortsPerPatch = 160 ! integer, parameter :: cohorts_per_gcell = 1600 ! This is the max number of individual items one can store per + ! each grid cell and effects the striding in the ED restart ! data as some fields are arrays where each array is ! associated with one cohort @@ -41,6 +43,7 @@ module EDTypesMod integer , parameter :: numpft_ed = 2 ! number of PFTs used in ED. integer , parameter :: maxPft = 79 ! max number of PFTs potentially used by CLM + ! SPITFIRE integer , parameter :: NLSC = 6 ! number carbon compartments in above ground litter array integer , parameter :: NFSC = 6 ! number fuel size classes @@ -62,7 +65,39 @@ module EDTypesMod integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI integer , parameter :: N_DBH_BINS = 5 ! no. of dbh bins used when comparing patches - character*4 yearchar + character*4 yearchar + + !the lower limit of the size classes of ED cohorts + !0-10,10-20... + integer, parameter :: nlevsclass_ed = 13 ! Number of dbh size classes for size structure analysis + ! |0-1,1-2,2-3,3-4,4-5,5-10,10-20,20-30,30-40,40-50,50-60,60-70,70-80,80-90,90-100,100+| +! real(r8), parameter, dimension(16) :: sclass_ed = (/0.0_r8,1.0_r8,2.0_r8,3.0_r8,4.0_r8,5.0_r8,10.0_r8,20.0_r8,30.0_r8,40.0_r8, & +! 50.0_r8,60.0_r8,70.0_r8,80.0_r8,90.0_r8,100.0_r8/) + + real(r8), parameter, dimension(13) :: sclass_ed = (/0.0_r8,5.0_r8,10.0_r8,15.0_r8,20.0_r8,30.0_r8,40.0_r8, & + 50.0_r8,60.0_r8,70.0_r8,80.0_r8,90.0_r8,100.0_r8/) + + + ! integer, parameter :: nlevsclass_ed = 17 + ! real(r8), parameter, dimension(17) :: sclass_ed = (/0.1_r8, 5.0_r8,10.0_r8,15.0_r8,20.0_r8,25.0_r8, & + ! 30.0_r8,35.0_r8,40.0_r8,45.0_r8,50.0_r8,55.0_r8, & + ! 60.0_r8,65.0_r8,70.0_r8,75.0_r8,80.0_r8/) + + integer, parameter :: nlevmclass_ed = 5 ! nlev "mortality" classes in ED + ! Number of ways to die + ! (background,hydraulic,carbon,impact,fire) + + character(len = 10), parameter,dimension(5) :: char_list = (/"background","hydraulic ","carbon ","impact ","fire "/) + + + ! These three vectors are used for history output mapping + real(r8) ,allocatable :: levsclass_ed(:) ! The lower bound on size classes for ED trees. This + ! is used really for IO into the + ! history tapes. It gets copied from + ! the parameter array sclass_ed. + integer , allocatable :: pft_levscpf_ed(:) + integer , allocatable :: scls_levscpf_ed(:) + !************************************ !** COHORT type structure ** @@ -102,6 +137,8 @@ module EDTypesMod real(r8) :: c_area ! areal extent of canopy (m2) real(r8) :: treelai ! lai of tree (total leaf area (m2) / canopy area (m2) real(r8) :: treesai ! stem area index of tree (total stem area (m2) / canopy area (m2) + logical :: isnew ! flag to signify a new cohort, new cohorts have not experienced + ! npp or mortality and should therefore not be fused or averaged ! CARBON FLUXES real(r8) :: gpp ! GPP: kgC/indiv/year @@ -114,6 +151,13 @@ module EDTypesMod real(r8) :: resp_acc ! Resp: kgC/indiv/day real(r8) :: resp_clm ! Resp: kgC/indiv/timestep + real(r8) :: npp_leaf ! NPP into leaves (includes replacement of turnover): KgC/indiv/day + real(r8) :: npp_froot ! NPP into fine roots (includes replacement of turnover): KgC/indiv/day + real(r8) :: npp_bsw ! NPP into sapwood: KgC/indiv/day + real(r8) :: npp_bdead ! NPP into deadwood (structure): KgC/indiv/day + real(r8) :: npp_bseed ! NPP into seeds: KgC/indiv/day + real(r8) :: npp_store ! NPP into storage: KgC/indiv/day + real(r8) :: ts_net_uptake(nlevcan_ed) ! Net uptake of leaf layers: kgC/m2/s real(r8) :: year_net_uptake(nlevcan_ed) ! Net uptake of leaf layers: kgC/m2/year @@ -137,6 +181,13 @@ module EDTypesMod !MORTALITY real(r8) :: dmort ! proportional mortality rate. (year-1) + ! Mortality Rate Partitions + real(r8) :: bmort ! background mortality rate n/year + real(r8) :: cmort ! carbon starvation mortality rate n/year + real(r8) :: hmort ! hydraulic failure mortality rate n/year + real(r8) :: imort ! mortality from impacts by others n/year + real(r8) :: fmort ! fire mortality n/year + ! NITROGEN POOLS real(r8) :: livestemn ! live stem nitrogen : KgN/invid real(r8) :: livecrootn ! live coarse root nitrogen: KgN/invid @@ -386,8 +437,42 @@ module EDTypesMod type(userdata), public, target :: udata !-------------------------------------------------------------------------------------! + public :: ed_hist_scpfmaps + contains + !-------------------------------------------------------------------------------------! + subroutine ed_hist_scpfmaps + ! This subroutine allocates and populates the variables + ! that define the mapping of variables in history files in the "scpf" format + ! back to + ! its respective size-class "sc" and pft "pf" + + integer :: i + integer :: isc + integer :: ipft + + allocate( levsclass_ed(1:nlevsclass_ed )) + allocate( pft_levscpf_ed(1:nlevsclass_ed*mxpft)) + allocate(scls_levscpf_ed(1:nlevsclass_ed*mxpft)) + + ! Fill the IO array of plant size classes + ! For some reason the history files did not like + ! a hard allocation of sclass_ed + levsclass_ed(:) = sclass_ed(:) + + ! Fill the IO arrays that match pft and size class to their combined array + i=0 + do ipft=1,mxpft + do isc=1,nlevsclass_ed + i=i+1 + pft_levscpf_ed(i) = ipft + scls_levscpf_ed(i) = isc + end do + end do + + end subroutine ed_hist_scpfmaps + !-------------------------------------------------------------------------------------! function map_clmpatch_to_edpatch(site, clmpatch_number) result(edpatch_pointer) ! diff --git a/components/clm/src/main/histFileMod.F90 b/components/clm/src/main/histFileMod.F90 index 1680f5bdef..1d017fb459 100644 --- a/components/clm/src/main/histFileMod.F90 +++ b/components/clm/src/main/histFileMod.F90 @@ -12,7 +12,7 @@ module histFileMod use shr_sys_mod , only : shr_sys_flush use spmdMod , only : masterproc use abortutils , only : endrun - use clm_varctl , only : iulog, use_vertsoilc + use clm_varctl , only : iulog, use_vertsoilc, use_ed use clm_varcon , only : spval, ispval, dzsoi_decomp use clm_varcon , only : grlnd, nameg, namel, namec, namep, nameCohort use decompMod , only : get_proc_bounds, get_proc_global, bounds_type @@ -21,6 +21,7 @@ module histFileMod use ColumnType , only : col use PatchType , only : patch use ncdio_pio + use EDtypesMod , only : nlevsclass_ed ! implicit none save @@ -1672,7 +1673,7 @@ subroutine htape_create (t, histrest) ! wrapper calls to define the history file contents. ! ! !USES: - use clm_varpar , only : nlevgrnd, nlevsno, nlevlak, nlevurb, numrad, nlevcan + use clm_varpar , only : nlevgrnd, nlevsno, nlevlak, nlevurb, numrad, nlevcan, mxpft use clm_varpar , only : natpft_size, cft_size, maxpatch_glcmec, nlevdecomp_full use landunit_varcon , only : max_lunit use clm_varctl , only : caseid, ctitle, fsurdat, finidat, paramfile @@ -1838,6 +1839,11 @@ subroutine htape_create (t, histrest) call ncd_defdim(lnfid, 'string_length', 8, strlen_dimid) call ncd_defdim( lnfid, 'levdcmp', nlevdecomp_full, dimid) + if(use_ed)then + call ncd_defdim(lnfid, 'levscls', nlevsclass_ed, dimid) + call ncd_defdim(lnfid, 'levscpf', nlevsclass_ed*mxpft, dimid) + end if + if ( .not. lhistrest )then call ncd_defdim(lnfid, 'hist_interval', 2, hist_interval_dimid) call ncd_defdim(lnfid, 'time', ncd_unlimited, time_dimid) @@ -2228,6 +2234,7 @@ subroutine htape_timeconst(t, mode) use domainMod , only : ldomain, lon1d, lat1d use clm_time_manager, only : get_nstep, get_curr_date, get_curr_time use clm_time_manager, only : get_ref_date, get_calendar, NO_LEAP_C, GREGORIAN_C + use EDTypesMod, only : levsclass_ed, pft_levscpf_ed, scls_levscpf_ed ! ! !ARGUMENTS: integer, intent(in) :: t ! tape index @@ -2277,6 +2284,16 @@ subroutine htape_timeconst(t, mode) long_name='coordinate lake levels', units='m', ncid=nfid(t)) call ncd_defvar(varname='levdcmp', xtype=tape(t)%ncprec, dim1name='levdcmp', & long_name='coordinate soil levels', units='m', ncid=nfid(t)) + + if(use_ed)then + call ncd_defvar(varname='levscls', xtype=tape(t)%ncprec, dim1name='levscls', & + long_name='diameter size class lower bound', units='cm', ncid=nfid(t)) + call ncd_defvar(varname='pft_levscpf',xtype=ncd_int, dim1name='levscpf', & + long_name='pft index of the combined pft-size class dimension', units='-', ncid=nfid(t)) + call ncd_defvar(varname='scls_levscpf',xtype=ncd_int, dim1name='levscpf', & + long_name='size index of the combined pft-size class dimension', units='-', ncid=nfid(t)) + end if + elseif (mode == 'write') then if ( masterproc ) write(iulog, *) ' zsoi:',zsoi call ncd_io(varname='levgrnd', data=zsoi, ncid=nfid(t), flag='write') @@ -2287,6 +2304,12 @@ subroutine htape_timeconst(t, mode) zsoi_1d(1) = 1._r8 call ncd_io(varname='levdcmp', data=zsoi_1d, ncid=nfid(t), flag='write') end if + if(use_ed)then + call ncd_io(varname='levscls',data=levsclass_ed, ncid=nfid(t), flag='write') + call ncd_io(varname='pft_levscpf',data=pft_levscpf_ed, ncid=nfid(t), flag='write') + call ncd_io(varname='scls_levscpf',data=scls_levscpf_ed, ncid=nfid(t), flag='write') + end if + endif endif @@ -4266,7 +4289,7 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out, ! ! !USES: use clm_varpar , only : nlevgrnd, nlevsno, nlevlak, numrad, nlevdecomp_full, nlevcan - use clm_varpar , only : natpft_size, cft_size, maxpatch_glcmec + use clm_varpar , only : natpft_size, cft_size, maxpatch_glcmec, mxpft use landunit_varcon , only : max_lunit ! ! !ARGUMENTS: @@ -4344,6 +4367,10 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out, num2d = numrad case ('levdcmp') num2d = nlevdecomp_full + case ('levscls') + num2d = nlevsclass_ed + case ('levscpf') + num2d = nlevsclass_ed*mxpft case('ltype') num2d = max_lunit case('natpft') diff --git a/components/clm/src/main/initVerticalMod.F90 b/components/clm/src/main/initVerticalMod.F90 index f771d6d904..2c91d6be3f 100644 --- a/components/clm/src/main/initVerticalMod.F90 +++ b/components/clm/src/main/initVerticalMod.F90 @@ -16,7 +16,8 @@ module initVerticalMod use clm_varpar , only : toplev_equalspace, nlev_equalspace use clm_varpar , only : nlevsoi, nlevsoifl, nlevurb use clm_varctl , only : fsurdat, iulog - use clm_varctl , only : use_vancouver, use_mexicocity, use_vertsoilc, use_extralakelayers + use clm_varctl , only : use_vancouver, use_mexicocity + use clm_varctl , only : use_vertsoilc, use_extralakelayers, use_ed use clm_varcon , only : zlak, dzlak, zsoi, dzsoi, zisoi, dzsoi_decomp, spval, grlnd use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv use landunit_varcon , only : istdlak, istice_mec @@ -24,6 +25,7 @@ module initVerticalMod use LandunitType , only : lun use ColumnType , only : col use SnowHydrologyMod , only : InitSnowLayers + use EDTypesMod , only : ed_hist_scpfmaps use ncdio_pio ! ! !PUBLIC TYPES: @@ -70,6 +72,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof) real(r8) :: depthratio ! ratio of lake depth to standard deep lake depth integer :: begc, endc integer :: begl, endl + integer :: ipft, isc !------------------------------------------------------------------------ begc = bounds%begc; endc= bounds%endc @@ -174,6 +177,11 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof) end if end if + if(use_ed)then + call ed_hist_scpfmaps + end if + + ! Column level initialization for urban wall and roof layers and interfaces do l = bounds%begl,bounds%endl