From fe08e916b238d3d28cbbae694bc0c999106b4fc5 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Thu, 18 Jan 2018 12:23:11 -0800 Subject: [PATCH 01/24] fixed bug in which patch and cohort numbers were overestimated by one --- main/FatesHistoryInterfaceMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 79e63357ca..33899517f5 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2578,12 +2578,12 @@ subroutine define_history_vars(this, initialize_variables) ! Site level counting variables call this%set_history_var(vname='ED_NPATCHES', units='none', & long='Total number of ED patches per site', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=1.0_r8, upfreq=1, & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_npatches_si) call this%set_history_var(vname='ED_NCOHORTS', units='none', & long='Total number of ED cohorts per site', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=1.0_r8, upfreq=1, & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_ncohorts_si) ! Patch variables From 05620557c5e384c19fa7a28d81a3882b30c5c619 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 22 Jan 2018 20:15:38 -0700 Subject: [PATCH 02/24] added logic and parameter to force fusion of old patches --- biogeochem/EDPatchDynamicsMod.F90 | 71 +++++++++++++++++++------------ 1 file changed, 43 insertions(+), 28 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 051baf689a..738d4ea95a 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1322,6 +1322,10 @@ subroutine fuse_patches( csite, bc_in ) integer :: nopatches !number of patches presently in gridcell integer :: iterate !switch of patch reduction iteration scheme. 1 to keep going, 0 to stop integer :: fuse_flag !do patches get fused (1) or not (0). + ! + ! parameters + real(r8), parameter :: patch_fusion_tolerance_relaxation_increment = 1.1_r8 + real(r8), parameter :: max_age_of_second_oldest_patch = 200._r8 ! age in years above which to combine all patches !--------------------------------------------------------------------- !maxpatch = 4 @@ -1372,33 +1376,44 @@ subroutine fuse_patches( csite, bc_in ) fuse_flag = 1 !the default is to fuse the patches if(currentPatch%patchno /= tpp%patchno) then !these should be the same patch - !---------------------------------------------------------------------! - ! Calculate the difference criteria for each pft and dbh class ! - !---------------------------------------------------------------------! - do ft = 1,numpft ! loop over pfts - do z = 1,n_dbh_bins ! loop over hgt bins - !is there biomass in this category? - if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then - norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& - &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) - !---------------------------------------------------------------------! - ! Look for differences in profile biomass, above the minimum biomass ! - !---------------------------------------------------------------------! - - if(norm > profiletol)then - !looking for differences between profile density. - if(currentPatch%pft_agb_profile(ft,z) > NTOL.or.tpp%pft_agb_profile(ft,z) > NTOL)then - fuse_flag = 0 !do not fuse - keep apart. - endif - endif ! profile tol - endif ! NTOL - enddo !ht bins - enddo ! PFT - - !---------------------------------------------------------------------! - ! Call the patch fusion routine if there is a meaningful difference ! - ! any of the pft x height categories ! - !---------------------------------------------------------------------! + !----------------------------------------------------------------------------------- + ! check to see if both patches are older than the age at which we force them to fuse + !----------------------------------------------------------------------------------- + + if ( tpp%age .le. max_age_of_second_oldest_patch .or. & + currentPatch%age .le. max_age_of_second_oldest_patch ) then + + !---------------------------------------------------------------------! + ! Calculate the difference criteria for each pft and dbh class ! + !---------------------------------------------------------------------! + do ft = 1,numpft ! loop over pfts + do z = 1,n_dbh_bins ! loop over hgt bins + !is there biomass in this category? + if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then + + norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& + &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) + + !---------------------------------------------------------------------! + ! Look for differences in profile biomass, above the minimum biomass ! + !---------------------------------------------------------------------! + + if(norm > profiletol)then + !looking for differences between profile density. + if(currentPatch%pft_agb_profile(ft,z) > NTOL.or.tpp%pft_agb_profile(ft,z) > NTOL)then + fuse_flag = 0 !do not fuse - keep apart. + endif + endif ! profile tol + endif ! NTOL + enddo !ht bins + enddo ! PFT + + endif ! maxage + !-------------------------------------------------------------------------! + ! Call the patch fusion routine if there is not a meaningful difference ! + ! any of the pft x height categories ! + ! or both are older than forced fusion age ! + !-------------------------------------------------------------------------! if(fuse_flag == 1)then tmpptr => currentPatch%older @@ -1434,7 +1449,7 @@ subroutine fuse_patches( csite, bc_in ) if(nopatches > maxpatch)then iterate = 1 - profiletol = profiletol * 1.1_r8 + profiletol = profiletol * patch_fusion_tolerance_relaxation_increment !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! From 798fd29bb16c8b100c76cdb3b0e36a7805a026fa Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 23 Jan 2018 10:31:43 -0700 Subject: [PATCH 03/24] changed bin spacing for patch fusion comparisons to resolve smaller trees better --- biogeochem/EDPatchDynamicsMod.F90 | 15 +++++---------- main/EDTypesMod.F90 | 8 +++++--- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 738d4ea95a..d3a1253649 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -7,7 +7,7 @@ module EDPatchDynamicsMod use FatesInterfaceMod , only : hlm_freq_day use EDPftvarcon , only : EDPftvarcon_inst use EDCohortDynamicsMod , only : fuse_cohorts, sort_cohorts, insert_cohort - use EDtypesMod , only : ncwd, n_dbh_bins, ntol, area, dbhmax + use EDtypesMod , only : ncwd, n_dbh_bins, ntol, area, patchfusion_dbhbin_loweredges use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : min_patch_area @@ -1767,20 +1767,15 @@ subroutine patch_pft_size_profile(cp_pnt) currentPatch => cp_pnt - delta_dbh = (DBHMAX/N_DBH_BINS) - currentPatch%pft_agb_profile(:,:) = 0.0_r8 do j = 1,N_DBH_BINS - if (j == 1) then - mind(j) = 0.0_r8 - maxd(j) = delta_dbh - else if (j == N_DBH_BINS) then - mind(j) = (j-1) * delta_dbh + if (j == N_DBH_BINS) then + mind(j) = patchfusion_dbhbin_loweredges(j-1) maxd(j) = gigantictrees else - mind(j) = (j-1) * delta_dbh - maxd(j) = (j)*delta_dbh + mind(j) = patchfusion_dbhbin_loweredges(j-1) + maxd(j) = patchfusion_dbhbin_loweredges(j) endif enddo diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 05e0775205..383097c072 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -82,12 +82,14 @@ module EDTypesMod real(r8), parameter :: fire_threshold = 50.0_r8 ! threshold for fires that spread or go out. KWm-2 (Pyne 1986) ! PATCH FUSION + integer , parameter :: N_DBH_BINS = 6 ! no. of dbh bins used when comparing patches + real(r8), parameter :: patchfusion_dbhbin_loweredges(N_DBH_BINS) = & + (/0._r8, 5._r8, 20._r8, 50._r8, 100._r8, 150._r8/) ! array of bin lower edges for comparing patches + + ! COHORT FUSION real(r8), parameter :: NTOL = 0.05_r8 ! min plant density for hgt bin to be used in height profile comparisons real(r8), parameter :: HITEMAX = 30.0_r8 ! max dbh value used in hgt profile comparison - real(r8), parameter :: DBHMAX = 150.0_r8 ! max dbh value used in hgt profile comparison integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI - integer , parameter :: N_DBH_BINS = 5 ! no. of dbh bins used when comparing patches - real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination From b7ad51b83061b1b82119c7f0545091c684353067 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 23 Jan 2018 15:03:24 -0700 Subject: [PATCH 04/24] reduced min biomass threshold for patch fusino logic to 5 gC/m2 --- main/EDTypesMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 383097c072..e9220251aa 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -82,12 +82,12 @@ module EDTypesMod real(r8), parameter :: fire_threshold = 50.0_r8 ! threshold for fires that spread or go out. KWm-2 (Pyne 1986) ! PATCH FUSION + real(r8), parameter :: NTOL = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches integer , parameter :: N_DBH_BINS = 6 ! no. of dbh bins used when comparing patches real(r8), parameter :: patchfusion_dbhbin_loweredges(N_DBH_BINS) = & (/0._r8, 5._r8, 20._r8, 50._r8, 100._r8, 150._r8/) ! array of bin lower edges for comparing patches ! COHORT FUSION - real(r8), parameter :: NTOL = 0.05_r8 ! min plant density for hgt bin to be used in height profile comparisons real(r8), parameter :: HITEMAX = 30.0_r8 ! max dbh value used in hgt profile comparison integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI From f5a7dc7caaa8401d90511ce89db29d174d010cd0 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 23 Jan 2018 21:43:28 -0700 Subject: [PATCH 05/24] some cleanup and addition of comments for patch fusion criteria --- biogeochem/EDPatchDynamicsMod.F90 | 52 ++++++++++++++++++++++--------- main/EDTypesMod.F90 | 3 ++ 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index d3a1253649..bcf2eb8c29 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1307,6 +1307,8 @@ subroutine fuse_patches( csite, bc_in ) ! ! !USES: use EDParamsMod , only : ED_val_patch_fusion_tol + use EDTypesMod , only : patch_fusion_tolerance_relaxation_increment + use EDTypesMod , only : max_age_of_second_oldest_patch ! ! !ARGUMENTS: type(ed_site_type), intent(inout), target :: csite @@ -1318,19 +1320,12 @@ subroutine fuse_patches( csite, bc_in ) integer :: ft,z !counters for pft and height class real(r8) :: norm !normalized difference between biomass profiles real(r8) :: profiletol !tolerance of patch fusion routine. Starts off high and is reduced if there are too many patches. - integer :: maxpatch !maximum number of allowed patches. FIX-RF. These should be namelist variables. integer :: nopatches !number of patches presently in gridcell integer :: iterate !switch of patch reduction iteration scheme. 1 to keep going, 0 to stop integer :: fuse_flag !do patches get fused (1) or not (0). ! - ! parameters - real(r8), parameter :: patch_fusion_tolerance_relaxation_increment = 1.1_r8 - real(r8), parameter :: max_age_of_second_oldest_patch = 200._r8 ! age in years above which to combine all patches !--------------------------------------------------------------------- - !maxpatch = 4 - maxpatch = maxPatchesPerSite - currentSite => csite profiletol = ED_val_patch_fusion_tol @@ -1347,7 +1342,7 @@ subroutine fuse_patches( csite, bc_in ) iterate = 1 !---------------------------------------------------------------------! - ! Keep doing this until nopatches >= maxpatch ! + ! Keep doing this until nopatches >= maxPatchesPerSite ! !---------------------------------------------------------------------! do while(iterate == 1) @@ -1373,7 +1368,17 @@ subroutine fuse_patches( csite, bc_in ) endif if(associated(tpp).and.associated(currentPatch))then - fuse_flag = 1 !the default is to fuse the patches + + !-------------------------------------------------------------------------------------------- + ! The default is to fuse the patches, unless some criteria is met which keeps them separated. + ! there are multiple criteria which all need to be met to keep them distinct: + ! (a) one of them is younger than the max age at which we force fusion. + ! (b) for at least one pft x size class, where there is biomass in at least one patch, + ! and the normalized differenve between the patches exceeds a threshold, and + ! (c) there is more than a threshold (tiny) amount of biomass in at least one of the patches + !-------------------------------------------------------------------------------------------- + + fuse_flag = 1 if(currentPatch%patchno /= tpp%patchno) then !these should be the same patch !----------------------------------------------------------------------------------- @@ -1386,11 +1391,20 @@ subroutine fuse_patches( csite, bc_in ) !---------------------------------------------------------------------! ! Calculate the difference criteria for each pft and dbh class ! !---------------------------------------------------------------------! + do ft = 1,numpft ! loop over pfts do z = 1,n_dbh_bins ! loop over hgt bins + + !---------------------------------- !is there biomass in this category? + !---------------------------------- + if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then + !------------------------------------------------------------------------------------- + ! what is the relative difference in biomass i nthis category between the two patches? + !------------------------------------------------------------------------------------- + norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) @@ -1399,16 +1413,24 @@ subroutine fuse_patches( csite, bc_in ) !---------------------------------------------------------------------! if(norm > profiletol)then - !looking for differences between profile density. + + !--------------------------------------------------------------------------------------------------------- + ! the next bit of logic forces fusion of two patches which both have tiny biomass densities. without this, + ! fates gives a bunch of really young patches which all have almost no biomass and so don't need to be + ! distinguished from each other. but if NTOL is too big, it takes too long for the youngest patch to build + ! up enough biomass to be its own distinct entity, which leads to large oscillations in the patch dynamics + ! and dependent variables. + !--------------------------------------------------------------------------------------------------------- + if(currentPatch%pft_agb_profile(ft,z) > NTOL.or.tpp%pft_agb_profile(ft,z) > NTOL)then fuse_flag = 0 !do not fuse - keep apart. - endif + endif ! biomass .gt. NTOL endif ! profile tol - endif ! NTOL + endif ! biomass .gt. 0 enddo !ht bins enddo ! PFT - endif ! maxage + !-------------------------------------------------------------------------! ! Call the patch fusion routine if there is not a meaningful difference ! ! any of the pft x height categories ! @@ -1447,7 +1469,7 @@ subroutine fuse_patches( csite, bc_in ) currentPatch => currentPatch%older enddo - if(nopatches > maxpatch)then + if(nopatches > maxPatchesPerSite)then iterate = 1 profiletol = profiletol * patch_fusion_tolerance_relaxation_increment @@ -1458,7 +1480,7 @@ subroutine fuse_patches( csite, bc_in ) iterate = 0 endif - enddo !do while nopatches>maxpatch + enddo !do while nopatches>maxPatchesPerSite end subroutine fuse_patches diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index e9220251aa..ecc86f0c9c 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -86,11 +86,14 @@ module EDTypesMod integer , parameter :: N_DBH_BINS = 6 ! no. of dbh bins used when comparing patches real(r8), parameter :: patchfusion_dbhbin_loweredges(N_DBH_BINS) = & (/0._r8, 5._r8, 20._r8, 50._r8, 100._r8, 150._r8/) ! array of bin lower edges for comparing patches + real(r8), parameter :: patch_fusion_tolerance_relaxation_increment = 1.1_r8 ! amount by which to increment patch fusion threshold + real(r8), parameter :: max_age_of_second_oldest_patch = 200._r8 ! age in years above which to combine all patches ! COHORT FUSION real(r8), parameter :: HITEMAX = 30.0_r8 ! max dbh value used in hgt profile comparison integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI + ! COHORT TERMINATION real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination real(r8), parameter :: min_nppatch = 1.0E-11_r8 ! minimum number of cohorts per patch (min_npm2*min_patch_area) From 77bc142e4bcea7f9fec0e6b756197940f04fb844 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 24 Jan 2018 09:37:07 -0700 Subject: [PATCH 06/24] reordered patch fusion logic to sum min biomass over PFT,Z and took out of those loops --- biogeochem/EDPatchDynamicsMod.F90 | 74 ++++++++++++++++--------------- 1 file changed, 39 insertions(+), 35 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index bcf2eb8c29..db72291927 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1372,10 +1372,10 @@ subroutine fuse_patches( csite, bc_in ) !-------------------------------------------------------------------------------------------- ! The default is to fuse the patches, unless some criteria is met which keeps them separated. ! there are multiple criteria which all need to be met to keep them distinct: - ! (a) one of them is younger than the max age at which we force fusion. - ! (b) for at least one pft x size class, where there is biomass in at least one patch, - ! and the normalized differenve between the patches exceeds a threshold, and - ! (c) there is more than a threshold (tiny) amount of biomass in at least one of the patches + ! (a) one of them is younger than the max age at which we force fusion; + ! (b) there is more than a threshold (tiny) amount of biomass in at least one of the patches; + ! (c) for at least one pft x size class, where there is biomass in that class in at least one patch, + ! and the normalized difference between the patches exceeds a threshold. !-------------------------------------------------------------------------------------------- fuse_flag = 1 @@ -1388,47 +1388,51 @@ subroutine fuse_patches( csite, bc_in ) if ( tpp%age .le. max_age_of_second_oldest_patch .or. & currentPatch%age .le. max_age_of_second_oldest_patch ) then - !---------------------------------------------------------------------! - ! Calculate the difference criteria for each pft and dbh class ! - !---------------------------------------------------------------------! + + !--------------------------------------------------------------------------------------------------------- + ! the next bit of logic forces fusion of two patches which both have tiny biomass densities. without this, + ! fates gives a bunch of really young patches which all have almost no biomass and so don't need to be + ! distinguished from each other. but if NTOL is too big, it takes too long for the youngest patch to build + ! up enough biomass to be its own distinct entity, which leads to large oscillations in the patch dynamics + ! and dependent variables. + !--------------------------------------------------------------------------------------------------------- + + if(sum(currentPatch%pft_agb_profile(:,:)) > NTOL .or. & + sum(tpp%pft_agb_profile(:,:)) > NTOL ) then - do ft = 1,numpft ! loop over pfts - do z = 1,n_dbh_bins ! loop over hgt bins + !---------------------------------------------------------------------! + ! Calculate the difference criteria for each pft and dbh class ! + !---------------------------------------------------------------------! - !---------------------------------- - !is there biomass in this category? - !---------------------------------- + do ft = 1,numpft ! loop over pfts + do z = 1,n_dbh_bins ! loop over hgt bins - if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then + !---------------------------------- + !is there biomass in this category? + !---------------------------------- - !------------------------------------------------------------------------------------- - ! what is the relative difference in biomass i nthis category between the two patches? - !------------------------------------------------------------------------------------- + if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then - norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& - &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) + !------------------------------------------------------------------------------------- + ! what is the relative difference in biomass i nthis category between the two patches? + !------------------------------------------------------------------------------------- - !---------------------------------------------------------------------! - ! Look for differences in profile biomass, above the minimum biomass ! - !---------------------------------------------------------------------! + norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& + &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) - if(norm > profiletol)then + !---------------------------------------------------------------------! + ! Look for differences in profile biomass, above the minimum biomass ! + !---------------------------------------------------------------------! - !--------------------------------------------------------------------------------------------------------- - ! the next bit of logic forces fusion of two patches which both have tiny biomass densities. without this, - ! fates gives a bunch of really young patches which all have almost no biomass and so don't need to be - ! distinguished from each other. but if NTOL is too big, it takes too long for the youngest patch to build - ! up enough biomass to be its own distinct entity, which leads to large oscillations in the patch dynamics - ! and dependent variables. - !--------------------------------------------------------------------------------------------------------- + if(norm > profiletol)then - if(currentPatch%pft_agb_profile(ft,z) > NTOL.or.tpp%pft_agb_profile(ft,z) > NTOL)then fuse_flag = 0 !do not fuse - keep apart. - endif ! biomass .gt. NTOL - endif ! profile tol - endif ! biomass .gt. 0 - enddo !ht bins - enddo ! PFT + + endif ! profile tol + endif ! biomass(ft,z) .gt. 0 + enddo !ht bins + enddo ! PFT + endif ! sum(biomass(:,:) .gt. NTOL endif ! maxage !-------------------------------------------------------------------------! From c8a76db47dbb5697da7e6f5829aaecd33b0177c5 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 24 Jan 2018 20:42:26 -0700 Subject: [PATCH 07/24] bugfix to dbhbinedges_indexing --- biogeochem/EDPatchDynamicsMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index db72291927..ef2ea36ef9 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1797,11 +1797,11 @@ subroutine patch_pft_size_profile(cp_pnt) do j = 1,N_DBH_BINS if (j == N_DBH_BINS) then - mind(j) = patchfusion_dbhbin_loweredges(j-1) + mind(j) = patchfusion_dbhbin_loweredges(j) maxd(j) = gigantictrees else - mind(j) = patchfusion_dbhbin_loweredges(j-1) - maxd(j) = patchfusion_dbhbin_loweredges(j) + mind(j) = patchfusion_dbhbin_loweredges(j) + maxd(j) = patchfusion_dbhbin_loweredges(j+1) endif enddo From 89eecbdfd1a00f4f6797114c117663da237d0ed2 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 5 Mar 2018 18:17:02 -0700 Subject: [PATCH 08/24] fixed some line length issues --- biogeochem/EDPatchDynamicsMod.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index e89df3705d..0da91d8877 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -118,7 +118,8 @@ subroutine disturbance_rates( site_in, bc_in) call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort) currentCohort%dmort = cmort+hmort+bmort+frmort - call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft,currentCohort%c_area) + call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft, & + currentCohort%c_area) ! Initialize diagnostic mortality rates currentCohort%cmort = cmort @@ -1425,13 +1426,15 @@ subroutine fuse_patches( csite, bc_in ) !is there biomass in this category? !---------------------------------- - if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then + if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8 .or. & + tpp%pft_agb_profile(ft,z) > 0.0_r8)then !------------------------------------------------------------------------------------- ! what is the relative difference in biomass i nthis category between the two patches? !------------------------------------------------------------------------------------- - norm = abs(currentPatch%pft_agb_profile(ft,z) - tpp%pft_agb_profile(ft,z))/(0.5_r8*& + norm = abs(currentPatch%pft_agb_profile(ft,z) - & + tpp%pft_agb_profile(ft,z))/(0.5_r8 * & &(currentPatch%pft_agb_profile(ft,z) + tpp%pft_agb_profile(ft,z))) !---------------------------------------------------------------------! From 8ddcec3f267ae74cbc7e2676293aa5aaeb912548 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 7 Mar 2018 15:29:37 -0700 Subject: [PATCH 09/24] changed NTOL to force_patchfuse_min_biomass --- biogeochem/EDPatchDynamicsMod.F90 | 15 ++++++++------- main/EDTypesMod.F90 | 2 +- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 0da91d8877..15a2b0aed4 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -7,7 +7,8 @@ module EDPatchDynamicsMod use FatesInterfaceMod , only : hlm_freq_day use EDPftvarcon , only : EDPftvarcon_inst use EDCohortDynamicsMod , only : fuse_cohorts, sort_cohorts, insert_cohort - use EDtypesMod , only : ncwd, n_dbh_bins, ntol, area, patchfusion_dbhbin_loweredges + use EDtypesMod , only : ncwd, n_dbh_bins, area, patchfusion_dbhbin_loweredges + use EDtypesMod , only : force_patchfuse_min_biomass use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : min_patch_area @@ -1407,13 +1408,13 @@ subroutine fuse_patches( csite, bc_in ) !--------------------------------------------------------------------------------------------------------- ! the next bit of logic forces fusion of two patches which both have tiny biomass densities. without this, ! fates gives a bunch of really young patches which all have almost no biomass and so don't need to be - ! distinguished from each other. but if NTOL is too big, it takes too long for the youngest patch to build - ! up enough biomass to be its own distinct entity, which leads to large oscillations in the patch dynamics - ! and dependent variables. + ! distinguished from each other. but if force_patchfuse_min_biomass is too big, it takes too long for the + ! youngest patch to build up enough biomass to be its own distinct entity, which leads to large oscillations + ! in the patch dynamics and dependent variables. !--------------------------------------------------------------------------------------------------------- - if(sum(currentPatch%pft_agb_profile(:,:)) > NTOL .or. & - sum(tpp%pft_agb_profile(:,:)) > NTOL ) then + if(sum(currentPatch%pft_agb_profile(:,:)) > force_patchfuse_min_biomass .or. & + sum(tpp%pft_agb_profile(:,:)) > force_patchfuse_min_biomass ) then !---------------------------------------------------------------------! ! Calculate the difference criteria for each pft and dbh class ! @@ -1449,7 +1450,7 @@ subroutine fuse_patches( csite, bc_in ) endif ! biomass(ft,z) .gt. 0 enddo !ht bins enddo ! PFT - endif ! sum(biomass(:,:) .gt. NTOL + endif ! sum(biomass(:,:) .gt. force_patchfuse_min_biomass endif ! maxage !-------------------------------------------------------------------------! diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 45e3ea1ea3..fafce06313 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -81,7 +81,7 @@ module EDTypesMod real(r8), parameter :: fire_threshold = 50.0_r8 ! threshold for fires that spread or go out. KWm-2 (Pyne 1986) ! PATCH FUSION - real(r8), parameter :: NTOL = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches + real(r8), parameter :: force_patchfuse_min_biomass = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches integer , parameter :: N_DBH_BINS = 6 ! no. of dbh bins used when comparing patches real(r8), parameter :: patchfusion_dbhbin_loweredges(N_DBH_BINS) = & (/0._r8, 5._r8, 20._r8, 50._r8, 100._r8, 150._r8/) ! array of bin lower edges for comparing patches From c8f2e8a86cd0a3a709604bea33021b449b10818c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Mar 2018 16:46:38 -0800 Subject: [PATCH 10/24] Code cleaning in Canopy structure. Removed redundant patch variable %canopy_area, which is also %total_canopy_area. Added endruns to various checks which should trigger a fail. --- biogeochem/EDCanopyStructureMod.F90 | 322 +++++++++++---------- biogeochem/EDPatchDynamicsMod.F90 | 5 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 2 +- main/EDTypesMod.F90 | 6 +- 4 files changed, 175 insertions(+), 160 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a2f241caa5..239f8dd31e 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1,9 +1,9 @@ module EDCanopyStructureMod - ! ============================================================================ - ! Code to determine whether the canopy is closed, and which plants are either in the understorey or overstorey - ! This is obviosuly far too complicated for it's own good and needs re-writing. - ! ============================================================================ + ! ===================================================================================== + ! Code to determine whether the canopy is closed, and which plants are either in the + ! understorey or overstorey. This is obviosuly far too complicated for it's own good + ! ===================================================================================== use FatesConstantsMod , only : r8 => fates_r8 use FatesGlobals , only : fates_log @@ -127,7 +127,7 @@ subroutine canopy_structure( currentSite , bc_in ) ! do while (associated(currentPatch)) ! Patch loop - + ! ------------------------------------------------------------------------------ ! Perform numerical checks on some cohort and patch structures ! ------------------------------------------------------------------------------ @@ -147,128 +147,140 @@ subroutine canopy_structure( currentSite , bc_in ) currentCohort => currentCohort%shorter enddo - if (currentPatch%area .gt. min_patch_area) then ! avoid numerical weirdness that shouldn't be happening anyway - ! Does any layer have excess area in it? Keep going until it does not... - patch_area_counter = 0 - area_not_balanced = .true. - - do while(area_not_balanced) + ! ------------------------------------------------------------------------------ + ! Check patch area to prevent numerical weirdness + ! ------------------------------------------------------------------------------ + + if (currentPatch%area .lt. min_patch_area) then + + write(fates_log(),*) 'An incredibly small patch exists that should' + write(fates_log(),*) 'had been fused or culled already' + write(fates_log(),*) 'currentPatch%area: ',currentPatch%area + write(fates_log(),*) 'min_patch_area: ',min_patch_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end if - ! --------------------------------------------------------------------------------------- - ! Demotion Phase: Identify upper layers that are too full, and demote them to - ! the layers below. - ! --------------------------------------------------------------------------------------- + ! Does any layer have excess area in it? Keep going until it does not... + patch_area_counter = 0 + area_not_balanced = .true. + + do while(area_not_balanced) + + ! --------------------------------------------------------------------------- + ! Demotion Phase: Identify upper layers that are too full, and demote them to + ! the layers below. + ! --------------------------------------------------------------------------- + + ! Calculate how many layers we have in this canopy + ! This also checks the understory to see if its crown + ! area is large enough to warrant a temporary sub-understory layer + z = NumPotentialCanopyLayers(currentPatch,include_substory=.true.) + + do i_lyr = 1,z ! Loop around the currently occupied canopy layers. + call DemoteFromLayer(currentSite, currentPatch, i_lyr) + end do + + ! Remove cohorts that are incredibly sparse + call terminate_cohorts(currentSite, currentPatch, 1) + + call fuse_cohorts(currentPatch, bc_in) + + ! Remove cohorts for various other reasons + call terminate_cohorts(currentSite, currentPatch, 2) - ! Calculate how many layers we have in this canopy - ! This also checks the understory to see if its crown - ! area is large enough to warrant a temporary sub-understory layer - z = NumPotentialCanopyLayers(currentPatch,include_substory=.true.) + + ! --------------------------------------------------------------------------------------- + ! Promotion Phase: Identify if any upper-layers are underful and layers below them + ! have cohorts that can be split and promoted to the layer above. + ! --------------------------------------------------------------------------------------- + + ! Re-calculate Number of layers without the false substory + z = NumPotentialCanopyLayers(currentPatch,include_substory=.false.) - do i_lyr = 1,z ! Loop around the currently occupied canopy layers. - call DemoteFromLayer(currentSite, currentPatch, i_lyr) + ! We only promote if we have at least two layers + if (z>1) then + + do i_lyr=1,z-1 + call PromoteIntoLayer(currentSite, currentPatch, i_lyr) end do - + ! Remove cohorts that are incredibly sparse call terminate_cohorts(currentSite, currentPatch, 1) - + call fuse_cohorts(currentPatch, bc_in) - + ! Remove cohorts for various other reasons call terminate_cohorts(currentSite, currentPatch, 2) - - - ! --------------------------------------------------------------------------------------- - ! Promotion Phase: Identify if any upper-layers are underful and layers below them - ! have cohorts that can be split and promoted to the layer above. - ! --------------------------------------------------------------------------------------- - - ! Re-calculate Number of layers without the false substory - z = NumPotentialCanopyLayers(currentPatch,include_substory=.false.) - - ! We only promote if we have at least two layers - if (z>1) then - - do i_lyr=1,z-1 - call PromoteIntoLayer(currentSite, currentPatch, i_lyr) - end do - - ! Remove cohorts that are incredibly sparse - call terminate_cohorts(currentSite, currentPatch, 1) - - call fuse_cohorts(currentPatch, bc_in) - - ! Remove cohorts for various other reasons - call terminate_cohorts(currentSite, currentPatch, 2) - - end if - - ! --------------------------------------------------------------------------------------- - ! Check on Layer Area (if the layer differences are not small - ! Continue trying to demote/promote. Its possible on the first pass through, - ! that cohort fusion has nudged the areas a little bit. - ! --------------------------------------------------------------------------------------- - - z = NumPotentialCanopyLayers(currentPatch,include_substory=.false.) - area_not_balanced = .false. - do i_lyr = 1,z - call CanopyLayerArea(currentPatch,i_lyr,arealayer(i_lyr)) - if( (arealayer(i_lyr)-currentPatch%area) > area_check_precision )then - area_not_balanced = .true. - endif - enddo - - ! --------------------------------------------------------------------------------------- - ! Gracefully exit if too many iterations have gone by - ! --------------------------------------------------------------------------------------- - - patch_area_counter = patch_area_counter + 1 - if(patch_area_counter > max_patch_iterations) then - write(fates_log(),*) 'PATCH AREA CHECK NOT CLOSING' - write(fates_log(),*) 'patch area:',currentpatch%area - write(fates_log(),*) 'lat:',currentpatch%siteptr%lat - write(fates_log(),*) 'lon:',currentpatch%siteptr%lon - currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - write(fates_log(),*) 'coh ilayer:',currentCohort%canopy_layer - write(fates_log(),*) 'coh dbh:',currentCohort%dbh - write(fates_log(),*) 'coh pft:',currentCohort%pft - write(fates_log(),*) 'coh n:',currentCohort%n - write(fates_log(),*) 'coh carea:',currentCohort%c_area - currentCohort => currentCohort%shorter - enddo - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - enddo ! do while(area_not_balanced) - - - ! Set current canopy layer occupancy indicator. - currentPatch%NCL_p = min(nclmax,z) - - ! ------------------------------------------------------------------------------------------- - ! if we are using "strict PPA", then calculate a z_star value as - ! the height of the smallest tree in the canopy - ! loop from top to bottom and locate the shortest cohort in level 1 whose shorter - ! neighbor is in level 2 set zstar as the ehight of that shortest level 1 cohort - ! ------------------------------------------------------------------------------------------- + + end if + + ! --------------------------------------------------------------------------------------- + ! Check on Layer Area (if the layer differences are not small + ! Continue trying to demote/promote. Its possible on the first pass through, + ! that cohort fusion has nudged the areas a little bit. + ! --------------------------------------------------------------------------------------- + + z = NumPotentialCanopyLayers(currentPatch,include_substory=.false.) + area_not_balanced = .false. + do i_lyr = 1,z + call CanopyLayerArea(currentPatch,i_lyr,arealayer(i_lyr)) + if( (arealayer(i_lyr)-currentPatch%area) > area_check_precision )then + area_not_balanced = .true. + endif + enddo + + ! --------------------------------------------------------------------------------------- + ! Gracefully exit if too many iterations have gone by + ! --------------------------------------------------------------------------------------- - if ( ED_val_comp_excln .lt. 0.0_r8) then - currentPatch%zstar = 0._r8 + patch_area_counter = patch_area_counter + 1 + if(patch_area_counter > max_patch_iterations) then + write(fates_log(),*) 'PATCH AREA CHECK NOT CLOSING' + write(fates_log(),*) 'patch area:',currentpatch%area + write(fates_log(),*) 'lat:',currentpatch%siteptr%lat + write(fates_log(),*) 'lon:',currentpatch%siteptr%lon currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - if(currentCohort%canopy_layer .eq. 2)then - if (associated(currentCohort%taller)) then - if (currentCohort%taller%canopy_layer .eq. 1 ) then - currentPatch%zstar = currentCohort%taller%hite - endif - endif - endif + do while (associated(currentCohort)) + write(fates_log(),*) 'coh ilayer:',currentCohort%canopy_layer + write(fates_log(),*) 'coh dbh:',currentCohort%dbh + write(fates_log(),*) 'coh pft:',currentCohort%pft + write(fates_log(),*) 'coh n:',currentCohort%n + write(fates_log(),*) 'coh carea:',currentCohort%c_area currentCohort => currentCohort%shorter enddo - endif - - end if !if (currentPatch%area .gt. min_patch_area) then + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + enddo ! do while(area_not_balanced) + + + ! Set current canopy layer occupancy indicator. + currentPatch%NCL_p = min(nclmax,z) + + ! ------------------------------------------------------------------------------------------- + ! if we are using "strict PPA", then calculate a z_star value as + ! the height of the smallest tree in the canopy + ! loop from top to bottom and locate the shortest cohort in level 1 whose shorter + ! neighbor is in level 2 set zstar as the ehight of that shortest level 1 cohort + ! ------------------------------------------------------------------------------------------- + + if ( ED_val_comp_excln .lt. 0.0_r8) then + currentPatch%zstar = 0._r8 + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer .eq. 2)then + if (associated(currentCohort%taller)) then + if (currentCohort%taller%canopy_layer .eq. 1 ) then + currentPatch%zstar = currentCohort%taller%hite + endif + endif + endif + currentCohort => currentCohort%shorter + enddo + endif + currentPatch => currentPatch%younger enddo !patch @@ -320,7 +332,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort => currentPatch%tallest do while (associated(currentCohort)) - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) + call carea_allom(currentCohort%dbh,currentCohort%n, & + currentSite%spread,currentCohort%pft,currentCohort%c_area) if(arealayer > currentPatch%area.and.currentCohort%canopy_layer == i_lyr)then if (ED_val_comp_excln .ge. 0.0_r8 ) then @@ -405,11 +418,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) !put the litter from the terminated cohorts into the fragmenting pools do i_cwd=1,ncwd - currentPatch%CWD_AG(i_cwd) = currentPatch%CWD_AG(i_cwd) + (currentCohort%bdead+currentCohort%bsw) * & + currentPatch%CWD_AG(i_cwd) = currentPatch%CWD_AG(i_cwd) + & + (currentCohort%bdead+currentCohort%bsw) * & EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) * & SF_val_CWD_frac(i_cwd)*currentCohort%n/currentPatch%area - currentPatch%CWD_BG(i_cwd) = currentPatch%CWD_BG(i_cwd) + (currentCohort%bdead+currentCohort%bsw) * & + currentPatch%CWD_BG(i_cwd) = currentPatch%CWD_BG(i_cwd) + & + (currentCohort%bdead+currentCohort%bsw) * & (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) * & SF_val_CWD_frac(i_cwd)*currentCohort%n/currentPatch%area !litter flux per m2. @@ -433,8 +448,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentSite%CWD_BG_diagnostic_input_carbonflux(i_cwd) = & currentSite%CWD_BG_diagnostic_input_carbonflux(i_cwd) & + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * & - SF_val_CWD_frac(i_cwd) * (1.0_r8 - EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) & - * hlm_days_per_year / AREA + SF_val_CWD_frac(i_cwd) * (1.0_r8 - & + EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) * hlm_days_per_year / AREA enddo currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) = & @@ -447,7 +462,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%n = 0.0_r8 currentCohort%c_area = 0._r8 else - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) + call carea_allom(currentCohort%dbh,currentCohort%n, & + currentSite%spread,currentCohort%pft,currentCohort%c_area) endif call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,copyc%c_area) @@ -929,7 +945,8 @@ subroutine canopy_summarization( nsites, sites, bc_in ) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) - call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,currentCohort%pft,currentCohort%c_area) + call carea_allom(currentCohort%dbh,currentCohort%n, & + sites(s)%spread,currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(currentCohort) canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area @@ -943,16 +960,19 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ! Check for erroneous zero values. if(currentCohort%dbh <= 0._r8 .or. currentCohort%n == 0._r8)then - write(fates_log(),*) 'ED: dbh or n is zero in canopy_summarization', & + write(fates_log(),*) 'FATES: dbh or n is zero in canopy_summarization', & currentCohort%dbh,currentCohort%n + call endrun(msg=errMsg(sourcefile, __LINE__)) endif if(currentCohort%pft == 0.or.currentCohort%canopy_trim <= 0._r8)then - write(fates_log(),*) 'ED: PFT or trim is zero in canopy_summarization', & + write(fates_log(),*) 'FATES: PFT or trim is zero in canopy_summarization', & currentCohort%pft,currentCohort%canopy_trim + call endrun(msg=errMsg(sourcefile, __LINE__)) endif if( (currentCohort%bsw + currentCohort%bl + currentCohort%br) <= 0._r8)then - write(fates_log(),*) 'ED: alive biomass is zero in canopy_summarization', & + write(fates_log(),*) 'FATES: alive biomass is zero in canopy_summarization', & currentCohort%bsw + currentCohort%bl + currentCohort%br + call endrun(msg=errMsg(sourcefile, __LINE__)) endif currentCohort => currentCohort%taller @@ -960,12 +980,12 @@ subroutine canopy_summarization( nsites, sites, bc_in ) enddo ! ends 'do while(associated(currentCohort)) if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then - write(fates_log(),*) 'ED: canopy area bigger than area', & + write(fates_log(),*) 'FATES: canopy area bigger than area', & currentPatch%total_canopy_area ,currentPatch%area + call endrun(msg=errMsg(sourcefile, __LINE__)) currentPatch%total_canopy_area = currentPatch%area endif - currentPatch => currentPatch%younger end do !patch loop @@ -1014,7 +1034,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) real(r8) :: max_chite ! top of cohort canopy (m) real(r8) :: lai ! summed lai for checking m2 m-2 real(r8) :: snow_depth_avg ! avg snow over whole site - integer :: NC ! number of cohorts, for bug fixing. !---------------------------------------------------------------------- @@ -1025,43 +1044,40 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! 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 ! ed patch + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) - - !Calculate tree and canopy areas. - currentPatch%canopy_area = 0._r8 - currentPatch%canopy_layer_lai(:) = 0._r8 - NC = 0 - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) - currentPatch%canopy_area = currentPatch%canopy_area + currentCohort%c_area - NC = NC+1 - currentCohort => currentCohort%taller - enddo - ! if plants take up all the tile, then so does the canopy. - currentPatch%canopy_area = min(currentPatch%canopy_area,currentPatch%area) - - !calculate tree lai and sai. - currentPatch%ncan(:,:) = 0 - currentPatch%nrad(:,:) = 0 - currentPatch%lai = 0._r8 + + ! -------------------------------------------------------------------------------- + ! Calculate tree and canopy areas. + ! calculate tree lai and sai. + ! -------------------------------------------------------------------------------- + + currentPatch%canopy_layer_tai(:) = 0._r8 + currentPatch%ncan(:,:) = 0 + currentPatch%nrad(:,:) = 0 + currentPatch%lai = 0._r8 + currentCohort => currentPatch%shortest do while(associated(currentCohort)) + + ft = currentCohort%pft + currentCohort%treelai = tree_lai(currentCohort) currentCohort%treesai = tree_sai(currentCohort) - currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%canopy_area - currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%canopy_area - !Calculate the LAI plus SAI in each canopy storey. + currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area + currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area + + ! Number of actual vegetation layers in this cohort's crown currentCohort%NV = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) - currentPatch%lai = currentPatch%lai +currentCohort%lai + currentPatch%ncan(currentCohort%canopy_layer,ft) = & + max(currentPatch%ncan(currentCohort%canopy_layer,ft),currentCohort%NV) + + currentPatch%lai = currentPatch%lai + currentCohort%lai do L = 1,nclmax-1 if(currentCohort%canopy_layer == L)then - currentPatch%canopy_layer_lai(L) = currentPatch%canopy_layer_lai(L) + currentCohort%lai + & + currentPatch%canopy_layer_tai(L) = currentPatch%canopy_layer_tai(L) + currentCohort%lai + & currentCohort%sai endif enddo diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 01ed393d7f..1191e41804 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1217,9 +1217,8 @@ subroutine zero_patch(cp_p) currentPatch%age = nan currentPatch%age_class = 1 currentPatch%area = nan - currentPatch%canopy_layer_lai(:) = nan + currentPatch%canopy_layer_tai(:) = nan currentPatch%total_canopy_area = nan - currentPatch%canopy_area = nan currentPatch%bare_frac_area = nan currentPatch%tlai_profile(:,:,:) = nan @@ -1300,7 +1299,7 @@ subroutine zero_patch(cp_p) currentPatch%burnt_frac_litter(:) = 0.0_r8 currentPatch%btran_ft(:) = 0.0_r8 - currentPatch%canopy_layer_lai(:) = 0.0_r8 + currentPatch%canopy_layer_tai(:) = 0.0_r8 currentPatch%seeds_in(:) = 0.0_r8 currentPatch%seed_decay(:) = 0.0_r8 diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 9ac33909df..648f384160 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -327,7 +327,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(cl==NCL_p)then !are we in the top canopy layer or a shaded layer? laican = 0._r8 else - laican = sum(currentPatch%canopy_layer_lai(cl+1:NCL_p)) + laican = sum(currentPatch%canopy_layer_tai(cl+1:NCL_p)) end if ! Loop over leaf-layers diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 492dcf6ac3..054d9aed2e 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -291,10 +291,11 @@ module EDTypesMod ! LEAF ORGANIZATION real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 - real(r8) :: canopy_layer_lai(nclmax) ! lai that is shading this canopy layer: m2/m2 + real(r8) :: canopy_layer_tai(nclmax) ! total area index of each canopy layer + ! used to determine attenuation of parameters during + ! photosynthesis m2 veg / m2 of canopy area (patch without bare ground) 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) :: canopy_area ! area that is covered by vegetation : m2 (is this different to total_canopy_area? real(r8) :: bare_frac_area ! bare soil in this patch expressed as a fraction of the total soil surface. real(r8) :: lai ! leaf area index of patch real(r8) :: zstar ! height of smallest canopy tree -- only meaningful in "strict PPA" mode @@ -677,7 +678,6 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%ncl_p = ',cpatch%ncl_p write(fates_log(),*) 'pa%total_canopy_area = ',cpatch%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',cpatch%total_tree_area - write(fates_log(),*) 'pa%canopy_area = ',cpatch%canopy_area write(fates_log(),*) 'pa%bare_frac_area = ',cpatch%bare_frac_area write(fates_log(),*) 'pa%lai = ',cpatch%lai write(fates_log(),*) 'pa%zstar = ',cpatch%zstar From 4a01c7d186421b37bfbf3bcb3872e04c01d18949 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Mar 2018 16:49:25 -0800 Subject: [PATCH 11/24] Forgot to remove one extra instance of %canopy_area --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 239f8dd31e..a40a52e96f 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1361,7 +1361,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) if(currentCohort%canopy_layer==1)then write(fates_log(), *) 'ED: cohorts',currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area,currentPatch%canopy_area + currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area endif From 19974649606c6cc9a28f54598bf9311b26ff9ff2 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Mar 2018 17:08:55 -0800 Subject: [PATCH 12/24] Changed indexing from L to cl in canopy lai profiling. --- biogeochem/EDCanopyStructureMod.F90 | 128 +++++++++++++++------------- 1 file changed, 68 insertions(+), 60 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a40a52e96f..cbe31f270d 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1020,8 +1020,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) real(r8) :: fleaf ! fraction of cohort incepting area that is leaves. integer :: ft ! Plant functional type index. integer :: iv ! Vertical leaf layer index - integer :: L ! Canopy layer index - integer :: p ! clm patch index + integer :: cl ! Canopy layer index real(r8) :: fraction_exposed ! how much of this layer is not covered by snow? real(r8) :: layer_top_hite ! notional top height of this canopy layer (m) real(r8) :: layer_bottom_hite ! notional bottom height of this canopy layer (m) @@ -1061,6 +1060,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) do while(associated(currentCohort)) ft = currentCohort%pft + cl = currentCohort%canopy_layer currentCohort%treelai = tree_lai(currentCohort) currentCohort%treesai = tree_sai(currentCohort) @@ -1070,17 +1070,12 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! Number of actual vegetation layers in this cohort's crown currentCohort%NV = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - currentPatch%ncan(currentCohort%canopy_layer,ft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,ft),currentCohort%NV) + currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) currentPatch%lai = currentPatch%lai + currentCohort%lai - do L = 1,nclmax-1 - if(currentCohort%canopy_layer == L)then - currentPatch%canopy_layer_tai(L) = currentPatch%canopy_layer_tai(L) + currentCohort%lai + & - currentCohort%sai - endif - enddo + currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & + currentCohort%lai + currentCohort%sai currentCohort => currentCohort%taller @@ -1088,8 +1083,13 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%nrad = currentPatch%ncan if(smooth_leaf_distribution == 1)then - ! we are going to ignore the concept of canopy layers, and put all of the leaf area into height banded bins. - ! using the same domains as we had before, except that CL always = 1 + + ! ----------------------------------------------------------------------------- + ! we are going to ignore the concept of canopy layers, and put all of the leaf + ! area into height banded bins. using the same domains as we had before, except + ! that CL always = 1 + ! ----------------------------------------------------------------------------- + currentPatch%tlai_profile = 0._r8 currentPatch%tsai_profile = 0._r8 currentPatch%elai_profile = 0._r8 @@ -1180,7 +1180,13 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) else ! smooth leaf distribution - !Go through all cohorts and add their leaf area and canopy area to the accumulators. + + ! ----------------------------------------------------------------------------- + ! Standard canopy layering model. + ! Go through all cohorts and add their leaf area + ! and canopy area to the accumulators. + + currentPatch%tlai_profile = 0._r8 currentPatch%tsai_profile = 0._r8 currentPatch%elai_profile = 0._r8 @@ -1192,8 +1198,10 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentPatch%shortest do while(associated(currentCohort)) - L = currentCohort%canopy_layer + ft = currentCohort%pft + cl = currentCohort%canopy_layer + !Calculate the number of layers of thickness dlai, including the last one. currentCohort%NV = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) !how much of each tree is stem area index? Assuming that there is @@ -1205,9 +1213,9 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort%treelai,currentCohort%treesai,currentCohort%dbh, & currentCohort%n,currentCohort%status_coh endif - currentPatch%ncan(L,ft) = max(currentPatch%ncan(L,ft),currentCohort%NV) - currentPatch%nrad(L,ft) = currentPatch%ncan(L,ft) !fudge - this needs to be altered for snow burial - if(currentCohort%NV > currentPatch%nrad(L,ft))then + currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) + currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) !fudge - this needs to be altered for snow burial + if(currentCohort%NV > currentPatch%nrad(cl,ft))then write(fates_log(), *) 'ED: issue with NV',currentCohort%NV,currentCohort%pft,currentCohort%canopy_layer endif @@ -1257,28 +1265,28 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) remainder = dinc_ed end if - currentPatch%tlai_profile(L,ft,iv) = currentPatch%tlai_profile(L,ft,iv) + & + currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area - currentPatch%elai_profile(L,ft,iv) = currentPatch%elai_profile(L,ft,iv) + & + 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(L,ft,iv) = currentPatch%tsai_profile(L,ft,iv) + & + 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(L,ft,iv) = currentPatch%esai_profile(L,ft,iv) + & + 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(L,ft,iv) = min(1.0_r8,currentPatch%canopy_area_profile(L,ft,iv) + & + currentPatch%canopy_area_profile(cl,ft,iv) = min(1.0_r8,currentPatch%canopy_area_profile(cl,ft,iv) + & currentCohort%c_area/currentPatch%total_canopy_area) - currentPatch%layer_height_profile(L,ft,iv) = currentPatch%layer_height_profile(L,ft,iv) + & + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. if ( DEBUG ) then write(fates_log(), *) 'calc snow 2', snow_depth_si , frac_sno_eff_si - write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(L,ft,iv) + write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(cl,ft,iv) write(fates_log(), *) 'leaf_area_profile 1229 ', currentPatch%elai_profile(1,ft,iv) end if @@ -1288,72 +1296,72 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo !cohort - do L = 1,currentPatch%NCL_p + do cl = 1,currentPatch%NCL_p do ft = 1,numpft - do iv = 1,currentPatch%nrad(L,ft) + do iv = 1,currentPatch%nrad(cl,ft) !account for total canopy area - if(currentPatch%canopy_area_profile(L,ft,iv) > tiny(currentPatch%canopy_area_profile(L,ft,iv)))then + if(currentPatch%canopy_area_profile(cl,ft,iv) > tiny(currentPatch%canopy_area_profile(cl,ft,iv)))then - currentPatch%tlai_profile(L,ft,iv) = currentPatch%tlai_profile(L,ft,iv) / & - currentPatch%canopy_area_profile(L,ft,iv) + currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - currentPatch%tsai_profile(L,ft,iv) = currentPatch%tsai_profile(L,ft,iv) / & - currentPatch%canopy_area_profile(L,ft,iv) + currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - currentPatch%elai_profile(L,ft,iv) = currentPatch%elai_profile(L,ft,iv) / & - currentPatch%canopy_area_profile(L,ft,iv) + currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - currentPatch%esai_profile(L,ft,iv) = currentPatch%esai_profile(L,ft,iv) / & - currentPatch%canopy_area_profile(L,ft,iv) + currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) end if - if(currentPatch%tlai_profile(L,ft,iv)>tiny(currentPatch%tlai_profile(L,ft,iv)))then - currentPatch%layer_height_profile(L,ft,iv) = currentPatch%layer_height_profile(L,ft,iv) & - /currentPatch%tlai_profile(L,ft,iv) + if(currentPatch%tlai_profile(cl,ft,iv)>tiny(currentPatch%tlai_profile(cl,ft,iv)))then + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & + /currentPatch%tlai_profile(cl,ft,iv) end if enddo - currentPatch%tlai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8 - currentPatch%tsai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8 - currentPatch%elai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8 - currentPatch%esai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8 + currentPatch%tlai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%tsai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%elai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%esai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 enddo enddo currentPatch%nrad = currentPatch%ncan - do L = 1,currentPatch%NCL_p + do cl = 1,currentPatch%NCL_p do ft = 1,numpft - if(currentPatch%nrad(L,ft) > 30)then + if(currentPatch%nrad(cl,ft) > 30)then write(fates_log(), *) 'ED: issue w/ nrad' endif - currentPatch%present(L,ft) = 0 - do iv = 1, currentPatch%nrad(L,ft); - if(currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then - currentPatch%present(L,ft) = 1 + currentPatch%present(cl,ft) = 0 + do iv = 1, currentPatch%nrad(cl,ft); + if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then + currentPatch%present(cl,ft) = 1 endif end do !iv enddo !ft - if ( L == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & + if ( cl == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & .and. currentPatch%NCL_p > 1 ) then write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) endif - if (L == 1 .and. currentPatch%NCL_p > 1 .and. & + if (cl == 1 .and. currentPatch%NCL_p > 1 .and. & abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999) then write(fates_log(), *) 'ED: not enough area in the top canopy', & - sum(currentPatch%canopy_area_profile(L,1:numpft,1)), & - currentPatch%canopy_area_profile(L,1:numpft,1) + sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & + currentPatch%canopy_area_profile(cl,1:numpft,1) endif - if(abs(sum(currentPatch%canopy_area_profile(L,1:numpft,1))) > 1.00001)then + if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then write(fates_log(), *) 'ED: canopy-area-profile wrong', & - sum(currentPatch%canopy_area_profile(L,1:numpft,1)), & - currentPatch%patchno, L - write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(L,1:numpft,1),currentPatch%patchno + sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & + currentPatch%patchno, cl + write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno currentCohort => currentPatch%shortest @@ -1372,11 +1380,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) endif enddo ! loop over L - do L = 1,currentPatch%NCL_p + do cl = 1,currentPatch%NCL_p do ft = 1,numpft - if(currentPatch%present(L,FT) > 1)then - write(fates_log(), *) 'ED: present issue',L,ft,currentPatch%present(L,FT) - currentPatch%present(L,ft) = 1 + if(currentPatch%present(cl,FT) > 1)then + write(fates_log(), *) 'ED: present issue',cl,ft,currentPatch%present(cl,FT) + currentPatch%present(cl,ft) = 1 endif enddo enddo From a65dbebedba272c9cd4e8b9ae7587e2901694626 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Mar 2018 17:35:52 -0800 Subject: [PATCH 13/24] Removed unnecessary variable patch%lai. --- biogeochem/EDCanopyStructureMod.F90 | 52 +++++++++++++++++++++-------- biogeochem/EDPatchDynamicsMod.F90 | 7 ++-- biogeophys/EDSurfaceAlbedoMod.F90 | 4 +-- main/EDTypesMod.F90 | 2 -- main/FatesHistoryInterfaceMod.F90 | 5 +-- 5 files changed, 47 insertions(+), 23 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index cbe31f270d..7214592444 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -930,7 +930,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 currentPatch%total_tree_area = 0.0_r8 - currentPatch%lai = 0.0_r8 canopy_leaf_area = 0.0_r8 !update cohort quantitie s @@ -999,9 +998,37 @@ end subroutine canopy_summarization ! ===================================================================================== subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) + + ! ----------------------------------------------------------------------------------- + ! This subroutine calculates how leaf and stem areas are distributed + ! in vertical and horizontal space. ! - ! !DESCRIPTION: + ! The following cohort level diagnostics are updated here: + ! + ! currentCohort%treelai ! LAI per unit crown area (m2/m2) + ! currentCohort%treesai ! SAI per unit crown area (m2/m2) + ! currentCohort%lai ! LAI per unit canopy area (m2/m2) + ! currentCohort%sai ! SAI per unit canopy area (m2/m2) + ! currentCohort%NV ! The number of discrete vegetation + ! ! layers needed to describe this crown + ! + ! The following patch level diagnostics are updated here: + ! + ! currentPatch%canopy_layer_tai(cl) ! TAI of each canopy layer + ! currentPatch%ncan(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 + ! ! CURRENTLY SAME AS NCAN + ! currentPatch%tlai_profile(cl,ft,iv) ! + ! currentPatch%elai_profile(cl,ft,iv) + ! currentPatch%tsai_profile(cl,ft,iv) + ! currentPatch%esai_profile(cl,ft,iv) + ! currentPatch%canopy_area_profile(cl,ft,iv) + ! currentPatch%layer_height_profile(cl,ft,iv) ! + ! ----------------------------------------------------------------------------------- + ! !USES: use EDtypesMod , only : area, dinc_ed, hitemax, n_hite_bins @@ -1026,6 +1053,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) real(r8) :: layer_bottom_hite ! notional bottom height of this canopy layer (m) integer :: smooth_leaf_distribution ! is the leaf distribution this option (1) or not (0) real(r8) :: frac_canopy(N_HITE_BINS) ! amount of canopy in each height class + real(r8) :: patch_lai ! LAI summed over the patch in m2/m2 of canopy area real(r8) :: minh(N_HITE_BINS) ! minimum height in height class (m) real(r8) :: maxh(N_HITE_BINS) ! maximum height in height class (m) real(r8) :: dh ! vertical detph of height class (m) @@ -1036,6 +1064,8 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) !---------------------------------------------------------------------- + + smooth_leaf_distribution = 0 ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft @@ -1054,7 +1084,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%canopy_layer_tai(:) = 0._r8 currentPatch%ncan(:,:) = 0 currentPatch%nrad(:,:) = 0 - currentPatch%lai = 0._r8 + patch_lai = 0._r8 currentCohort => currentPatch%shortest do while(associated(currentCohort)) @@ -1072,7 +1102,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) - currentPatch%lai = currentPatch%lai + currentCohort%lai + patch_lai = patch_lai + currentCohort%lai currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & currentCohort%lai + currentCohort%sai @@ -1162,20 +1192,16 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo !currentCohort - !check - currentPatch%lai = 0._r8 - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - currentPatch%lai = currentPatch%lai +currentCohort%lai - currentCohort => currentCohort%taller - enddo !currentCohort + ! ----------------------------------------------------------------------------- + ! Perform a leaf area conservation check on the LAI profile lai = 0.0_r8 do ft = 1,numpft lai = lai+ sum(currentPatch%tlai_profile(1,ft,:)) enddo - if(lai > currentPatch%lai)then - write(fates_log(), *) 'ED: problem with lai assignments' + if(lai > patch_lai)then + write(fates_log(), *) 'FATES: problem with lai assignments' + call endrun(msg=errMsg(sourcefile, __LINE__)) endif diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 1191e41804..6baf5bce9a 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1248,7 +1248,6 @@ subroutine zero_patch(cp_p) currentPatch%present(:,:) = 999 ! is there any of this pft in this layer? currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft - currentPatch%lai = nan ! leaf area index of patch currentPatch%pft_agb_profile(:,:) = nan ! DISTURBANCE @@ -1667,13 +1666,13 @@ subroutine terminate_patches(cs_pnt) ! This is only really meant for very old patches. if(associated(currentPatch%older) )then write(fates_log(),*) 'fusing to older patch because this one is too small',& - currentPatch%area, currentPatch%lai, & - currentPatch%older%area,currentPatch%older%lai + currentPatch%area, & + currentPatch%older%area call fuse_2_patches(currentPatch%older, currentPatch) write(fates_log(),*) 'after fusion to older patch',currentPatch%area else write(fates_log(),*) 'fusing to younger patch because oldest one is too small',& - currentPatch%area, currentPatch%lai + currentPatch%area tmpptr => currentPatch%younger call fuse_2_patches(currentPatch, currentPatch%younger) write(fates_log(),*) 'after fusion to younger patch' diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 7ee743cba1..ad085a31f2 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -799,7 +799,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) currentPatch%tr_soil_dir(ib)* & (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%lai + (1.0_r8-bc_in(s)%albgr_dir_rb(ib)) do ft =1,3 iv = currentPatch%nrad(1,ft) + 1 @@ -1003,7 +1003,7 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) if(bc_out(s)%fsun_pa(ifp) > 1._r8)then write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & - cpatch%lai,sunlai,shalai + sunlai,shalai endif elai = calc_areaindex(cpatch,'elai') diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 054d9aed2e..cf8089a469 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -297,7 +297,6 @@ module EDTypesMod 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) :: bare_frac_area ! bare soil in this patch expressed as a fraction of the total soil surface. - real(r8) :: lai ! leaf area index of patch real(r8) :: zstar ! height of smallest canopy tree -- only meaningful in "strict PPA" mode real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. m2/m2 @@ -679,7 +678,6 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%total_canopy_area = ',cpatch%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',cpatch%total_tree_area write(fates_log(),*) 'pa%bare_frac_area = ',cpatch%bare_frac_area - write(fates_log(),*) 'pa%lai = ',cpatch%lai write(fates_log(),*) 'pa%zstar = ',cpatch%zstar write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate write(fates_log(),*) '----------------------------------------' diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index d80a0cee5d..f869c2309c 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -1368,9 +1368,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! Increment some patch-age-resolved diagnostics hio_lai_si_age(io_si,cpatch%age_class) = hio_lai_si_age(io_si,cpatch%age_class) & - + cpatch%lai * cpatch%area + + sum(cpatch%tlai_profile(:,:,:)) * cpatch%area + hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & - + cpatch%ncl_p * cpatch%area + + cpatch%ncl_p * cpatch%area hio_npatches_si_age(io_si,cpatch%age_class) = hio_npatches_si_age(io_si,cpatch%age_class) + 1._r8 if ( ED_val_comp_excln .lt. 0._r8 ) then ! only valid when "strict ppa" enabled hio_zstar_si_age(io_si,cpatch%age_class) = hio_zstar_si_age(io_si,cpatch%age_class) & From e7dc96ffba6de572dfa908f31cd53c7e0f644491 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Mar 2018 18:50:54 -0800 Subject: [PATCH 14/24] more cleaning of EDCanopyStructure, trying to clarify normalization of lai profiles --- biogeochem/EDCanopyStructureMod.F90 | 385 +++++++++++++++------------- main/EDTypesMod.F90 | 19 +- 2 files changed, 220 insertions(+), 184 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 7214592444..1700fda549 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1020,11 +1020,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include ! ! layers occluded by snow ! ! CURRENTLY SAME AS NCAN - ! currentPatch%tlai_profile(cl,ft,iv) ! - ! currentPatch%elai_profile(cl,ft,iv) - ! currentPatch%tsai_profile(cl,ft,iv) - ! currentPatch%esai_profile(cl,ft,iv) - ! currentPatch%canopy_area_profile(cl,ft,iv) + ! currentPatch%tlai_profile(cl,ft,iv) ! m2 of leaves per m2 of canopy area + ! currentPatch%elai_profile(cl,ft,iv) ! non-snow covered m2 of leaves per m2 of canopy area + ! currentPatch%tsai_profile(cl,ft,iv) ! m2 of stems per m2 of canopy area + ! currentPatch%esai_profile(cl,ft,iv) ! non-snow covered m2 of stems per m2 of canopy area + ! currentPatch%canopy_area_profile(cl,ft,iv) ! currentPatch%layer_height_profile(cl,ft,iv) ! ! ----------------------------------------------------------------------------------- @@ -1099,7 +1099,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! Number of actual vegetation layers in this cohort's crown currentCohort%NV = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - + currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) patch_lai = patch_lai + currentCohort%lai @@ -1211,7 +1211,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! Standard canopy layering model. ! Go through all cohorts and add their leaf area ! and canopy area to the accumulators. - + ! ----------------------------------------------------------------------------- currentPatch%tlai_profile = 0._r8 currentPatch%tsai_profile = 0._r8 @@ -1219,210 +1219,237 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%esai_profile = 0._r8 currentPatch%layer_height_profile = 0._r8 currentPatch%canopy_area_profile(:,:,:) = 0._r8 - currentPatch%ncan(:,:) = 0 - currentPatch%nrad(:,:) = 0 - currentCohort => currentPatch%shortest + + ! ------------------------------------------------------------------------------ + ! It is remotely possible that in deserts we will not have any canopy + ! area, ie not plants at all... + ! ------------------------------------------------------------------------------ - do while(associated(currentCohort)) - - ft = currentCohort%pft - cl = currentCohort%canopy_layer - - !Calculate the number of layers of thickness dlai, including the last one. - currentCohort%NV = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - !how much of each tree is stem area index? Assuming that there is - if(currentCohort%treelai+currentCohort%treesai > 0._r8)then - fleaf = currentCohort%lai / (currentCohort%lai + currentCohort%sai) - else - fleaf = 0._r8 - write(fates_log(), *) 'ED: no stem or leaf area' ,currentCohort%pft,currentCohort%bl, & - currentCohort%treelai,currentCohort%treesai,currentCohort%dbh, & - currentCohort%n,currentCohort%status_coh - endif - currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) - currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) !fudge - this needs to be altered for snow burial - if(currentCohort%NV > currentPatch%nrad(cl,ft))then - write(fates_log(), *) 'ED: issue with NV',currentCohort%NV,currentCohort%pft,currentCohort%canopy_layer - endif + if (currentPatch%total_canopy_area < tiny(currentPatch%total_canopy_area)) then - !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. FIX(RF,032414)- for debugging jan 2012 - - do iv = 1,currentCohort%NV + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) - ! 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. - ! (RGK 03-01-2018 : we are not occulding any vegetation from snow right now) - - layer_top_hite = currentCohort%hite - & - ( dble(iv-1.0)/currentCohort%NV * currentCohort%hite * EDPftvarcon_inst%crown(currentCohort%pft) ) + ft = currentCohort%pft + cl = currentCohort%canopy_layer - layer_bottom_hite = currentCohort%hite - & - ( dble(iv)/currentCohort%NV * currentCohort%hite * EDPftvarcon_inst%crown(currentCohort%pft) ) + ! ---------------------------------------------------------------- + ! How much of each tree is stem area index? Assuming that there is + ! This may indeed be zero if there is a sensecent grass + ! ---------------------------------------------------------------- - fraction_exposed = 1.0_r8 - snow_depth_avg = snow_depth_si * frac_sno_eff_si - if(snow_depth_avg > layer_top_hite)then - fraction_exposed = 0._r8 + if( (currentCohort%treelai+currentCohort%treesai) > 0._r8)then + fleaf = currentCohort%lai / (currentCohort%lai + currentCohort%sai) + else + fleaf = 0._r8 endif - if(snow_depth_avg < layer_bottom_hite)then - fraction_exposed = 1._r8 + + ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX + ! SNOW BURIAL IS CURRENTLY TURNED OFF + ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX + + currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) + + ! -------------------------------------------------------------------------- + ! 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 - endif - if(snow_depth_avg>= layer_bottom_hite.and.snow_depth_avg <= layer_top_hite)then !only partly hidden... - fraction_exposed = max(0._r8,(min(1.0_r8,(snow_depth_avg-layer_bottom_hite)/ & + ! 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_hite = currentCohort%hite - & + ( dble(iv-1.0)/currentCohort%NV * currentCohort%hite * & + EDPftvarcon_inst%crown(currentCohort%pft) ) + + layer_bottom_hite = currentCohort%hite - & + ( dble(iv)/currentCohort%NV * currentCohort%hite * & + EDPftvarcon_inst%crown(currentCohort%pft) ) + + fraction_exposed = 1.0_r8 + snow_depth_avg = snow_depth_si * frac_sno_eff_si + if(snow_depth_avg > layer_top_hite)then + fraction_exposed = 0._r8 + endif + if(snow_depth_avg < layer_bottom_hite)then + fraction_exposed = 1._r8 + endif + if( snow_depth_avg>= layer_bottom_hite .and. & + snow_depth_avg <= layer_top_hite) then !only partly hidden... + fraction_exposed = max(0._r8,(min(1.0_r8,(snow_depth_avg-layer_bottom_hite)/ & (layer_top_hite-layer_bottom_hite )))) - endif - - ! =========== OVER-WRITE ================= - fraction_exposed= 1.0_r8 - ! =========== OVER-WRITE ================= - - if(iv==currentCohort%NV) then - remainder = (currentCohort%treelai + currentCohort%treesai) - (dinc_ed*dble(currentCohort%NV-1.0_r8)) - if(remainder > dinc_ed )then - write(fates_log(), *)'ED: issue with remainder',currentCohort%treelai,currentCohort%treesai,dinc_ed, & - currentCohort%NV,remainder - call endrun(msg=errMsg(sourcefile, __LINE__)) endif - else - remainder = dinc_ed - 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 + + ! =========== OVER-WRITE ================= + fraction_exposed= 1.0_r8 + ! =========== OVER-WRITE ================= + + if(iv==currentCohort%NV) then + remainder = (currentCohort%treelai + currentCohort%treesai) - & + (dinc_ed*dble(currentCohort%NV-1.0_r8)) + if(remainder > dinc_ed )then + write(fates_log(), *)'ED: issue with remainder', & + currentCohort%treelai,currentCohort%treesai,dinc_ed, & + currentCohort%NV,remainder + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + else + remainder = dinc_ed + 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 + + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & + (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & + (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. + + end do - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & - remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area + currentCohort => currentCohort%taller - 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 + enddo !cohort + + do cl = 1,currentPatch%NCL_p + + ! Check the canopy area profile + if( sum(currentPatch%canopy_area_profile(cl,:,:)) > 1.00000001_r8 ) then + write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' + write(fates_log(), *) 'cl: ',cl + write(fates_log(), *) 'ft: ',ft + write(fates_log(), *) 'iv: ',iv + write(fates_log(), *) 'cpatch%canopy_area_profile(cl,ft,iv): ', & + currentPatch%canopy_area_profile(cl,ft,iv) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - currentPatch%canopy_area_profile(cl,ft,iv) = min(1.0_r8,currentPatch%canopy_area_profile(cl,ft,iv) + & - currentCohort%c_area/currentPatch%total_canopy_area) - currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & - (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & - (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. - - if ( DEBUG ) then - write(fates_log(), *) 'calc snow 2', snow_depth_si , frac_sno_eff_si - write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(cl,ft,iv) - write(fates_log(), *) 'leaf_area_profile 1229 ', currentPatch%elai_profile(1,ft,iv) - end if + ! NOTES: IT SEEMS LIKE THESE FOLLOWING SHOULD BE m2 leaf / m2 of group footprint + ! IT SEEMS LIKE EACH LAYER-PFT SHOULD HAVE ITS OWN PARTIAL AREA - end do - - currentCohort => currentCohort%taller - - enddo !cohort - - do cl = 1,currentPatch%NCL_p - do ft = 1,numpft - do iv = 1,currentPatch%nrad(cl,ft) - !account for total canopy area - if(currentPatch%canopy_area_profile(cl,ft,iv) > tiny(currentPatch%canopy_area_profile(cl,ft,iv)))then - - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & - currentPatch%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) + + do ft = 1,numpft + do iv = 1,currentPatch%nrad(cl,ft) + + !account for total canopy area ? + if( currentPatch%canopy_area_profile(cl,ft,iv) > & + tiny(currentPatch%canopy_area_profile(cl,ft,iv)) )then + + currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & + currentPatch%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) + + currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & + currentPatch%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) + end if - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + if(currentPatch%tlai_profile(cl,ft,iv)>tiny(currentPatch%tlai_profile(cl,ft,iv)))then + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & + /currentPatch%tlai_profile(cl,ft,iv) + end if - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) - end if + enddo + + currentPatch%tlai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%tsai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%elai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 + currentPatch%esai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - if(currentPatch%tlai_profile(cl,ft,iv)>tiny(currentPatch%tlai_profile(cl,ft,iv)))then - currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & - /currentPatch%tlai_profile(cl,ft,iv) - end if - enddo - - currentPatch%tlai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%tsai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%elai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%esai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - enddo - enddo - - currentPatch%nrad = currentPatch%ncan - do cl = 1,currentPatch%NCL_p - do ft = 1,numpft - if(currentPatch%nrad(cl,ft) > 30)then - write(fates_log(), *) 'ED: issue w/ nrad' - endif - currentPatch%present(cl,ft) = 0 - do iv = 1, currentPatch%nrad(cl,ft); - if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%present(cl,ft) = 1 - endif - end do !iv - enddo !ft - if ( cl == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & - .and. currentPatch%NCL_p > 1 ) then - write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) - write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) - endif - - if (cl == 1 .and. currentPatch%NCL_p > 1 .and. & - abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999) then - write(fates_log(), *) 'ED: not enough area in the top canopy', & - sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & - currentPatch%canopy_area_profile(cl,1:numpft,1) - endif - - if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then - write(fates_log(), *) 'ED: canopy-area-profile wrong', & - sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & - currentPatch%patchno, cl - write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno + currentPatch%nrad = currentPatch%ncan + do cl = 1,currentPatch%NCL_p + do ft = 1,numpft + if(currentPatch%nrad(cl,ft) > 30)then + write(fates_log(), *) 'ED: issue w/ nrad' + endif + currentPatch%present(cl,ft) = 0 + do iv = 1, currentPatch%nrad(cl,ft); + if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then + currentPatch%present(cl,ft) = 1 + endif + end do !iv + enddo !ft - currentCohort => currentPatch%shortest + if ( cl == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & + .and. currentPatch%NCL_p > 1 ) then + write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) + write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) + endif - do while(associated(currentCohort)) - - if(currentCohort%canopy_layer==1)then - write(fates_log(), *) 'ED: cohorts',currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area - write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area - endif + if (cl == 1 .and. currentPatch%NCL_p > 1 .and. & + abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999) then + write(fates_log(), *) 'ED: not enough area in the top canopy', & + sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & + currentPatch%canopy_area_profile(cl,1:numpft,1) + endif + + if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then + write(fates_log(), *) 'ED: canopy-area-profile wrong', & + sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & + currentPatch%patchno, cl + write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno - currentCohort => currentCohort%taller + currentCohort => currentPatch%shortest - enddo !currentCohort - endif - enddo ! loop over L - - do cl = 1,currentPatch%NCL_p - do ft = 1,numpft - if(currentPatch%present(cl,FT) > 1)then - write(fates_log(), *) 'ED: present issue',cl,ft,currentPatch%present(cl,FT) - currentPatch%present(cl,ft) = 1 + do while(associated(currentCohort)) + + if(currentCohort%canopy_layer==1)then + write(fates_log(), *) 'ED: cohorts',currentCohort%dbh,currentCohort%c_area, & + currentPatch%total_canopy_area,currentPatch%area + write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & + currentCohort%c_area/currentPatch%total_canopy_area + endif + + currentCohort => currentCohort%taller + + enddo !currentCohort endif + enddo ! loop over cl + + do cl = 1,currentPatch%NCL_p + do ft = 1,numpft + if(currentPatch%present(cl,FT) > 1)then + write(fates_log(), *) 'ED: present issue',cl,ft,currentPatch%present(cl,FT) + currentPatch%present(cl,ft) = 1 + endif + enddo enddo - enddo + + endif !leaf distribution - endif !leaf distribution + end if currentPatch => currentPatch%younger enddo !patch - + return - end subroutine leaf_area_profile +end subroutine leaf_area_profile ! ====================================================================================== diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index cf8089a469..7c93a311d6 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -299,12 +299,21 @@ module EDTypesMod real(r8) :: bare_frac_area ! bare soil in this patch expressed as a fraction of the total soil surface. real(r8) :: zstar ! height of smallest canopy tree -- only meaningful in "strict PPA" mode - real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer. m2/m2 + + ! UNITS for the ai profiles + ! [ m2 leaf / m2 contributing crown footprints] + real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. + real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer + real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer + real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer + real(r8) :: layer_height_profile(nclmax,maxpft,nlevleaf) - real(r8) :: canopy_area_profile(nclmax,maxpft,nlevleaf) ! fraction of canopy in each canopy + 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 + + ! layer, pft, and leaf layer:- integer :: present(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 From 5704ed810e31fc127769bb671552c2ec74f02b3c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 12 Mar 2018 12:12:24 -0700 Subject: [PATCH 15/24] Light refactors in lai_profile calculations. --- biogeochem/EDCanopyStructureMod.F90 | 74 +++++++++++----------- biogeochem/EDPatchDynamicsMod.F90 | 2 +- biogeophys/EDSurfaceAlbedoMod.F90 | 24 +++---- biogeophys/FatesPlantRespPhotosynthMod.F90 | 14 ++-- main/EDTypesMod.F90 | 2 +- 5 files changed, 58 insertions(+), 58 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 1700fda549..dbac3a7ec0 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1110,7 +1110,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentCohort%taller enddo !currentCohort - currentPatch%nrad = currentPatch%ncan if(smooth_leaf_distribution == 1)then @@ -1246,10 +1245,23 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ! SNOW BURIAL IS CURRENTLY TURNED OFF + ! WHEN IT IS TURNED ON, IT WILL HAVE TO BE COMPARED + ! WITH SNOW HEIGHTS CALCULATED BELOW. ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 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 + + ! -------------------------------------------------------------------------- ! 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. @@ -1327,6 +1339,15 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentCohort%taller enddo !cohort + + ! -------------------------------------------------------------------------- + ! In the following loop we are now normalizing the effective and + ! total area profiles to convert from units of leaf/stem area per vegetated + ! canopy area, into leaf/stem area per area of their own radiative column + ! which is typically the footprint of all cohorts contained in the canopy + ! layer x pft bins. + ! Also perform some checks on area normalization. + ! -------------------------------------------------------------------------- do cl = 1,currentPatch%NCL_p @@ -1341,15 +1362,9 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - - ! NOTES: IT SEEMS LIKE THESE FOLLOWING SHOULD BE m2 leaf / m2 of group footprint - ! IT SEEMS LIKE EACH LAYER-PFT SHOULD HAVE ITS OWN PARTIAL AREA - - do ft = 1,numpft - do iv = 1,currentPatch%nrad(cl,ft) + do iv = 1,currentPatch%ncan(cl,ft) - !account for total canopy area ? if( currentPatch%canopy_area_profile(cl,ft,iv) > & tiny(currentPatch%canopy_area_profile(cl,ft,iv)) )then @@ -1381,16 +1396,20 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo enddo - currentPatch%nrad = currentPatch%ncan + ! -------------------------------------------------------------------------- + ! 1) Set the mask that identifies which PFT x can-layer combinations have + ! scattering elements in them. + ! 2) Run some final checks to see if canopy_area_profiles have summed up to + ! expected ranges + ! -------------------------------------------------------------------------- + do cl = 1,currentPatch%NCL_p + do ft = 1,numpft - if(currentPatch%nrad(cl,ft) > 30)then - write(fates_log(), *) 'ED: issue w/ nrad' - endif - currentPatch%present(cl,ft) = 0 + currentPatch%canopy_mask(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft); if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%present(cl,ft) = 1 + currentPatch%canopy_mask(cl,ft) = 1 endif end do !iv enddo !ft @@ -1399,13 +1418,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) .and. currentPatch%NCL_p > 1 ) then write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) - endif - - if (cl == 1 .and. currentPatch%NCL_p > 1 .and. & - abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999) then - write(fates_log(), *) 'ED: not enough area in the top canopy', & - sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & - currentPatch%canopy_area_profile(cl,1:numpft,1) + call endrun(msg=errMsg(sourcefile, __LINE__)) endif if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then @@ -1413,33 +1426,20 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & currentPatch%patchno, cl write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer==1)then write(fates_log(), *) 'ED: cohorts',currentCohort%dbh,currentCohort%c_area, & currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area endif - currentCohort => currentCohort%taller - enddo !currentCohort + call endrun(msg=errMsg(sourcefile, __LINE__)) endif enddo ! loop over cl - do cl = 1,currentPatch%NCL_p - do ft = 1,numpft - if(currentPatch%present(cl,FT) > 1)then - write(fates_log(), *) 'ED: present issue',cl,ft,currentPatch%present(cl,FT) - currentPatch%present(cl,ft) = 1 - endif - enddo - enddo - endif !leaf distribution end if @@ -1449,7 +1449,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo !patch return -end subroutine leaf_area_profile + end subroutine leaf_area_profile ! ====================================================================================== diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 6baf5bce9a..a467aa15b3 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1245,7 +1245,7 @@ subroutine zero_patch(cp_p) currentPatch%fabd(:) = nan ! fraction of incoming direct radiation that is absorbed by the canopy currentPatch%fabi(:) = nan ! fraction of incoming diffuse radiation that is absorbed by the canopy - currentPatch%present(:,:) = 999 ! is there any of this pft in this layer? + currentPatch%canopy_mask(:,:) = 999 ! is there any of this pft in this layer? currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft currentPatch%pft_agb_profile(:,:) = nan diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index ad085a31f2..396f54e008 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -197,10 +197,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) ! Is this pft/canopy layer combination present in this patch? do L = 1,nclmax do ft = 1,numpft - currentPatch%present(L,ft) = 0 + currentPatch%canopy_mask(L,ft) = 0 do iv = 1, currentPatch%nrad(L,ft) if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then - currentPatch%present(L,ft) = 1 + currentPatch%canopy_mask(L,ft) = 1 !I think 'present' is only used here... endif end do !iv @@ -260,7 +260,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8 !Each canopy layer (canopy, understorey) has multiple 'parallel' pft's do ft =1,numpft - if (currentPatch%present(L,ft) == 1)then !only do calculation if there are the appropriate leaves. + if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves. !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. ! Estimated do nine sky angles in increments of 10 degrees @@ -392,7 +392,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. do ft = 1,numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then !==============================================================================! ! Iterative solution do scattering !==============================================================================! @@ -440,7 +440,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) dif_ratio(L,ft,1,ib) * ftweight(L,ft,1) !instance where the first layer ftweight is used a proxy for the whole column. FTWA end do!hlm_numSWb - endif ! currentPatch%present + endif ! currentPatch%canopy_mask end do!ft end do!L @@ -450,7 +450,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = 1, currentPatch%NCL_p !work down from the top of the canopy. weighted_dif_down(L) = 0._r8 do ft = 1, numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! First estimates do downward and upward diffuse flux ! @@ -506,7 +506,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. weighted_dif_up(L) = 0._r8 do ft = 1, numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then !Bounce diffuse radiation off soil surface. iv = currentPatch%nrad(L,ft) + 1 if (L==currentPatch%NCL_p)then !is this the bottom layer ? @@ -562,7 +562,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = 1,currentPatch%NCL_p !working from the top down weighted_dif_down(L) = 0._r8 do ft =1,numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then ! forward diffuse flux within the canopy and at soil, working forward through canopy ! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy. ! Note: down = forward flux onto next layer @@ -618,7 +618,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = 1, currentPatch%NCL_p ! working from the top down. weighted_dif_up(L) = 0._r8 do ft =1,numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam) iv = currentPatch%nrad(L,ft) + 1 if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil @@ -670,7 +670,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) abs_dir_z(:,:) = 0._r8 abs_dif_z(:,:) = 0._r8 do ft =1,numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then !==============================================================================! ! Compute absorbed flux densities !==============================================================================! @@ -826,7 +826,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) lai_reduction(:) = 0.0_r8 do L = 1, currentPatch%NCL_p do ft =1,numpft - if (currentPatch%present(L,ft) == 1)then + if (currentPatch%canopy_mask(L,ft) == 1)then do iv = 1, currentPatch%nrad(L,ft) if (lai_change(L,ft,iv) > 0.0_r8)then lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv)) @@ -879,7 +879,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) write(fates_log(),*) 'bc_in(s)%albgr_dif_rb(ib)',bc_in(s)%albgr_dif_rb(ib) write(fates_log(),*) 'rhol',rhol(1:numpft,:) write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1) - write(fates_log(),*) 'present',currentPatch%present(1,1:numpft) + write(fates_log(),*) 'present',currentPatch%canopy_mask(1,1:numpft) write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1) bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 648f384160..35286333d4 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -242,7 +242,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! And then identify which layer/pft combinations have things in them. ! Output: ! currentPatch%ncan(:,:) - ! currentPatch%present(:,:) + ! currentPatch%canopy_mask(:,:) call UpdateCanopyNCanNRadPresent(currentPatch) @@ -322,7 +322,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! are there any leaves of this pft in this layer? - if(currentPatch%present(cl,ft) == 1)then + if(currentPatch%canopy_mask(cl,ft) == 1)then if(cl==NCL_p)then !are we in the top canopy layer or a shaded layer? laican = 0._r8 @@ -484,7 +484,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) currentCohort%gscan = 0.0_r8 currentCohort%ts_net_uptake(:) = 0.0_r8 - end if ! if(currentPatch%present(cl,ft) == 1)then + end if ! if(currentPatch%canopy_mask(cl,ft) == 1)then ! ------------------------------------------------------------------ @@ -1281,7 +1281,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! --------------------------------------------------------------------------------- ! This subroutine calculates two patch level quanities: ! currentPatch%ncan and - ! currentPatch%present + ! currentPatch%canopy_mask ! ! currentPatch%ncan(:,:) is a two dimensional array that indicates ! the total number of leaf layers (including those that are not exposed to light) @@ -1291,7 +1291,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! 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... ! - ! currentPatch%present(:,:) has the same dimensions, is binary, and + ! currentPatch%canopy_mask(:,:) has the same dimensions, is binary, and ! indicates whether or not leaf layers are present (by evaluating the canopy area ! profile). ! --------------------------------------------------------------------------------- @@ -1334,10 +1334,10 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! Now loop through and identify which layer and pft combo has scattering elements do cl = 1,nclmax do ft = 1,numpft - currentPatch%present(cl,ft) = 0 + currentPatch%canopy_mask(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft); if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%present(cl,ft) = 1 + currentPatch%canopy_mask(cl,ft) = 1 end if end do !iv enddo !ft diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 7c93a311d6..9f8b36d745 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -315,7 +315,7 @@ module EDTypesMod ! layer, pft, and leaf layer:- - integer :: present(nclmax,maxpft) ! is there any of this pft in this canopy layer? + 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 From 585c99bbc3920c869ada49708b096a7137c4ce6d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 12 Mar 2018 13:30:40 -0700 Subject: [PATCH 16/24] Minor syntax cleaning. --- biogeochem/EDCanopyStructureMod.F90 | 46 ++++++++++++----------------- 1 file changed, 19 insertions(+), 27 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index dbac3a7ec0..253de2f1be 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1020,12 +1020,13 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include ! ! layers occluded by snow ! ! CURRENTLY SAME AS NCAN - ! currentPatch%tlai_profile(cl,ft,iv) ! m2 of leaves per m2 of canopy area - ! currentPatch%elai_profile(cl,ft,iv) ! non-snow covered m2 of leaves per m2 of canopy area - ! currentPatch%tsai_profile(cl,ft,iv) ! m2 of stems per m2 of canopy area - ! currentPatch%esai_profile(cl,ft,iv) ! non-snow covered m2 of stems per m2 of canopy area - ! currentPatch%canopy_area_profile(cl,ft,iv) - ! currentPatch%layer_height_profile(cl,ft,iv) + ! currentPatch%tlai_profile(cl,ft,iv) ! m2 of leaves per m2 of the PFT's footprint + ! currentPatch%elai_profile(cl,ft,iv) ! non-snow covered m2 of leaves per m2 of PFT footprint + ! currentPatch%tsai_profile(cl,ft,iv) ! m2 of stems per m2 of PFT footprint + ! currentPatch%esai_profile(cl,ft,iv) ! non-snow covered m2 of stems per m2 of PFT footprint + ! currentPatch%canopy_area_profile(cl,ft,iv) ! Fractional area of leaf layer + ! ! relative to vegetated area + ! currentPatch%layer_height_profile(cl,ft,iv) ! Elevation of layer in m ! ! ----------------------------------------------------------------------------------- @@ -1081,10 +1082,16 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! calculate tree lai and sai. ! -------------------------------------------------------------------------------- - currentPatch%canopy_layer_tai(:) = 0._r8 - currentPatch%ncan(:,:) = 0 - currentPatch%nrad(:,:) = 0 - patch_lai = 0._r8 + currentPatch%canopy_layer_tai(:) = 0._r8 + currentPatch%ncan(:,:) = 0 + currentPatch%nrad(:,:) = 0 + patch_lai = 0._r8 + currentPatch%tlai_profile(:,:,:) = 0._r8 + currentPatch%tsai_profile(:,:,:) = 0._r8 + currentPatch%elai_profile(:,:,:) = 0._r8 + currentPatch%esai_profile(:,:,:) = 0._r8 + currentPatch%layer_height_profile(:,:,:) = 0._r8 + currentPatch%canopy_area_profile(:,:,:) = 0._r8 currentCohort => currentPatch%shortest do while(associated(currentCohort)) @@ -1118,11 +1125,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! area into height banded bins. using the same domains as we had before, except ! that CL always = 1 ! ----------------------------------------------------------------------------- - - currentPatch%tlai_profile = 0._r8 - currentPatch%tsai_profile = 0._r8 - currentPatch%elai_profile = 0._r8 - currentPatch%esai_profile = 0._r8 ! this is a crude way of dividing up the bins. Should it be a function of actual maximum height? dh = 1.0_r8*(HITEMAX/N_HITE_BINS) @@ -1212,12 +1214,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! and canopy area to the accumulators. ! ----------------------------------------------------------------------------- - currentPatch%tlai_profile = 0._r8 - currentPatch%tsai_profile = 0._r8 - currentPatch%elai_profile = 0._r8 - currentPatch%esai_profile = 0._r8 - currentPatch%layer_height_profile = 0._r8 - currentPatch%canopy_area_profile(:,:,:) = 0._r8 + ! ------------------------------------------------------------------------------ ! It is remotely possible that in deserts we will not have any canopy @@ -1388,11 +1385,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo - currentPatch%tlai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%tsai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%elai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - currentPatch%esai_profile(cl,ft,currentPatch%nrad(cl,ft)+1: nlevleaf) = 0._r8 - enddo enddo @@ -1407,7 +1399,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) do ft = 1,numpft currentPatch%canopy_mask(cl,ft) = 0 - do iv = 1, currentPatch%nrad(cl,ft); + do iv = 1, currentPatch%nrad(cl,ft) if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then currentPatch%canopy_mask(cl,ft) = 1 endif From 1fadfa4bff54ef8d52cb9e9a1d1cfe0feb57c3ce Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 12 Mar 2018 17:23:59 -0700 Subject: [PATCH 17/24] Bugfix on lai profile cleaning: a filter that checked if any canopy area existed was excluding when it should had been including. --- biogeochem/EDCanopyStructureMod.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 253de2f1be..58e63d0787 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1110,10 +1110,17 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) patch_lai = patch_lai + currentCohort%lai - + currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & currentCohort%lai + currentCohort%sai - + +! do cl = 1,nclmax-1 +! if(currentCohort%canopy_layer == cl)then +! currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & +! currentCohort%lai + currentCohort%sai +! endif +! enddo + currentCohort => currentCohort%taller enddo !currentCohort @@ -1214,14 +1221,12 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! and canopy area to the accumulators. ! ----------------------------------------------------------------------------- - - ! ------------------------------------------------------------------------------ ! It is remotely possible that in deserts we will not have any canopy ! area, ie not plants at all... ! ------------------------------------------------------------------------------ - if (currentPatch%total_canopy_area < tiny(currentPatch%total_canopy_area)) then + if (currentPatch%total_canopy_area > tiny(currentPatch%total_canopy_area)) then currentCohort => currentPatch%shortest do while(associated(currentCohort)) From 78cac81beb761e3198d12c946c3db8e9c39dc6f0 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 13 Mar 2018 15:04:42 -0700 Subject: [PATCH 18/24] Fixed a check on canopy area profile. It was summing over wrong indices and failing when it shouldnt. --- biogeochem/EDCanopyStructureMod.F90 | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 58e63d0787..09ea35ca56 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1354,16 +1354,17 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) do cl = 1,currentPatch%NCL_p ! Check the canopy area profile - if( sum(currentPatch%canopy_area_profile(cl,:,:)) > 1.00000001_r8 ) then - write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' - write(fates_log(), *) 'cl: ',cl - write(fates_log(), *) 'ft: ',ft - write(fates_log(), *) 'iv: ',iv - write(fates_log(), *) 'cpatch%canopy_area_profile(cl,ft,iv): ', & - currentPatch%canopy_area_profile(cl,ft,iv) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - + do iv = 1,currentPatch%ncan(cl,ft) + if( sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.00000001_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)) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + do ft = 1,numpft do iv = 1,currentPatch%ncan(cl,ft) From 64e0d799a77321bb49cd372b2ddbc4dc9bbeef7e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Mar 2018 11:35:52 -0700 Subject: [PATCH 19/24] Added some descriptive text and better zeroing for canopy_mask. --- biogeochem/EDCanopyStructureMod.F90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 09ea35ca56..8259a203ce 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1020,6 +1020,8 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include ! ! layers occluded by snow ! ! CURRENTLY SAME AS NCAN + ! currentPatch%canopy_mask(cl,ft) ! are there canopy elements in this pft-layer? + ! ! (This is redundant with nrad though...) ! currentPatch%tlai_profile(cl,ft,iv) ! m2 of leaves per m2 of the PFT's footprint ! currentPatch%elai_profile(cl,ft,iv) ! non-snow covered m2 of leaves per m2 of PFT footprint ! currentPatch%tsai_profile(cl,ft,iv) ! m2 of stems per m2 of PFT footprint @@ -1092,6 +1094,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%esai_profile(:,:,:) = 0._r8 currentPatch%layer_height_profile(:,:,:) = 0._r8 currentPatch%canopy_area_profile(:,:,:) = 0._r8 + currentPatch%canopy_mask(:,:) = 0 currentCohort => currentPatch%shortest do while(associated(currentCohort)) @@ -1404,7 +1407,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) do cl = 1,currentPatch%NCL_p do ft = 1,numpft - currentPatch%canopy_mask(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft) if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then currentPatch%canopy_mask(cl,ft) = 1 @@ -1414,20 +1416,20 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) if ( cl == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & .and. currentPatch%NCL_p > 1 ) then - write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) - write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) + write(fates_log(), *) 'FATES: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) + write(fates_log(), *) 'FATES: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) call endrun(msg=errMsg(sourcefile, __LINE__)) endif if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then - write(fates_log(), *) 'ED: canopy-area-profile wrong', & + write(fates_log(), *) 'FATES: canopy-area-profile wrong', & sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & currentPatch%patchno, cl - write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno + write(fates_log(), *) 'FATES: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==1)then - write(fates_log(), *) 'ED: cohorts',currentCohort%dbh,currentCohort%c_area, & + write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, & currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area From b81036d7f42917b3b4b12ea13c5916344ffb9370 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Mar 2018 16:00:22 -0700 Subject: [PATCH 20/24] relaxed run-fail on canopy area check against patch area. --- biogeochem/EDCanopyStructureMod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 8259a203ce..cc50a602ac 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -979,9 +979,11 @@ subroutine canopy_summarization( nsites, sites, bc_in ) enddo ! ends 'do while(associated(currentCohort)) if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then - write(fates_log(),*) 'FATES: canopy area bigger than area', & - currentPatch%total_canopy_area ,currentPatch%area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if ( currentPatch%total_canopy_area-currentPatch%area > 0.001_r8 ) then + write(fates_log(),*) 'FATES: canopy area bigger than area', & + currentPatch%total_canopy_area ,currentPatch%area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if currentPatch%total_canopy_area = currentPatch%area endif From db1ae64770dfaab2cdb40b59326b1de42af219e9 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Mar 2018 17:29:37 -0700 Subject: [PATCH 21/24] increased error tolerance on canopy layer area check. --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index cc50a602ac..70a1290705 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1360,7 +1360,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! Check the canopy area profile do iv = 1,currentPatch%ncan(cl,ft) - if( sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.00000001_r8 ) then + if( sum(currentPatch%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 From c4848a3656d0721b3eeda04a0a156526779244b1 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Mar 2018 19:09:57 -0700 Subject: [PATCH 22/24] For b4b comparability, reverting method of calculating canopy_layer_tai. --- biogeochem/EDCanopyStructureMod.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 70a1290705..c60b758737 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1116,15 +1116,15 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) patch_lai = patch_lai + currentCohort%lai - currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & - currentCohort%lai + currentCohort%sai - -! do cl = 1,nclmax-1 -! if(currentCohort%canopy_layer == cl)then -! currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & -! currentCohort%lai + currentCohort%sai -! endif -! enddo +! currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & +! currentCohort%lai + currentCohort%sai + + do cl = 1,nclmax-1 + if(currentCohort%canopy_layer == cl)then + currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & + currentCohort%lai + currentCohort%sai + endif + enddo currentCohort => currentCohort%taller From e89343cbcec4226cdc6f7a97c5f293a07678a1ae Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 22 Mar 2018 13:01:07 -0700 Subject: [PATCH 23/24] Improved some checks on canopy area, better logging, grouped together. --- biogeochem/EDCanopyStructureMod.F90 | 81 +++++++++++++++++------------ 1 file changed, 47 insertions(+), 34 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index c60b758737..38d0b64cc0 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1346,6 +1346,34 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentCohort%taller enddo !cohort + + ! -------------------------------------------------------------------------- + + ! If there is an upper-story, the top canopy layer + ! should have a value of exactly 1.0 in its top leaf layer + ! -------------------------------------------------------------------------- + + if ( (curentPatch%NLC_p > 1) .and. & + (sum(currentPatch%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 + 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 + write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & + currentCohort%c_area/currentPatch%total_canopy_area + endif + currentCohort => currentCohort%taller + enddo !currentCohort + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end if + ! -------------------------------------------------------------------------- ! In the following loop we are now normalizing the effective and @@ -1354,20 +1382,33 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! which is typically the footprint of all cohorts contained in the canopy ! layer x pft bins. ! Also perform some checks on area normalization. + ! Check the area of each leaf layer, across pfts. + ! It should never be larger than 1 or less than 0. ! -------------------------------------------------------------------------- - - do cl = 1,currentPatch%NCL_p - ! Check the canopy area profile + do cl = 1,currentPatch%NCL_p do iv = 1,currentPatch%ncan(cl,ft) + if( sum(currentPatch%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)) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + currentCohort => currentPatch%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 + write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & + currentCohort%c_area/currentPatch%total_canopy_area + endif + currentCohort => currentCohort%taller + enddo !currentCohort + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end do do ft = 1,numpft @@ -1400,14 +1441,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo ! -------------------------------------------------------------------------- - ! 1) Set the mask that identifies which PFT x can-layer combinations have + ! Set the mask that identifies which PFT x can-layer combinations have ! scattering elements in them. - ! 2) Run some final checks to see if canopy_area_profiles have summed up to - ! expected ranges ! -------------------------------------------------------------------------- do cl = 1,currentPatch%NCL_p - do ft = 1,numpft do iv = 1, currentPatch%nrad(cl,ft) if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then @@ -1415,31 +1453,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) endif end do !iv enddo !ft - - if ( cl == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & - .and. currentPatch%NCL_p > 1 ) then - write(fates_log(), *) 'FATES: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) - write(fates_log(), *) 'FATES: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - - if(abs(sum(currentPatch%canopy_area_profile(cl,1:numpft,1))) > 1.00001)then - write(fates_log(), *) 'FATES: canopy-area-profile wrong', & - sum(currentPatch%canopy_area_profile(cl,1:numpft,1)), & - currentPatch%patchno, cl - write(fates_log(), *) 'FATES: areas',currentPatch%canopy_area_profile(cl,1:numpft,1),currentPatch%patchno - currentCohort => currentPatch%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 - write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area - endif - currentCohort => currentCohort%taller - enddo !currentCohort - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif enddo ! loop over cl endif !leaf distribution From f5d5a7c9a3e7d7c81aa8f7b890c69aa3571b1725 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 22 Mar 2018 16:18:42 -0400 Subject: [PATCH 24/24] bug fix on logging and checks in canopy layer area calculations. --- biogeochem/EDCanopyStructureMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 38d0b64cc0..82121cae64 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1353,7 +1353,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! should have a value of exactly 1.0 in its top leaf layer ! -------------------------------------------------------------------------- - if ( (curentPatch%NLC_p > 1) .and. & + if ( (currentPatch%NCL_p > 1) .and. & (sum(currentPatch%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 @@ -1399,7 +1399,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==cl)then - write(fates_log(), *) 'FATES: cohorts in layer cl = 'cl, & + write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, & currentCohort%dbh,currentCohort%c_area, & currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, &