From e1fb8190fd288f8e0b611077b408a61ea221b360 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 10 May 2024 11:01:24 -0400 Subject: [PATCH 01/19] removing dependency on small nclmax --- biogeochem/EDCanopyStructureMod.F90 | 174 +++++---------------- biogeochem/EDPatchDynamicsMod.F90 | 1 - biogeochem/EDPhysiologyMod.F90 | 2 +- biogeochem/FatesAllometryMod.F90 | 5 +- biogeochem/FatesCohortMod.F90 | 3 +- biogeochem/FatesPatchMod.F90 | 50 +++--- biogeophys/FatesPlantRespPhotosynthMod.F90 | 13 +- main/EDParamsMod.F90 | 3 +- main/FatesRestartInterfaceMod.F90 | 2 +- radiation/FatesNormanRadMod.F90 | 2 +- radiation/FatesTwoStreamUtilsMod.F90 | 3 - 11 files changed, 87 insertions(+), 171 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 3fbd7b78bd..b38c371c18 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -85,8 +85,6 @@ module EDCanopyStructureMod real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm ! 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 ! ============================================================================ @@ -1546,6 +1544,21 @@ subroutine leaf_area_profile( currentSite ) currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) + ! If the large patch arrays are not new, deallocate them + if(associated(currentPatch%tlai_profile)) then + deallocate(currentPatch%tlai_profile,currentPatch%tsai_profile) + deallocate(currentPatch%elai_profile,currentPatch%esai_profile) + end if + + ! Evaluate required array sizes and possibly re-allocate + + max_leaf_layer = maxval(currentPatch%nrad(:,:)) + allocate(currentPatch%tlai_profile(currentPatch%NCL_p,numpft,max_leaf_layer)) + allocate(currentPatch%tsai_profile(currentPatch%NCL_p,numpft,max_leaf_layer)) + allocate(currentPatch%elai_profile(currentPatch%NCL_p,numpft,max_leaf_layer)) + allocate(currentPatch%esai_profile(currentPatch%NCL_p,numpft,max_leaf_layer)) + allocate(currentPatch%canopy_area_profile(currentPatch%NCL_p,numpft,max_leaf_layer)) + ! -------------------------------------------------------------------------------- ! Calculate tree and canopy areas. ! calculate tree lai and sai. @@ -1587,136 +1600,35 @@ 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 - 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 - 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) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - !---~--- - ! Find current crown depth using the allometric function. - !---~--- - 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 - endif - if(currentSite%snow_depth < layer_bottom_height)then - fraction_exposed = 1._r8 - endif - if(currentSite%snow_depth >= layer_bottom_height .and. & - currentSite%snow_depth <= layer_top_height) then !only partly hidden... - 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)) - if(remainder > dinc_vai(iv) )then - 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 * & - 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 * & - fraction_exposed - - currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & - currentCohort%c_area/currentPatch%total_canopy_area - - - end do - - else !if_preserve_b4b - - do iv = 1,currentCohort%NV - - call VegAreaLayer(currentCohort%treelai, & - currentCohort%treesai, & - currentCohort%height, & - iv,currentCohort%nv,currentCohort%pft, & - 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 - - end do - - end if if_preserve_b4b - + do iv = 1,currentCohort%NV + + call VegAreaLayer(currentCohort%treelai, & + currentCohort%treesai, & + currentCohort%height, & + iv,currentCohort%nv,currentCohort%pft, & + 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 + + end do + currentCohort => currentCohort%taller enddo !cohort diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 140c108d66..f9e684022d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -29,7 +29,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/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index e25a722f7a..b13e8e4ee4 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -3301,7 +3301,7 @@ subroutine UpdateRecruitL2FR(csite) end do ! Find the daily mean for each PFT weighted by number and add it to the running mean - do cl = 1,nclmax + do cl = 1,cpatch%ncl_p do ft = 1,numpft if(rec_n(ft,cl)>nearzero)then rec_l2fr0(ft,cl) = rec_l2fr0(ft,cl) / rec_n(ft,cl) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 1aa1cc4525..7aaf415719 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 @@ -670,7 +669,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 @@ -803,7 +802,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 d86e5c5d51..773b6c409a 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -99,21 +99,27 @@ 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 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) @@ -129,20 +135,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] @@ -321,7 +327,7 @@ subroutine NanValues(this) this%c_stomata = nan this%c_lblayer = nan - this%psn_z(:,:,:) = nan + !this%psn_z(:,:,:) = nan this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = nan this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = nan @@ -412,7 +418,7 @@ subroutine ZeroValues(this) this%elai_profile(:,:,:) = 0.0_r8 this%c_stomata = 0.0_r8 this%c_lblayer = 0.0_r8 - this%psn_z(:,:,:) = 0.0_r8 + ! this%psn_z(:,:,:) = 0.0_r8 this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0.0_r8 this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0.0_r8 diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 0aad6eb977..673e19bce6 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -186,13 +186,16 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! leaf maintenance (dark) respiration [umol CO2/m**2/s] real(r8) :: lmr_z(nlevleaf,maxpft,nclmax) - + ! stomatal resistance [s/m] real(r8) :: rs_z(nlevleaf,maxpft,nclmax) ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s] real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax) + ! Photosynthesis [umol /m2 /s] + real(r8) :: psn_z(nlevleaf,maxpft,nclmax) + ! Mask used to determine which leaf-layer biophysical rates have been ! used already logical :: rate_mask_z(nlevleaf,maxpft,nclmax) @@ -723,7 +726,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) lmr_z(iv,ft,cl), & ! in leaf_psi, & ! in bc_in(s)%rb_pa(ifp), & ! in - currentPatch%psn_z(cl,ft,iv), & ! out + psn_z(iv,ft,cl), & ! out rs_z(iv,ft,cl), & ! out anet_av_z(iv,ft,cl), & ! out c13disc_z(cl,ft,iv)) ! out @@ -755,7 +758,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(radiation_model.eq.norman_solver) then call ScaleLeafLayerFluxToCohort(nv, & !in - currentPatch%psn_z(cl,ft,1:nv), & !in + psn_z(1:nv,ft,cl), & !in lmr_z(1:nv,ft,cl), & !in rs_z(1:nv,ft,cl), & !in currentPatch%elai_profile(cl,ft,1:nv), & !in @@ -773,7 +776,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else call ScaleLeafLayerFluxToCohort(nv, & !in - currentPatch%psn_z(cl,ft,1:nv), & !in + psn_z(1:nv,ft,cl), & !in lmr_z(1:nv,ft,cl), & !in rs_z(1:nv,ft,cl), & !in cohort_layer_elai(1:nv), & !in @@ -1991,7 +1994,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) currentPatch%nrad = currentPatch%ncan ! 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 1389e1c489..4ce1babcdb 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -104,7 +104,8 @@ 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 = 8 ! Maximum number of canopy layers (used only for scratch 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/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index ffef39eb6f..cb808cf379 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -3688,7 +3688,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/FatesTwoStreamUtilsMod.F90 b/radiation/FatesTwoStreamUtilsMod.F90 index 5a87ff24b0..3abe0ce849 100644 --- a/radiation/FatesTwoStreamUtilsMod.F90 +++ b/radiation/FatesTwoStreamUtilsMod.F90 @@ -87,9 +87,6 @@ subroutine FatesConstructRadElements(site,fcansno_pa,coszen_pa) !real(r8), parameter :: init_max_vai_diff_per_elem = 0.2_r8 !type(fates_cohort_type), pointer :: elem_co_ptrs(ncl*max_el_per_layer,100) - - - max_elements = -1 ifp=0 patch => site%oldest_patch From 5fb8e974802e6647cf178192cde8d7b6c5b45465 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 10 May 2024 13:25:50 -0400 Subject: [PATCH 02/19] final pass at the first version of dynamic patch arrays --- biogeochem/EDCanopyStructureMod.F90 | 223 ++++++++++++++++++--------- biogeochem/FatesPatchMod.F90 | 4 +- main/EDParamsMod.F90 | 8 +- radiation/FatesRadiationDriveMod.F90 | 2 +- 4 files changed, 157 insertions(+), 80 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index b38c371c18..1b3474e961 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 @@ -1514,13 +1515,19 @@ 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 + 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 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) @@ -1541,49 +1548,113 @@ 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)) - + cpatch => currentSite%oldest_patch + do while(associated(cpatch)) + + cpatch%ncan(:,:) = 0 + ! This routine updates the %ncan array + call UpdatePatchLAI(cpatch) + + + ! 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. + ! ------------------------------------------------------------------------ + + ncan = cpatch%ncl_p + nveg = maxval(cpatch%ncan(:,:)) + + ! 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(associated(currentPatch%tlai_profile)) then - deallocate(currentPatch%tlai_profile,currentPatch%tsai_profile) - deallocate(currentPatch%elai_profile,currentPatch%esai_profile) + if(associated(cpatch%tlai_profile)) then + + prev_ncan = ubound(cpatch%tlai_profile,1) + prev_nveg = ubound(cpatch%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 nearzero ) then - - call UpdatePatchLAI(currentPatch) - - currentPatch%nrad(:,:) = currentPatch%ncan(:,:) + if_any_canopy_area: if (cpatch%total_canopy_area > nearzero ) then ! ----------------------------------------------------------------------------- ! Standard canopy layering model. @@ -1591,7 +1662,7 @@ subroutine leaf_area_profile( currentSite ) ! and canopy area to the accumulators. ! ----------------------------------------------------------------------------- - currentCohort => currentPatch%shortest + currentCohort => cpatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft cl = currentCohort%canopy_layer @@ -1612,20 +1683,20 @@ subroutine leaf_area_profile( currentSite ) 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 + cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + & + tlai_layer * currentCohort%c_area/cpatch%total_canopy_area - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & - elai_layer * currentCohort%c_area/currentPatch%total_canopy_area + cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + & + elai_layer * currentCohort%c_area/cpatch%total_canopy_area - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & - tsai_layer * currentCohort%c_area/currentPatch%total_canopy_area + cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + & + tsai_layer * currentCohort%c_area/cpatch%total_canopy_area - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & - esai_layer * currentCohort%c_area/currentPatch%total_canopy_area + cpatch%esai_profile(cl,ft,iv) = cpatch%esai_profile(cl,ft,iv) + & + esai_layer * currentCohort%c_area/cpatch%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%canopy_area_profile(cl,ft,iv) = cpatch%canopy_area_profile(cl,ft,iv) + & + currentCohort%c_area/cpatch%total_canopy_area end do @@ -1639,20 +1710,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 @@ -1672,24 +1743,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%ncan(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 @@ -1698,21 +1769,21 @@ subroutine leaf_area_profile( currentSite ) end do do ft = 1,numpft - do iv = 1,currentPatch%ncan(cl,ft) + do iv = 1,cpatch%ncan(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 @@ -1728,23 +1799,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 @@ -1752,7 +1823,7 @@ subroutine leaf_area_profile( currentSite ) end if if_any_canopy_area - currentPatch => currentPatch%younger + cpatch => cpatch%younger enddo !patch @@ -2099,21 +2170,25 @@ subroutine UpdatePatchLAI(currentPatch) type(fates_cohort_type), pointer :: currentCohort integer :: cl ! Canopy layer index integer :: ft ! Plant functional type index - + + integer :: ncl_p + ! 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. - + + currentPatch%ncl_p = 0 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 + currentPatch%ncl_p = max(currentPatch%ncl_p,cl) ! Update the cohort level lai and related variables call UpdateCohortLAI(currentCohort,currentPatch%canopy_layer_tlai, & currentPatch%total_canopy_area) diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 773b6c409a..4eed088dc8 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -120,8 +120,8 @@ module FatesPatchMod 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) :: 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) !--------------------------------------------------------------------------- diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 4ce1babcdb..69d48ed53d 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -104,11 +104,13 @@ 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 = 8 ! Maximum number of canopy layers (used only for scratch arrays) - ! For large arrays at patch level we use dynamic allocation + integer, parameter, public :: nclmax = 5 ! 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 + integer, parameter, public :: nlevleaf = 50 ! number of leaf+stem layers in each canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins diff --git a/radiation/FatesRadiationDriveMod.F90 b/radiation/FatesRadiationDriveMod.F90 index cb642c1289..e8170a74a9 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%ncan(cl,ft)! nlevleaf if(area_vlpfcl(iv,ft,cl) Date: Fri, 10 May 2024 13:40:50 -0400 Subject: [PATCH 03/19] Added ncl_p protections --- biogeochem/EDCanopyStructureMod.F90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 1b3474e961..b1d5b203ec 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -306,9 +306,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(),*) 'FATES generated more canopy layers than scratch-space allows' + write(fates_log(),*) 'Predicted: ',z + write(fates_log(),*) 'Scratch space (nclmax, see EDParamsMod.F90): ',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 From d9f25a003686c2702443dd393150e44e2d59e228 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 10 May 2024 14:03:12 -0400 Subject: [PATCH 04/19] Reverted some changes where we update ncl_p, but it doesn't need to be updated in that spot --- biogeochem/EDCanopyStructureMod.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index b1d5b203ec..fec424299f 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -2181,15 +2181,12 @@ subroutine UpdatePatchLAI(currentPatch) integer :: cl ! Canopy layer index integer :: ft ! Plant functional type index - integer :: ncl_p - ! 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. - currentPatch%ncl_p = 0 canopyloop: do cl = 1,nclmax currentCohort => currentPatch%tallest cohortloop: do while(associated(currentCohort)) @@ -2198,7 +2195,6 @@ subroutine UpdatePatchLAI(currentPatch) if (currentCohort%canopy_layer .eq. cl) then ft = currentCohort%pft - currentPatch%ncl_p = max(currentPatch%ncl_p,cl) ! Update the cohort level lai and related variables call UpdateCohortLAI(currentCohort,currentPatch%canopy_layer_tlai, & currentPatch%total_canopy_area) From 6dc21d1904049de29d4097afb07bd698689c57d7 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 10 May 2024 14:12:16 -0400 Subject: [PATCH 05/19] Revert some removed b4b code, will re-remove later under its own changeset --- biogeochem/EDCanopyStructureMod.F90 | 143 ++++++++++++++++++++++++---- 1 file changed, 124 insertions(+), 19 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index fec424299f..4e1a5f1a05 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -86,6 +86,8 @@ module EDCanopyStructureMod real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm ! 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 ! ============================================================================ @@ -1682,33 +1684,136 @@ subroutine leaf_area_profile( currentSite ) ! This may indeed be zero if there is a sensecent grass ! ---------------------------------------------------------------- - do iv = 1,currentCohort%NV + ! 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 - call VegAreaLayer(currentCohort%treelai, & - currentCohort%treesai, & - currentCohort%height, & - iv,currentCohort%nv,currentCohort%pft, & - currentSite%snow_depth, & - vai_top,vai_bot, & - elai_layer,esai_layer,tlai_layer,tsai_layer) + 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 - cpatch%tlai_profile(cl,ft,iv) = cpatch%tlai_profile(cl,ft,iv) + & - tlai_layer * currentCohort%c_area/cpatch%total_canopy_area + cpatch%nrad(cl,ft) = cpatch%ncan(cl,ft) - cpatch%elai_profile(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) + & - elai_layer * currentCohort%c_area/cpatch%total_canopy_area + 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(), *) ' cpatch%nrad(cl,ft): ', cpatch%nrad(cl,ft) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - cpatch%tsai_profile(cl,ft,iv) = cpatch%tsai_profile(cl,ft,iv) + & - tsai_layer * currentCohort%c_area/cpatch%total_canopy_area + !---~--- + ! Find current crown depth using the allometric function. + !---~--- + call CrownDepth(currentCohort%height,currentCohort%pft,crown_depth) + !---~--- - 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 + ! -------------------------------------------------------------------------- + ! 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. + ! -------------------------------------------------------------------------- - end do + 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 + endif + if(currentSite%snow_depth < layer_bottom_height)then + fraction_exposed = 1._r8 + endif + if(currentSite%snow_depth >= layer_bottom_height .and. & + currentSite%snow_depth <= layer_top_height) then !only partly hidden... + 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)) + if(remainder > dinc_vai(iv) )then + 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 + + 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 + + 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 + + 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, & + iv,currentCohort%nv,currentCohort%pft, & + currentSite%snow_depth, & + vai_top,vai_bot, & + elai_layer,esai_layer,tlai_layer,tsai_layer) + + + 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 From 556edd48452459b7ffbe248adfb921b45205af53 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 16 May 2024 12:55:42 -0600 Subject: [PATCH 06/19] nclmax nlevleaf: reverting for b4b testing --- biogeochem/EDCanopyStructureMod.F90 | 24 +++++++++++++++--------- biogeochem/FatesPatchMod.F90 | 2 +- main/EDParamsMod.F90 | 4 ++-- 3 files changed, 18 insertions(+), 12 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 4e1a5f1a05..a0facdbd7c 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1599,17 +1599,24 @@ subroutine leaf_area_profile( currentSite ) if( prev_ncan .ne. ncan .or. (nveg>prev_nveg) .or. (nveg Date: Wed, 22 May 2024 11:58:16 -0600 Subject: [PATCH 07/19] added memory cleanup to new patch allocated arrays --- biogeochem/FatesPatchMod.F90 | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 8901bfdd41..a8cd18322f 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -316,6 +316,7 @@ subroutine NanValues(this) this%total_canopy_area = nan this%total_tree_area = nan this%zstar = nan + this%elai_profile(:,:,:) = nan this%esai_profile(:,:,:) = nan this%tlai_profile(:,:,:) = nan @@ -679,6 +680,33 @@ subroutine FreeMemory(this, regeneration_model, numpft) this%fragmentation_scaler, & stat=istat, errmsg=smsg) + ! These arrays are allocated in EDCanopyStructureMod + ! while determining how many canopy and leaf layers + ! the patch has. Its possible that patches may + ! be spawned and destroyed before ever reaching + ! that routine, thus we must check to see + ! if the they are already allocated. + if(associated(this%tlai_profile)) then + deallocate(this%tlai_profile) + deallocate(this%tsai_profile) + deallocate(this%elai_profile) + deallocate(this%esai_profile) + deallocate(this%f_sun) + deallocate(this%fabd_sun_z) + deallocate(this%fabd_sha_z) + deallocate(this%fabi_sun_z) + deallocate(this%fabi_sha_z) + deallocate(this%nrmlzd_parprof_pft_dir_z) + deallocate(this%nrmlzd_parprof_pft_dif_z) + deallocate(this%ed_parsun_z) + deallocate(this%ed_parsha_z) + deallocate(this%ed_laisun_z) + deallocate(this%ed_laisha_z) + deallocate(this%parprof_pft_dir_z) + deallocate(this%parprof_pft_dif_z) + deallocate(this%canopy_area_profile) + end if + if (istat/=0) then write(fates_log(),*) 'dealloc009: fail on deallocate patch vectors:'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) From 40e67cfc5386220a1df22496c2180b429a92b4b5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 23 May 2024 11:24:13 -0400 Subject: [PATCH 08/19] Moved patch dynamic reallocations, zeroing and nanning into FatesPatchMod --- biogeochem/EDCanopyStructureMod.F90 | 94 ++----------- biogeochem/FatesPatchMod.F90 | 204 ++++++++++++++++++++++------ 2 files changed, 173 insertions(+), 125 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a0facdbd7c..d944017198 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1534,12 +1534,7 @@ subroutine leaf_area_profile( currentSite ) integer :: ft ! Plant functional type index. integer :: iv ! Vertical leaf layer index integer :: cl ! Canopy layer index - 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 + 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) @@ -1568,86 +1563,15 @@ subroutine leaf_area_profile( currentSite ) call UpdatePatchLAI(cpatch) - ! 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. - ! ------------------------------------------------------------------------ - - ncan = cpatch%ncl_p - nveg = maxval(cpatch%ncan(:,:)) - - ! 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(associated(cpatch%tlai_profile)) then - - prev_ncan = ubound(cpatch%tlai_profile,1) - prev_nveg = ubound(cpatch%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 shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg @@ -105,7 +106,7 @@ module FatesPatchMod ! 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), pointer :: tlai_profile(:,:,:) ! nclmax,maxpft,nlevleaf) + 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) @@ -231,8 +232,11 @@ module FatesPatchMod contains procedure :: Init - procedure :: NanValues + procedure :: NanValues + procedure :: NanDynamics procedure :: ZeroValues + procedure :: ZeroDynamics + procedure :: ReAllocateDynamics procedure :: InitRunningMeans procedure :: InitLitter procedure :: Create @@ -278,6 +282,137 @@ 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%ncan(:,:)) + + ! 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 Date: Thu, 23 May 2024 15:23:55 -0400 Subject: [PATCH 09/19] Updated nclmax and nlevleaf --- biogeochem/EDCanopyStructureMod.F90 | 4 ++-- main/EDParamsMod.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index d944017198..a101065cf6 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -311,9 +311,9 @@ subroutine canopy_structure( currentSite , bc_in ) ! Save number of canopy layers to the patch structure if(z > nclmax) then - write(fates_log(),*) 'FATES generated more canopy layers than scratch-space allows' + write(fates_log(),*) 'Termination should have ensured number of canopy layers was not larger than nclmax' write(fates_log(),*) 'Predicted: ',z - write(fates_log(),*) 'Scratch space (nclmax, see EDParamsMod.F90): ',nclmax + 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__)) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index e89257a3f6..a125c81413 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -104,13 +104,13 @@ 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 !5 ! Maximum number of canopy layers (used only for scratch arrays) + integer, parameter, public :: nclmax = 3 ! 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 !50 ! number of leaf+stem layers in each canopy layer + integer, parameter, public :: nlevleaf = 60 ! number of leaf+stem layers in each canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins From a2b3eada362248de87116f93da31c7337c700ce7 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 9 Jul 2024 14:06:20 -0400 Subject: [PATCH 10/19] changed ncan to nleaf --- biogeochem/EDCanopyStructureMod.F90 | 16 ++++++++-------- biogeochem/FatesPatchMod.F90 | 8 ++++---- biogeophys/FatesPlantRespPhotosynthMod.F90 | 16 ++++++++-------- main/FatesHistoryInterfaceMod.F90 | 2 +- main/FatesInterfaceMod.F90 | 2 +- radiation/FatesRadiationDriveMod.F90 | 2 +- 6 files changed, 23 insertions(+), 23 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a101065cf6..7f2db09cd2 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1500,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 @@ -1558,8 +1558,8 @@ subroutine leaf_area_profile( currentSite ) cpatch => currentSite%oldest_patch do while(associated(cpatch)) - cpatch%ncan(:,:) = 0 - ! This routine updates the %ncan array + cpatch%nleaf(:,:) = 0 + ! This routine updates the %nleaf array call UpdatePatchLAI(cpatch) @@ -1589,7 +1589,7 @@ subroutine leaf_area_profile( currentSite ) ! 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%ncan(:,:) + cpatch%nrad(:,:) = cpatch%nleaf(:,:) ! ------------------------------------------------------------------------------ ! It is remotely possible that in deserts we will not have any canopy @@ -1629,7 +1629,7 @@ subroutine leaf_area_profile( currentSite ) fleaf = 0._r8 endif - cpatch%nrad(cl,ft) = cpatch%ncan(cl,ft) + 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' @@ -1789,7 +1789,7 @@ subroutine leaf_area_profile( currentSite ) ! -------------------------------------------------------------------------- do cl = 1,cpatch%NCL_p - do iv = 1,cpatch%ncan(cl,ft) + do iv = 1,cpatch%nleaf(cl,ft) if( debug .and. sum(cpatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then @@ -1814,7 +1814,7 @@ subroutine leaf_area_profile( currentSite ) end do do ft = 1,numpft - do iv = 1,cpatch%ncan(cl,ft) + do iv = 1,cpatch%nleaf(cl,ft) if( cpatch%canopy_area_profile(cl,ft,iv) > nearzero ) then @@ -2235,7 +2235,7 @@ subroutine UpdatePatchLAI(currentPatch) 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/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index a660c725cd..952f5a7a0d 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -117,8 +117,8 @@ module FatesPatchMod ! 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] @@ -307,7 +307,7 @@ subroutine ReAllocateDynamics(this) ! as defined in the allocation space ncan = this%ncl_p - nveg = maxval(this%ncan(:,:)) + nveg = maxval(this%nleaf(:,:)) ! Assume we will need to allocate, unless the ! arrays already are allocated and require the same size @@ -456,7 +456,7 @@ subroutine NanValues(this) this%canopy_mask(:,:) = fates_unset_int this%nrad(:,:) = fates_unset_int - this%ncan(:,:) = fates_unset_int + this%nleaf(:,:) = fates_unset_int this%c_stomata = nan this%c_lblayer = nan diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 16f242a0f2..b338ab5b51 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -372,7 +372,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Part III. Calculate the number of sublayers for each pft and layer. ! And then identify which layer/pft combinations have things in them. ! Output: - ! currentPatch%ncan(:,:) + ! currentPatch%nleaf(:,:) ! currentPatch%canopy_mask(:,:) call UpdateCanopyNCanNRadPresent(currentPatch) @@ -1978,16 +1978,16 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! --------------------------------------------------------------------------------- ! This subroutine calculates two patch level quanities: - ! currentPatch%ncan and + ! currentPatch%nleaf and ! currentPatch%canopy_mask ! - ! currentPatch%ncan(:,:) is a two dimensional array that indicates + ! currentPatch%nleaf(:,:) is a two dimensional array that indicates ! the total number of leaf layers (including those that are not exposed to light) ! in each canopy layer and for each functional type. ! ! currentPatch%nrad(:,:) is a two dimensional array that indicates ! the total number of EXPOSED leaf layers, but for all intents and purposes - ! in the photosynthesis routine, this appears to be the same as %ncan... + ! in the photosynthesis routine, this appears to be the same as %nleaf... ! ! currentPatch%canopy_mask(:,:) has the same dimensions, is binary, and ! indicates whether or not leaf layers are present (by evaluating the canopy area @@ -2008,14 +2008,14 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! of the layer/pft index it is in ! --------------------------------------------------------------------------------- - currentPatch%ncan(:,:) = 0 + currentPatch%nleaf(:,:) = 0 ! redo the canopy structure algorithm to get round a ! bug that is happening for site 125, FT13. currentCohort => 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 @@ -2023,7 +2023,7 @@ 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,currentPatch%ncl_p 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/radiation/FatesRadiationDriveMod.F90 b/radiation/FatesRadiationDriveMod.F90 index e8170a74a9..eee478446c 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, cpatch%ncan(cl,ft)! nlevleaf + do_iv: do iv = 1, cpatch%nleaf(cl,ft)! nlevleaf if(area_vlpfcl(iv,ft,cl) Date: Wed, 10 Jul 2024 15:07:30 -0600 Subject: [PATCH 11/19] fix on recruit l2fr --- biogeochem/EDPhysiologyMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index daeffc98ff..0f080127c9 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -3320,7 +3320,7 @@ subroutine UpdateRecruitL2FR(csite) end do ! Find the daily mean for each PFT weighted by number and add it to the running mean - do cl = 1,cpatch%ncl_p + do cl = 1,nclmax do ft = 1,numpft if(rec_n(ft,cl)>nearzero)then rec_l2fr0(ft,cl) = rec_l2fr0(ft,cl) / rec_n(ft,cl) From d71c9f02ce2685a80812301c2c2785bffa16586c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 23 Jul 2024 19:42:07 -0600 Subject: [PATCH 12/19] reverting nclmax to 3 for tests --- main/EDParamsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 2a57fe8ace..7defaa7f60 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -110,7 +110,7 @@ module EDParamsMod ! 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 = 60 ! number of leaf+stem layers in each canopy layer + integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in each canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins From cccfac6a8edbbc6644eb592984d890acf04e15b3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 25 Jul 2024 09:04:11 -0600 Subject: [PATCH 13/19] More zeroing of dynamic arrays --- biogeochem/FatesPatchMod.F90 | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 952f5a7a0d..a1260517e7 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -535,6 +535,21 @@ subroutine ZeroDynamics(this) this%fabi_sha_z(:,:,:) = 0._r8 this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8 this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8 + + ! Added + this%elai_profile(:,:,:) = 0._r8 + this%esai_profile(:,:,:) = 0._r8 + this%tlai_profile(:,:,:) = 0._r8 + this%tsai_profile(:,:,:) = 0._r8 + this%canopy_area_profile(:,:,:) = 0._r8 + + this%ed_laisun_z(:,:,:) = 0._r8 + this%ed_laisha_z(:,:,:) = 0._r8 + this%ed_parsun_z(:,:,:) = 0._r8 + this%ed_parsha_z(:,:,:) = 0._r8 + + this%parprof_pft_dir_z(:,:,:) = 0._r8 + this%parprof_pft_dif_z(:,:,:) = 0._r8 end subroutine ZeroDynamics From 248595ca5e6ac0294928138645e7a2c010b95d06 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 26 Jul 2024 11:02:25 -0600 Subject: [PATCH 14/19] resetting nclmax to 2, since there is a bug that is preventing nclmax=3 in main --- main/EDParamsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 7defaa7f60..cc906fecef 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -104,7 +104,7 @@ 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 = 3 ! Maximum number of canopy layers (used only for scratch arrays) + 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 From 0897ba1cbcff1c9375b7126c5f2c0fa3e2ce6413 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Jul 2024 10:59:48 -0600 Subject: [PATCH 15/19] attempts to get b4b, unsuccesfull so far --- biogeochem/EDCanopyStructureMod.F90 | 17 ++++++++++++----- biogeochem/FatesPatchMod.F90 | 12 ++++++------ radiation/FatesNormanRadMod.F90 | 2 +- radiation/FatesRadiationDriveMod.F90 | 2 +- 4 files changed, 20 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 7f2db09cd2..d88f6057fe 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1557,12 +1557,12 @@ subroutine leaf_area_profile( currentSite ) cpatch => currentSite%oldest_patch do while(associated(cpatch)) - + cpatch%nleaf(:,:) = 0 ! 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 @@ -1589,14 +1589,21 @@ subroutine leaf_area_profile( currentSite ) ! 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(:,:) + !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 (cpatch%total_canopy_area > nearzero ) then + if_any_canopy_area: if (cpatch%total_canopy_area <= nearzero ) then + + cpatch%nleaf(:,:) = 0 + cpatch%nrad(:,:) = 0 + + else + + cpatch%nrad(:,:) = cpatch%nleaf(:,:) ! ----------------------------------------------------------------------------- ! Standard canopy layering model. diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index a1260517e7..0e3e8427f9 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -394,10 +394,11 @@ subroutine NanDynamics(this) this%elai_profile(:,:,:) = nan this%esai_profile(:,:,:) = nan this%tlai_profile(:,:,:) = nan - this%tsai_profile(:,:,:) = nan + this%tsai_profile(:,:,:) = nan this%canopy_area_profile(:,:,:) = nan this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = nan this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = nan + this%fabd_sun_z(:,:,:) = nan this%fabd_sha_z(:,:,:) = nan this%fabi_sun_z(:,:,:) = nan @@ -538,16 +539,15 @@ subroutine ZeroDynamics(this) ! Added this%elai_profile(:,:,:) = 0._r8 - this%esai_profile(:,:,:) = 0._r8 - this%tlai_profile(:,:,:) = 0._r8 - this%tsai_profile(:,:,:) = 0._r8 - this%canopy_area_profile(:,:,:) = 0._r8 + !this%esai_profile(:,:,:) = 0._r8 + !this%tlai_profile(:,:,:) = 0._r8 + !this%tsai_profile(:,:,:) = 0._r8 + !this%canopy_area_profile(:,:,:) = 0._r8 this%ed_laisun_z(:,:,:) = 0._r8 this%ed_laisha_z(:,:,:) = 0._r8 this%ed_parsun_z(:,:,:) = 0._r8 this%ed_parsha_z(:,:,:) = 0._r8 - this%parprof_pft_dir_z(:,:,:) = 0._r8 this%parprof_pft_dif_z(:,:,:) = 0._r8 diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index ae73748a42..2310f5ca54 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,currentPatch%ncl_p + do L = 1,nclmax !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 eee478446c..8718cf43b7 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, cpatch%nleaf(cl,ft)! nlevleaf + do_iv: do iv = 1, nlevleaf !cpatch%nleaf(cl,ft) if(area_vlpfcl(iv,ft,cl) Date: Tue, 30 Jul 2024 12:00:18 -0600 Subject: [PATCH 16/19] zeroing psn_z --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index b338ab5b51..b37e3630a0 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -356,6 +356,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) bc_out(s)%rssun_pa(ifp) = 0._r8 bc_out(s)%rssha_pa(ifp) = 0._r8 + psn_z(:,:,:) = 0._r8 + g_sb_leaves = 0._r8 patch_la = 0._r8 From ed62e8203e9b3cd2d7a174d73ea97e75e9ad7427 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Jul 2024 13:57:10 -0600 Subject: [PATCH 17/19] remove canopy closure clause --- biogeochem/EDCanopyStructureMod.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index d88f6057fe..da8c33860a 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1596,12 +1596,12 @@ subroutine leaf_area_profile( currentSite ) ! area, ie not plants at all... ! ------------------------------------------------------------------------------ - if_any_canopy_area: if (cpatch%total_canopy_area <= nearzero ) then +! if_any_canopy_area: if (cpatch%total_canopy_area <= nearzero ) then - cpatch%nleaf(:,:) = 0 - cpatch%nrad(:,:) = 0 +! cpatch%nleaf(:,:) = 0 +! cpatch%nrad(:,:) = 0 - else +! else cpatch%nrad(:,:) = cpatch%nleaf(:,:) @@ -1873,7 +1873,7 @@ subroutine leaf_area_profile( currentSite ) end if - end if if_any_canopy_area + !end if if_any_canopy_area cpatch => cpatch%younger From 9b6b33611a1f88bb7e18040d195b7214e0f34127 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 31 Jul 2024 09:00:53 -0600 Subject: [PATCH 18/19] fixes to summed lai --- biogeochem/EDCanopyStructureMod.F90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index da8c33860a..42402ce566 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1559,6 +1559,7 @@ subroutine leaf_area_profile( currentSite ) do while(associated(cpatch)) cpatch%nleaf(:,:) = 0 + cpatch%canopy_layer_tlai(:) = 0._r8 ! This routine updates the %nleaf array call UpdatePatchLAI(cpatch) @@ -1577,7 +1578,6 @@ subroutine leaf_area_profile( currentSite ) ! calculate tree lai and sai. ! -------------------------------------------------------------------------------- - cpatch%canopy_layer_tlai(:) = 0._r8 cpatch%tlai_profile(:,:,:) = 0._r8 cpatch%tsai_profile(:,:,:) = 0._r8 cpatch%elai_profile(:,:,:) = 0._r8 @@ -1589,22 +1589,15 @@ subroutine leaf_area_profile( currentSite ) ! 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(:,:) + 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 (cpatch%total_canopy_area <= nearzero ) then + if_any_canopy_area: if (cpatch%total_canopy_area > nearzero ) then -! cpatch%nleaf(:,:) = 0 -! cpatch%nrad(:,:) = 0 - -! else - - cpatch%nrad(:,:) = cpatch%nleaf(:,:) - ! ----------------------------------------------------------------------------- ! Standard canopy layering model. ! Go through all cohorts and add their leaf area @@ -1873,7 +1866,7 @@ subroutine leaf_area_profile( currentSite ) end if - !end if if_any_canopy_area + end if if_any_canopy_area cpatch => cpatch%younger From 8da0669aac08dbb84c1a7d8762ff81597e1fa5a8 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 31 Jul 2024 09:40:56 -0600 Subject: [PATCH 19/19] small additions to patch arrays --- biogeochem/EDCanopyStructureMod.F90 | 12 +----------- biogeochem/FatesPatchMod.F90 | 8 ++++---- radiation/FatesNormanRadMod.F90 | 2 +- radiation/FatesRadiationDriveMod.F90 | 2 +- 4 files changed, 7 insertions(+), 17 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 42402ce566..e859a7339f 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1573,17 +1573,7 @@ subroutine leaf_area_profile( currentSite ) call cpatch%NanDynamics() call cpatch%ZeroDynamics() - ! -------------------------------------------------------------------------------- - ! Calculate tree and canopy areas. - ! calculate tree lai and sai. - ! -------------------------------------------------------------------------------- - - cpatch%tlai_profile(:,:,:) = 0._r8 - cpatch%tsai_profile(:,:,:) = 0._r8 - cpatch%elai_profile(:,:,:) = 0._r8 - cpatch%esai_profile(:,:,:) = 0._r8 - cpatch%canopy_area_profile(:,:,:) = 0._r8 - cpatch%canopy_mask(:,:) = 0 + cpatch%canopy_mask(:,:) = 0 ! TO-DO: NRAD HYPOTHETICALLY WOULDNT INCLUDE LAYERS ! UNDER THE SNOW, BUT WE DONT REALLY USE IT TO FILTER diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 0e3e8427f9..cc6c77ded1 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -539,10 +539,10 @@ subroutine ZeroDynamics(this) ! Added this%elai_profile(:,:,:) = 0._r8 - !this%esai_profile(:,:,:) = 0._r8 - !this%tlai_profile(:,:,:) = 0._r8 - !this%tsai_profile(:,:,:) = 0._r8 - !this%canopy_area_profile(:,:,:) = 0._r8 + this%esai_profile(:,:,:) = 0._r8 + this%tlai_profile(:,:,:) = 0._r8 + this%tsai_profile(:,:,:) = 0._r8 + this%canopy_area_profile(:,:,:) = 0._r8 this%ed_laisun_z(:,:,:) = 0._r8 this%ed_laisha_z(:,:,:) = 0._r8 diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index 2310f5ca54..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 !currentPatch%ncl_p + 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 8718cf43b7..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 !cpatch%nleaf(cl,ft) + do_iv: do iv = 1,cpatch%nleaf(cl,ft) if(area_vlpfcl(iv,ft,cl)