diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 3fbd7b78bd..e859a7339f 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -49,6 +49,7 @@ module EDCanopyStructureMod use PRTGenericMod, only : carbon12_element use FatesTwoStreamUtilsMod, only : FatesConstructRadElements use FatesRadiationMemMod , only : twostr_solver + use FatesRadiationMemMod , only : num_rad_stream_types ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -86,7 +87,7 @@ module EDCanopyStructureMod ! can be roughly considered the same right? logical, parameter :: preserve_b4b = .true. - + ! 10/30/09: Created by Rosie Fisher ! 2017/2018: Modifications and updates by Ryan Knox ! ============================================================================ @@ -307,9 +308,19 @@ subroutine canopy_structure( currentSite , bc_in ) enddo ! do while(area_not_balanced) - ! Set current canopy layer occupancy indicator. - currentPatch%NCL_p = min(nclmax,z) - + ! Save number of canopy layers to the patch structure + + if(z > nclmax) then + write(fates_log(),*) 'Termination should have ensured number of canopy layers was not larger than nclmax' + write(fates_log(),*) 'Predicted: ',z + write(fates_log(),*) 'nclmax: ',nclmax + write(fates_log(),*) 'Consider increasing nclmax if this value is to low' + write(fates_log(),*) 'and you think this number of canopy layers is reasonable.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + else + currentPatch%NCL_p = z + end if + ! ------------------------------------------------------------------------------------------- ! if we are using "strict PPA", then calculate a z_star value as ! the height of the smallest tree in the canopy @@ -1489,7 +1500,7 @@ subroutine leaf_area_profile( currentSite ) ! The following patch level diagnostics are updated here: ! ! currentPatch%canopy_layer_tlai(cl) ! total leaf area index of canopy layer - ! currentPatch%ncan(cl,ft) ! number of vegetation layers needed + ! currentPatch%nleaf(cl,ft) ! number of vegetation layers needed ! ! in this patch's pft/canopy-layer ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include ! ! layers occluded by snow @@ -1516,13 +1527,14 @@ subroutine leaf_area_profile( currentSite ) ! ! !LOCAL VARIABLES: - type (fates_patch_type) , pointer :: currentPatch + type (fates_patch_type) , pointer :: cpatch type (fates_cohort_type) , pointer :: currentCohort real(r8) :: remainder !Thickness of layer at bottom of canopy. real(r8) :: fleaf ! fraction of cohort incepting area that is leaves. integer :: ft ! Plant functional type index. integer :: iv ! Vertical leaf layer index integer :: cl ! Canopy layer index + real(r8) :: fraction_exposed ! how much of this layer is not covered by snow? real(r8) :: frac_canopy(N_HEIGHT_BINS) ! amount of canopy in each height class real(r8) :: minh(N_HEIGHT_BINS) ! minimum height in height class (m) @@ -1543,42 +1555,46 @@ subroutine leaf_area_profile( currentSite ) ! leaf area index above it, irrespective of PFT identity... ! Each leaf is defined by how deep in the canopy it is, in terms of LAI units. (FIX(RF,032414), GB) - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) - - ! -------------------------------------------------------------------------------- - ! Calculate tree and canopy areas. - ! calculate tree lai and sai. - ! -------------------------------------------------------------------------------- + cpatch => currentSite%oldest_patch + do while(associated(cpatch)) + + cpatch%nleaf(:,:) = 0 + cpatch%canopy_layer_tlai(:) = 0._r8 + ! This routine updates the %nleaf array + call UpdatePatchLAI(cpatch) + + + ! This call assesses if the large, dynamically allocated + ! patch arrays need to be allocated for the first time, + ! or resized + call cpatch%ReAllocateDynamics() + + ! These calls NaN and zero the above mentioned arrays + call cpatch%NanDynamics() + call cpatch%ZeroDynamics() + + cpatch%canopy_mask(:,:) = 0 - currentPatch%canopy_layer_tlai(:) = 0._r8 - currentPatch%ncan(:,:) = 0 - currentPatch%nrad(:,:) = 0 - currentPatch%tlai_profile(:,:,:) = 0._r8 - currentPatch%tsai_profile(:,:,:) = 0._r8 - currentPatch%elai_profile(:,:,:) = 0._r8 - currentPatch%esai_profile(:,:,:) = 0._r8 - currentPatch%canopy_area_profile(:,:,:) = 0._r8 - currentPatch%canopy_mask(:,:) = 0 + ! TO-DO: NRAD HYPOTHETICALLY WOULDNT INCLUDE LAYERS + ! UNDER THE SNOW, BUT WE DONT REALLY USE IT TO FILTER + ! THEM OUT. CHECK THE CODE AND CONSIDER REMOVING NRAD + ! ALTOGETHER (RGK 05-2024) + cpatch%nrad(:,:) = cpatch%nleaf(:,:) ! ------------------------------------------------------------------------------ ! It is remotely possible that in deserts we will not have any canopy ! area, ie not plants at all... ! ------------------------------------------------------------------------------ - if_any_canopy_area: if (currentPatch%total_canopy_area > nearzero ) then + if_any_canopy_area: if (cpatch%total_canopy_area > nearzero ) then - call UpdatePatchLAI(currentPatch) - - currentPatch%nrad(:,:) = currentPatch%ncan(:,:) - ! ----------------------------------------------------------------------------- ! Standard canopy layering model. ! Go through all cohorts and add their leaf area ! and canopy area to the accumulators. ! ----------------------------------------------------------------------------- - currentCohort => currentPatch%shortest + currentCohort => cpatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft cl = currentCohort%canopy_layer @@ -1587,29 +1603,31 @@ subroutine leaf_area_profile( currentSite ) ! How much of each tree is stem area index? Assuming that there is ! This may indeed be zero if there is a sensecent grass ! ---------------------------------------------------------------- + ! preserve_b4b will be removed soon. This is kept here to prevent ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) if_preserve_b4b: if(preserve_b4b) then - lai = currentCohort%treelai * currentCohort%c_area/currentPatch%total_canopy_area - sai = currentCohort%treesai * currentCohort%c_area/currentPatch%total_canopy_area + + lai = currentCohort%treelai * currentCohort%c_area/cpatch%total_canopy_area + sai = currentCohort%treesai * currentCohort%c_area/cpatch%total_canopy_area if( (currentCohort%treelai+currentCohort%treesai) > nearzero)then - + ! See issue: https://github.com/NGEET/fates/issues/899 ! fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) fleaf = lai / (lai+sai) else fleaf = 0._r8 endif - - currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) - - if (currentPatch%nrad(cl,ft) > nlevleaf ) then + + cpatch%nrad(cl,ft) = cpatch%nleaf(cl,ft) + + if (cpatch%nrad(cl,ft) > nlevleaf ) then write(fates_log(), *) 'Number of radiative leaf layers is larger' write(fates_log(), *) ' than the maximum allowed.' write(fates_log(), *) ' cl: ',cl write(fates_log(), *) ' ft: ',ft write(fates_log(), *) ' nlevleaf: ',nlevleaf - write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft) + write(fates_log(), *) ' cpatch%nrad(cl,ft): ', cpatch%nrad(cl,ft) call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -1618,26 +1636,26 @@ subroutine leaf_area_profile( currentSite ) !---~--- call CrownDepth(currentCohort%height,currentCohort%pft,crown_depth) !---~--- - - + + ! -------------------------------------------------------------------------- ! Whole layers. Make a weighted average of the leaf area in each layer ! before dividing it by the total area. Fill up layer for whole layers. ! -------------------------------------------------------------------------- - + do iv = 1,currentCohort%NV - + ! This loop builds the arrays that define the effective (not snow covered) ! and total (includes snow covered) area indices for leaves and stems ! We calculate the absolute elevation of each layer to help determine if the layer ! is obscured by snow. - + layer_top_height = currentCohort%height - & ( real(iv-1,r8)/currentCohort%NV * crown_depth ) - + layer_bottom_height = currentCohort%height - & ( real(iv,r8)/currentCohort%NV * crown_depth ) - + fraction_exposed = 1.0_r8 if(currentSite%snow_depth > layer_top_height)then fraction_exposed = 0._r8 @@ -1650,7 +1668,7 @@ subroutine leaf_area_profile( currentSite ) fraction_exposed = 1._r8 - max(0._r8,(min(1.0_r8,(currentSite%snow_depth -layer_bottom_height)/ & (layer_top_height-layer_bottom_height )))) endif - + if(iv==currentCohort%NV) then remainder = (currentCohort%treelai + currentCohort%treesai) - & (dlower_vai(iv) - dinc_vai(iv)) @@ -1658,37 +1676,37 @@ subroutine leaf_area_profile( currentSite ) write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & currentCohort%NV,remainder - + call endrun(msg=errMsg(sourcefile, __LINE__)) endif else remainder = dinc_vai(iv) end if - - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & - remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & - remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & + + cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + & + remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + & + remainder * fleaf * currentCohort%c_area/cpatch%total_canopy_area * & fraction_exposed - - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & - remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & - remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area * & + + cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + & + remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + & + remainder * (1._r8 - fleaf) * currentCohort%c_area/cpatch%total_canopy_area * & fraction_exposed - - currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & - currentCohort%c_area/currentPatch%total_canopy_area - - + + cpatch%canopy_area_profile(cl,ft,iv) = cpatch%canopy_area_profile(cl,ft,iv) + & + currentCohort%c_area/cpatch%total_canopy_area + + end do - + else !if_preserve_b4b - + do iv = 1,currentCohort%NV - + call VegAreaLayer(currentCohort%treelai, & currentCohort%treesai, & currentCohort%height, & @@ -1696,27 +1714,27 @@ subroutine leaf_area_profile( currentSite ) currentSite%snow_depth, & vai_top,vai_bot, & elai_layer,esai_layer,tlai_layer,tsai_layer) - - - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & - tlai_layer * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & - elai_layer * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & - tsai_layer * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & - esai_layer * currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & - currentCohort%c_area/currentPatch%total_canopy_area - + + + cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + & + tlai_layer * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + & + elai_layer * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + & + tsai_layer * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + & + esai_layer * currentCohort%c_area/cpatch%total_canopy_area + + cpatch%canopy_area_profile(cl,ft,iv) = cpatch%canopy_area_profile(cl,ft,iv) + & + currentCohort%c_area/cpatch%total_canopy_area + end do - + end if if_preserve_b4b - + currentCohort => currentCohort%taller enddo !cohort @@ -1727,20 +1745,20 @@ subroutine leaf_area_profile( currentSite ) ! should have a value of exactly 1.0 in its top leaf layer ! -------------------------------------------------------------------------- - if ( (currentPatch%NCL_p > 1) .and. & - (sum(currentPatch%canopy_area_profile(1,:,1)) < 0.9999 )) then + if ( (cpatch%NCL_p > 1) .and. & + (sum(cpatch%canopy_area_profile(1,:,1)) < 0.9999 )) then write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top' write(fates_log(), *) 'cl: ',1 write(fates_log(), *) 'iv: ',1 write(fates_log(), *) 'sum(cpatch%canopy_area_profile(1,:,1)): ', & - sum(currentPatch%canopy_area_profile(1,:,1)) - currentCohort => currentPatch%shortest + sum(cpatch%canopy_area_profile(1,:,1)) + currentCohort => cpatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==1)then write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area + cpatch%total_canopy_area,cpatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area + currentCohort%c_area/cpatch%total_canopy_area endif currentCohort => currentCohort%taller enddo !currentCohort @@ -1760,24 +1778,24 @@ subroutine leaf_area_profile( currentSite ) ! It should never be larger than 1 or less than 0. ! -------------------------------------------------------------------------- - do cl = 1,currentPatch%NCL_p - do iv = 1,currentPatch%ncan(cl,ft) + do cl = 1,cpatch%NCL_p + do iv = 1,cpatch%nleaf(cl,ft) - if( debug .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then + if( debug .and. sum(cpatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' write(fates_log(), *) 'cl: ',cl write(fates_log(), *) 'iv: ',iv write(fates_log(), *) 'sum(cpatch%canopy_area_profile(cl,:,iv)): ', & - sum(currentPatch%canopy_area_profile(cl,:,iv)) - currentCohort => currentPatch%shortest + sum(cpatch%canopy_area_profile(cl,:,iv)) + currentCohort => cpatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==cl)then write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, & currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area + cpatch%total_canopy_area,cpatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area + currentCohort%c_area/cpatch%total_canopy_area endif currentCohort => currentCohort%taller enddo !currentCohort @@ -1786,21 +1804,21 @@ subroutine leaf_area_profile( currentSite ) end do do ft = 1,numpft - do iv = 1,currentPatch%ncan(cl,ft) + do iv = 1,cpatch%nleaf(cl,ft) - if( currentPatch%canopy_area_profile(cl,ft,iv) > nearzero ) then + if( cpatch%canopy_area_profile(cl,ft,iv) > nearzero ) then - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) / & + cpatch%canopy_area_profile(cl,ft,iv) - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) / & + cpatch%canopy_area_profile(cl,ft,iv) - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) / & + cpatch%canopy_area_profile(cl,ft,iv) - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) / & + cpatch%canopy_area_profile(cl,ft,iv) end if enddo @@ -1816,23 +1834,23 @@ subroutine leaf_area_profile( currentSite ) ! Leaving this for the time being. ! -------------------------------------------------------------------------- - currentPatch%canopy_mask(:,:) = 0 + cpatch%canopy_mask(:,:) = 0 ! preserve_b4b will be removed soon. This is kept here to prevent ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) if(preserve_b4b) then - do cl = 1,currentPatch%NCL_p + do cl = 1,cpatch%NCL_p do ft = 1,numpft - do iv = 1, currentPatch%nrad(cl,ft) - if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%canopy_mask(cl,ft) = 1 + do iv = 1, cpatch%nrad(cl,ft) + if(cpatch%canopy_area_profile(cl,ft,iv) > 0._r8)then + cpatch%canopy_mask(cl,ft) = 1 endif end do !iv end do end do else - do cl = 1,currentPatch%NCL_p + do cl = 1,cpatch%NCL_p do ft = 1,numpft - if(currentPatch%canopy_area_profile(cl,ft,1) > 0._r8 ) currentPatch%canopy_mask(cl,ft) = 1 + if(cpatch%canopy_area_profile(cl,ft,1) > 0._r8 ) cpatch%canopy_mask(cl,ft) = 1 end do end do end if @@ -1840,7 +1858,7 @@ subroutine leaf_area_profile( currentSite ) end if if_any_canopy_area - currentPatch => currentPatch%younger + cpatch => cpatch%younger enddo !patch @@ -2187,27 +2205,27 @@ subroutine UpdatePatchLAI(currentPatch) type(fates_cohort_type), pointer :: currentCohort integer :: cl ! Canopy layer index integer :: ft ! Plant functional type index - + ! Calculate LAI of layers above. Because it is possible for some understory cohorts ! to be taller than cohorts in the top canopy layer, we must iterate through the ! patch by canopy layer first. Given that canopy_layer_tlai is a patch level variable ! we could iterate through each cohort in any direction as long as we go down through ! the canopy layers. - + canopyloop: do cl = 1,nclmax currentCohort => currentPatch%tallest cohortloop: do while(associated(currentCohort)) ! Only update the current cohort tree lai if lai of the above layers have been calculated if (currentCohort%canopy_layer .eq. cl) then - ft = currentCohort%pft + ft = currentCohort%pft ! Update the cohort level lai and related variables call UpdateCohortLAI(currentCohort,currentPatch%canopy_layer_tlai, & currentPatch%total_canopy_area) ! Update the number of number of vegetation layers - currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) + currentPatch%nleaf(cl,ft) = max(currentPatch%nleaf(cl,ft),currentCohort%NV) ! Update the patch canopy layer tlai (LAI per canopy area) currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + & diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 586c3a1e3b..d75dbaa211 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -28,7 +28,6 @@ module EDPatchDynamicsMod use EDTypesMod , only : site_fluxdiags_type use EDTypesMod , only : min_patch_area use EDTypesMod , only : min_patch_area_forced - use EDParamsMod , only : nclmax use EDParamsMod , only : regeneration_model use FatesInterfaceTypesMod, only : numpft use FatesConstantsMod , only : dtype_ifall diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 2339df3b1a..3cf1ab36a4 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -98,7 +98,6 @@ module FatesAllometryMod use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : FatesWarn,N2S,A2S,I2S use EDParamsMod , only : nlevleaf,dinc_vai,dlower_vai - use EDParamsMod , only : nclmax use DamageMainMod , only : GetCrownReduction implicit none @@ -676,7 +675,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 real(r8), intent(in) :: c_area ! areal extent of canopy (m2) real(r8), intent(in) :: nplant ! number of individuals in cohort per ha integer, intent(in) :: cl ! canopy layer index - real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of + real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of ! each canopy layer real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at canopy ! top, ref 25C @@ -809,7 +808,7 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_ar real(r8), intent(in) :: c_area ! crown area (m2) real(r8), intent(in) :: nplant ! number of plants integer, intent(in) :: cl ! canopy layer index - real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of + real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of ! each canopy layer real(r8), intent(in) :: treelai ! tree LAI for checking purposes only real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90 index e98c4a345a..6dc6f515da 100644 --- a/biogeochem/FatesCohortMod.F90 +++ b/biogeochem/FatesCohortMod.F90 @@ -6,7 +6,6 @@ module FatesCohortMod use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : ican_upper, ican_ustory use EDParamsMod, only : nlevleaf - use EDParamsMod, only : nclmax use FatesGlobals, only : endrun => fates_endrun use FatesGlobals, only : fates_log use PRTGenericMod, only : max_nleafage @@ -560,7 +559,7 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & real(r8), intent(in) :: ctrim ! fraction of the maximum leaf biomass real(r8), intent(in) :: spread ! how spread crowns are in horizontal space real(r8), intent(in) :: carea ! area of cohort, for SP mode [m2] - real(r8), intent(in) :: can_tlai(nclmax) ! patch-level total LAI of each leaf layer + real(r8), intent(in) :: can_tlai(:) ! patch-level total LAI of each canopy layer real(r8), intent(in) :: elongf_leaf ! leaf elongation factor [fraction] real(r8), intent(in) :: elongf_fnrt ! fine-root "elongation factor" [fraction] real(r8), intent(in) :: elongf_stem ! stem "elongation factor" [fraction] diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index b5489f6e82..cc6c77ded1 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -26,6 +26,7 @@ module FatesPatchMod use FatesRadiationMemMod,only : num_swb use FatesRadiationMemMod,only : num_rad_stream_types use FatesInterfaceTypesMod,only : hlm_hio_ignore_val + use FatesInterfaceTypesMod, only : numpft use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg @@ -100,23 +101,29 @@ module FatesPatchMod real(r8) :: total_canopy_area ! area that is covered by vegetation [m2] real(r8) :: total_tree_area ! area that is covered by woody vegetation [m2] real(r8) :: zstar ! height of smallest canopy tree, only meaningful in "strict PPA" mode [m] - real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer [m2 leaf/m2 contributing crown area] - real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer [m2 leaf/m2 contributing crown area] - real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) - real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) - real(r8) :: canopy_area_profile(nclmax,maxpft,nlevleaf) ! fraction of crown area per canopy area in each layer - ! they will sum to 1.0 in the fully closed canopy layers - ! but only in leaf-layers that contain contributions - ! from all cohorts that donate to canopy_area + + ! exposed leaf area in each canopy layer, pft, and leaf layer [m2 leaf/m2 contributing crown area] + real(r8), allocatable :: elai_profile(:,:,:) ! nclmax,maxpft,nlevleaf) + ! exposed stem area in each canopy layer, pft, and leaf layer [m2 leaf/m2 contributing crown area] + real(r8), allocatable :: esai_profile(:,:,:) ! nclmax,maxpft,nlevleaf) + ! total leaf area (includes that which is under snow-pack) + real(r8), allocatable :: tlai_profile(:,:,:) ! nclmax,maxpft,nlevleaf) + ! total stem area (includes that which is under snow-pack) + real(r8), allocatable :: tsai_profile(:,:,:) ! nclmax,maxpft,nlevleaf) + + real(r8), allocatable :: canopy_area_profile(:,:,:) ! nclmax,maxpft,nlevleaf) ! fraction of crown area per canopy area in each layer + ! they will sum to 1.0 in the fully closed canopy layers + ! but only in leaf-layers that contain contributions + ! from all cohorts that donate to canopy_area + integer :: canopy_mask(nclmax,maxpft) ! is there any of this pft in this canopy layer? - integer :: nrad(nclmax,maxpft) ! number of exposed leaf layers for each canopy layer and pft - integer :: ncan(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft + integer :: nrad(nclmax,maxpft) ! number of exposed vegetation layers for each canopy layer and pft + integer :: nleaf(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft real(r8) :: c_stomata ! mean stomatal conductance of all leaves in the patch [umol/m2/s] real(r8) :: c_lblayer ! mean boundary layer conductance of all leaves in the patch [umol/m2/s] - real(r8) :: psn_z(nclmax,maxpft,nlevleaf) - real(r8) :: nrmlzd_parprof_pft_dir_z(num_rad_stream_types,nclmax,maxpft,nlevleaf) - real(r8) :: nrmlzd_parprof_pft_dif_z(num_rad_stream_types,nclmax,maxpft,nlevleaf) + real(r8),allocatable :: nrmlzd_parprof_pft_dir_z(:,:,:,:) !num_rad_stream_types,nclmax,maxpft,nlevleaf) + real(r8),allocatable :: nrmlzd_parprof_pft_dif_z(:,:,:,:) !num_rad_stream_types,nclmax,maxpft,nlevleaf) !--------------------------------------------------------------------------- @@ -130,20 +137,20 @@ module FatesPatchMod ! organized by canopy layer, pft, and leaf layer - real(r8) :: fabd_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of direct light absorbed [0-1] - real(r8) :: fabd_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of direct light absorbed [0-1] - real(r8) :: fabi_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of indirect light absorbed [0-1] - real(r8) :: fabi_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of indirect light absorbed [0-1] - real(r8) :: ed_parsun_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the sun [W/m2] - real(r8) :: ed_parsha_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the shade [W/m2] - real(r8) :: f_sun(nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun [0-1] - real(r8) :: ed_laisun_z(nclmax,maxpft,nlevleaf) - real(r8) :: ed_laisha_z(nclmax,maxpft,nlevleaf) + real(r8),allocatable :: fabd_sun_z(:,:,:) !nclmax,maxpft,nlevleaf) ! sun fraction of direct light absorbed [0-1] + real(r8),allocatable :: fabd_sha_z(:,:,:) !nclmax,maxpft,nlevleaf) ! shade fraction of direct light absorbed [0-1] + real(r8),allocatable :: fabi_sun_z(:,:,:) !nclmax,maxpft,nlevleaf) ! sun fraction of indirect light absorbed [0-1] + real(r8),allocatable :: fabi_sha_z(:,:,:) !nclmax,maxpft,nlevleaf) ! shade fraction of indirect light absorbed [0-1] + real(r8),allocatable :: ed_parsun_z(:,:,:) !nclmax,maxpft,nlevleaf) ! PAR absorbed in the sun [W/m2] + real(r8),allocatable :: ed_parsha_z(:,:,:) !nclmax,maxpft,nlevleaf) ! PAR absorbed in the shade [W/m2] + real(r8),allocatable :: f_sun(:,:,:) !nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun [0-1] + real(r8),allocatable :: ed_laisun_z(:,:,:) !nclmax,maxpft,nlevleaf) + real(r8),allocatable :: ed_laisha_z(:,:,:) !nclmax,maxpft,nlevleaf) ! radiation profiles for comparison against observations - real(r8) :: parprof_pft_dir_z(nclmax,maxpft,nlevleaf) ! direct-beam PAR profile through canopy, by canopy, PFT, leaf level [W/m2] - real(r8) :: parprof_pft_dif_z(nclmax,maxpft,nlevleaf) ! diffuse PAR profile through canopy, by canopy, PFT, leaf level [W/m2] + real(r8),allocatable :: parprof_pft_dir_z(:,:,:) !nclmax,maxpft,nlevleaf) ! direct-beam PAR profile through canopy, by canopy, PFT, leaf level [W/m2] + real(r8),allocatable :: parprof_pft_dif_z(:,:,:) !nclmax,maxpft,nlevleaf) ! diffuse PAR profile through canopy, by canopy, PFT, leaf level [W/m2] real(r8), allocatable :: tr_soil_dir(:) ! fraction of incoming direct radiation transmitted to the soil as direct, by numSWB [0-1] real(r8), allocatable :: tr_soil_dif(:) ! fraction of incoming diffuse radiation that is transmitted to the soil as diffuse [0-1] @@ -226,8 +233,11 @@ module FatesPatchMod contains procedure :: Init - procedure :: NanValues + procedure :: NanValues + procedure :: NanDynamics procedure :: ZeroValues + procedure :: ZeroDynamics + procedure :: ReAllocateDynamics procedure :: InitRunningMeans procedure :: InitLitter procedure :: Create @@ -273,6 +283,138 @@ end subroutine Init !=========================================================================== + subroutine ReAllocateDynamics(this) + + ! ------------------------------------------------------------------------ + ! Perform allocations and re-allocations of potentially large patch arrays + ! + ! Note that this routine is called at the very end of dynamics, after trimming + ! and after canopy structure. This is important because we want to allocate + ! the array spaces for leaf area and radiation scattering based on the most + ! updated values. Note, that this routine is also called before we re-initialize + ! radiation scattering during restarts. + ! ------------------------------------------------------------------------ + + ! arguments + class(fates_patch_type), intent(inout) :: this ! patch object + + ! locals + logical :: re_allocate ! Should we re-allocate the patch arrays? + integer :: prev_nveg ! Previous number of vegetation layers + integer :: nveg ! Number of vegetation layers + integer :: ncan ! Number of canopy layers + integer :: prev_ncan ! Number of canopy layers previously + ! as defined in the allocation space + + ncan = this%ncl_p + nveg = maxval(this%nleaf(:,:)) + + ! Assume we will need to allocate, unless the + ! arrays already are allocated and require the same size + re_allocate = .true. + + ! If the large patch arrays are not new, deallocate them + if(allocated(this%elai_profile)) then + + prev_ncan = ubound(this%tlai_profile,1) + prev_nveg = ubound(this%tlai_profile,3) + + ! We re-allocate if the number of canopy layers has changed, or + ! if the number of vegetation layers is larger than previously. + ! However, we also re-allocate if the number of vegetation layers + ! is not just smaller than previously allocated, but A GOOD BIT smaller + ! than previously allocated. Why? + ! We do this so that we are not always re-allocating. + + if( prev_ncan .ne. ncan .or. (nveg>prev_nveg) .or. (nveg currentPatch%tallest do while(associated(currentCohort)) - currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft), & + currentPatch%nleaf(currentCohort%canopy_layer,currentCohort%pft) = & + max(currentPatch%nleaf(currentCohort%canopy_layer,currentCohort%pft), & currentCohort%NV) currentCohort => currentCohort%shorter @@ -2020,10 +2025,10 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) enddo !cohort ! NRAD = NCAN ... - currentPatch%nrad = currentPatch%ncan + currentPatch%nrad = currentPatch%nleaf ! Now loop through and identify which layer and pft combo has scattering elements - do cl = 1,nclmax + do cl = 1,currentPatch%ncl_p do ft = 1,numpft currentPatch%canopy_mask(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft); diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index c3d8a76273..cc906fecef 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -104,7 +104,10 @@ module EDParamsMod real(r8), parameter, public :: soil_tfrz_thresh = -2.0_r8 ! Soil temperature threshold below which hydraulic failure mortality is off (non-hydro only) in degrees C - integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers + integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers (used only for scratch arrays) + ! We would make this even higher, but making this + ! a little lower keeps the size down on some output arrays + ! For large arrays at patch level we use dynamic allocation ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in each canopy layer diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index c5b48735f6..fd3c09dfd5 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -5395,7 +5395,7 @@ subroutine update_history_hifrq2(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) do_pft1: do ipft=1,numpft do_canlev1: do ican=1,cpatch%ncl_p - do_leaflev1: do ileaf=1,cpatch%ncan(ican,ipft) + do_leaflev1: do ileaf=1,cpatch%nleaf(ican,ipft) ! calculate where we are on multiplexed dimensions clllpf_indx = ileaf + (ican-1) * nlevleaf + (ipft-1) * nlevleaf * nclmax diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index eea65eaa0a..12eb83bcc0 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -2240,7 +2240,7 @@ subroutine SeedlingParPatch(cpatch, & cl_par = 0._r8 cl_area = 0._r8 do ipft = 1,numpft - iv = cpatch%ncan(cl,ipft) + iv = cpatch%nleaf(cl,ipft) ! Avoid calculating when there are no leaf layers for the given pft in the current canopy layer if (iv .ne. 0) then cl_par = cl_par + cpatch%canopy_area_profile(cl,ipft,1)* & diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 3f7e8375fe..8e8ac61ee5 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -3769,7 +3769,7 @@ subroutine update_3dpatch_radiation(this, nsites, sites, bc_out) currentpatch => sites(s)%oldest_patch do while (associated(currentpatch)) ifp = ifp+1 - + currentPatch%f_sun (:,:,:) = 0._r8 currentPatch%fabd_sun_z (:,:,:) = 0._r8 currentPatch%fabd_sha_z (:,:,:) = 0._r8 diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index 24cf38fbb7..ae73748a42 100644 --- a/radiation/FatesNormanRadMod.F90 +++ b/radiation/FatesNormanRadMod.F90 @@ -189,7 +189,7 @@ subroutine PatchNormanRadiation (currentPatch, & tau_layer(:,:,:,:)=0.0_r8 f_abs(:,:,:,:)=0.0_r8 f_abs_leaf(:,:,:,:)=0._r8 - do L = 1,nclmax + do L = 1,currentPatch%ncl_p do ft = 1,numpft currentPatch%canopy_mask(L,ft) = 0 do iv = 1, currentPatch%nrad(L,ft) diff --git a/radiation/FatesRadiationDriveMod.F90 b/radiation/FatesRadiationDriveMod.F90 index cb642c1289..b3e36c39b0 100644 --- a/radiation/FatesRadiationDriveMod.F90 +++ b/radiation/FatesRadiationDriveMod.F90 @@ -471,7 +471,7 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) end do do_icol do ft = 1,numpft - do_iv: do iv = 1, nlevleaf + do_iv: do iv = 1,cpatch%nleaf(cl,ft) if(area_vlpfcl(iv,ft,cl) site%oldest_patch