diff --git a/LICENSE.txt b/LICENSE.txt index 6c74023f..12f41fb6 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,34 +1,59 @@ -Copyright (c) 2013-2015, University Corporation for Atmospheric Research (UCAR) -All rights reserved. - -Developed by: - University Corporation for Atmospheric Research - National Center for Atmospheric Research - https://www2.cesm.ucar.edu/working-groups/sewg - -Permission is hereby granted, free of charge, to any person obtaining -a copy of this software and associated documentation files (the "Software"), -to deal with the Software without restriction, including without limitation -the rights to use, copy, modify, merge, publish, distribute, sublicense, -and/or sell copies of the Software, and to permit persons to whom -the Software is furnished to do so, subject to the following conditions: - - - Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimers. - - Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimers in the documentation - and/or other materials provided with the distribution. - - Neither the names of UCAR, or NCAR, - nor the names of its contributors may be used to endorse or promote - products derived from this Software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -POSSIBILITY OF SUCH DAMAGE. +Functionally Assembled Terrestrial Ecosystem Simulator (“FATES”) + +Copyright (c) 2016-2017, The Regents of the University of California, through Lawrence +Berkeley National Laboratory, University Corporation for Atmospheric Research, Los Alamos +National Security, LLC (LANS), as operator of Los Alamos National Laboratory (LANL), and +President and Fellows of Harvard College. All rights reserved. + +Copyright (c) 2013-2015, University Corporation for Atmospheric Research (UCAR). All +rights reserved. + +If you have questions about your rights to use or distribute this software, please contact +Berkeley Lab's Innovation & Partnerships Office at IPO@lbl.gov. + +NOTICE. This Software was developed under funding from the U.S. Department of Energy and +the U.S. Government consequently retains certain rights. As such, the U.S. Government has +been granted for itself and others acting on its behalf a paid-up, nonexclusive, +irrevocable, worldwide license in the Software to reproduce, distribute copies to the +public, prepare derivative works, and perform publicly and display publicly, and to permit +other to do so. + + Redistribution and use in source and binary forms, with or without modification, are + permitted provided that the following conditions are met: + + (1) Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + (2) Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + + (3) Neither the name of the University of California, Lawrence Berkeley National + Laboratory, University Corporation for Atmospheric Research, Los Alamos National + Security, LLC (LANS), as operator of Los Alamos National Laboratory (LANL), President and + Fellows of Harvard College, or the U.S. Dept. of Energy nor the names of its contributors + may be used to endorse or promote products derived from this software without specific + prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR +TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +You are under no obligation whatsoever to provide any bug fixes, patches, or upgrades to +the features, functionality or performance of the source code ("Enhancements") to anyone; +however, if you choose to make your Enhancements available either publicly, or directly to +Lawrence Berkeley National Laboratory, without imposing a separate written license +agreement for such Enhancements, then you hereby grant the following license to Lawrence +Berkeley National Laboratory, University Corporation for Atmospheric Research, Los Alamos +National Security, LLC (LANS), as operator of Los Alamos National Laboratory (LANL), +President and Fellows of Harvard College, and the U.S. Dept. of Energy: a non-exclusive, +royalty-free perpetual license to install, use, modify, prepare derivative works, +incorporate into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + diff --git a/README.md b/README.md index 50b3fc7e..f95e973e 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,22 @@ -# FATES - Functionally Assembled Terrestrial Ecosystem Simulator +# NGEE-T fates repository +------------------------------ -This is a READ-ONLY mirror of the NGEE-Tropics FATES public -release. Please contact the [NGEE-T](https://github.com/NGEET) group -for details, support and on-going development. +This is the developer repository of the Next Generation Ecosystem Experiment Tropics’ (NGEE-T) model: the Functionally Assembled Terrestrial Ecosystem Simulator (FATES). -FATES compsets are experimental and under active development. +For more information on the FATES model, see our wiki: https://github.com/NGEET/fates/wiki +## Important: +------------------------------ + +**Most users should not need to directly clone this repository. FATES needs to be run through a host model, and all supported host-models are in charge of cloning and loading the fates software.** + +FATES has support to be run via the Accelerated Climate Model for Energy (ACME) and the Community Earth System Model (CESM). + +https://climatemodeling.science.energy.gov/projects/accelerated-climate-modeling-energy + +http://www.cesm.ucar.edu/ + +The NGEE-T project maintains a mirror of CLM. That software system will automatically pull in the FATES software, and is where most users should go to clone the code: + +https://github.com/NGEET/fates-clm diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 old mode 100644 new mode 100755 index e4453a62..4b8cbd95 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -383,10 +383,15 @@ subroutine canopy_structure( currentSite , bc_in ) enddo currentPatch%ncl_p = min(z,nclmax) - enddo !is there still excess area in any layer? + enddo !is there still excess area in any layer? + + ! Remove cohorts that are incredibly sparse + call terminate_cohorts(currentSite, currentPatch, 1) call fuse_cohorts(currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch) + + ! Remove cohorts for various other reasons + call terminate_cohorts(currentSite, currentPatch, 2) ! ----------- Check cohort area ------------------------------! do i = 1,z @@ -605,8 +610,13 @@ subroutine canopy_structure( currentSite , bc_in ) endif enddo !is there still not enough canopy area in any layer? + ! remove cohorts that are extremely sparse + call terminate_cohorts(currentSite, currentPatch, 1) + call fuse_cohorts(currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch) + + ! remove cohorts for various other reasons + call terminate_cohorts(currentSite, currentPatch, 2) if(promswitch == 1)then !write(fates_log(),*) 'going into cohort check' @@ -1223,7 +1233,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) write(fates_log(), *) 'ED: canopy-area-profile wrong', & sum(currentPatch%canopy_area_profile(L,1:numpft_ed,1)), & currentPatch%patchno, L - write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(L,1:2,1),currentPatch%patchno + write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(L,1:numpft_ed,1),currentPatch%patchno currentCohort => currentPatch%shortest diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 old mode 100644 new mode 100755 index 0c6aaca1..46f297da --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -117,7 +117,6 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & call sizetype_class_index(new_cohort%dbh,new_cohort%pft, & new_cohort%size_class,new_cohort%size_by_pft_class) - if ( DEBUG ) write(fates_log(),*) 'EDCohortDyn I ',bstore ! This routine may be called during restarts, and at this point in the call sequence ! the actual cohort data is unknown, as this is really only used for allocation @@ -476,7 +475,7 @@ subroutine zero_cohort(cc_p) currentcohort%resp_acc_hold = 0._r8 currentcohort%carbon_balance = 0._r8 currentcohort%leaf_litter = 0._r8 - currentcohort%year_net_uptake(:) = 999 ! this needs to be 999, or trimming of new cohorts will break. + currentcohort%year_net_uptake(:) = 999._r8 ! this needs to be 999, or trimming of new cohorts will break. currentcohort%ts_net_uptake(:) = 0._r8 currentcohort%seed_prod = 0._r8 currentcohort%cfa = 0._r8 @@ -500,7 +499,7 @@ subroutine zero_cohort(cc_p) end subroutine zero_cohort !-------------------------------------------------------------------------------------! - subroutine terminate_cohorts( currentSite, patchptr ) + subroutine terminate_cohorts( currentSite, patchptr, level ) ! ! !DESCRIPTION: ! terminates cohorts when they get too small @@ -512,6 +511,16 @@ subroutine terminate_cohorts( currentSite, patchptr ) ! !ARGUMENTS type (ed_site_type) , intent(inout), target :: currentSite type (ed_patch_type), intent(inout), target :: patchptr + integer , intent(in) :: level + + ! Important point regarding termination levels. Termination is typically + ! called after fusion. We do this so that we can re-capture the biomass that would + ! otherwise be lost from termination. The biomass of a fused plant remains in the + ! live pool. However, some plant number densities can be so low that they + ! can cause numerical instabilities. Thus, we call terminate_cohorts at level=1 + ! before fusion to get rid of these cohorts that are so incredibly sparse, and then + ! terminate the remainder at level 2 for various other reasons. + ! ! !LOCAL VARIABLES: type (ed_patch_type) , pointer :: currentPatch @@ -529,16 +538,16 @@ subroutine terminate_cohorts( currentSite, patchptr ) nextc => currentCohort%shorter terminate = 0 - ! Check if number density is so low is breaks math - if (currentcohort%n < min_n_safemath) then + ! Check if number density is so low is breaks math (level 1) + if (currentcohort%n < min_n_safemath .and. level == 1) then terminate = 1 if ( DEBUG ) then write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh endif endif - ! The rest of these are only allowed if we are not dealing with a recruit - if (.not.currentCohort%isnew) then + ! The rest of these are only allowed if we are not dealing with a recruit (level 2) + if (.not.currentCohort%isnew .and. level == 2) then ! Not enough n or dbh if (currentCohort%n/currentPatch%area <= min_npm2 .or. & ! @@ -577,10 +586,10 @@ subroutine terminate_cohorts( currentSite, patchptr ) currentCohort%bstore, currentCohort%n endif - endif - endif + endif + endif ! if (.not.currentCohort%isnew .and. level == 2) then - if (terminate == 1) then + if (terminate == 1) then ! preserve a record of the to-be-terminated cohort for mortality accounting if (currentCohort%canopy_layer .eq. 1) then levcan = 1 @@ -649,244 +658,303 @@ end subroutine terminate_cohorts !-------------------------------------------------------------------------------------! subroutine fuse_cohorts(patchptr, bc_in) - ! - ! !DESCRIPTION: - ! Join similar cohorts to reduce total number - ! - ! !USES: - use EDTypesMod , only : nlevleaf - use EDParamsMod , only : ED_val_cohort_fusion_tol - ! - ! !ARGUMENTS - type (ed_patch_type), intent(inout), target :: patchptr - type (bc_in_type), intent(in) :: bc_in - ! - ! !LOCAL VARIABLES: - type (ed_patch_type) , pointer :: currentPatch - type (ed_cohort_type) , pointer :: currentCohort, nextc, nextnextc - integer :: i - integer :: fusion_took_place - integer :: maxcohorts !maximum total no of cohorts. Needs to be >numpft_edx2 - integer :: iterate !do we need to keep fusing to get below maxcohorts? - integer :: nocohorts - real(r8) :: newn - real(r8) :: diff - real(r8) :: dynamic_fusion_tolerance - !---------------------------------------------------------------------- - - !set initial fusion tolerance - dynamic_fusion_tolerance = ED_val_cohort_fusion_tol - - !This needs to be a function of the canopy layer, because otherwise, at canopy closure - !the number of cohorts doubles and very dissimilar cohorts are fused together - !because c_area and biomass are non-linear with dbh, this causes several mass inconsistancies - !in theory, all of this routine therefore causes minor losses of C and area, but these are below - !detection limit normally. - iterate = 1 - fusion_took_place = 0 - currentPatch => patchptr - maxcohorts = maxCohortsPerPatch - - !---------------------------------------------------------------------! - ! Keep doing this until nocohorts <= maxcohorts ! - !---------------------------------------------------------------------! - if (associated(currentPatch%shortest)) then - do while(iterate == 1) - - currentCohort => currentPatch%tallest - - ! The following logic continues the loop while the current cohort is not the shortest cohort - ! if they point to the same target (ie equivalence), then the loop ends. - ! This loop is different than the simple "continue while associated" loop in that - ! it omits the last cohort (because it has already been compared by that point) - - do while ( .not.associated(currentCohort,currentPatch%shortest) ) - - nextc => currentPatch%tallest - - do while (associated(nextc)) - nextnextc => nextc%shorter - diff = abs((currentCohort%dbh - nextc%dbh)/(0.5*(currentCohort%dbh + nextc%dbh))) - - !Criteria used to divide up the height continuum into different cohorts. - - if (diff < dynamic_fusion_tolerance) then - - ! Don't fuse a cohort with itself! - if (.not.associated(currentCohort,nextc) ) then - - if (currentCohort%pft == nextc%pft) then - - ! check cohorts in same c. layer. before fusing - - if (currentCohort%canopy_layer == nextc%canopy_layer) then - - ! Note: because newly recruited cohorts that have not experienced - ! a day yet will have un-known flux quantities or change rates - ! we don't want them fusing with non-new cohorts. We allow them - ! to fuse with other new cohorts to keep the total number of cohorts - ! down. - - if( .not.(currentCohort%isnew) .and. .not.(nextc%isnew) ) then - - newn = currentCohort%n + nextc%n - fusion_took_place = 1 - - - currentCohort%balive = (currentCohort%n*currentCohort%balive + nextc%n*nextc%balive)/newn - currentCohort%bdead = (currentCohort%n*currentCohort%bdead + nextc%n*nextc%bdead)/newn - - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn I ',currentCohort%bstore - - currentCohort%bstore = (currentCohort%n*currentCohort%bstore + nextc%n*nextc%bstore)/newn - - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn II ',currentCohort%bstore - - currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + nextc%n*nextc%seed_prod)/newn - currentCohort%root_md = (currentCohort%n*currentCohort%root_md + nextc%n*nextc%root_md)/newn - currentCohort%leaf_md = (currentCohort%n*currentCohort%leaf_md + nextc%n*nextc%leaf_md)/newn - currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory + nextc%n*nextc%laimemory)/newn - currentCohort%md = (currentCohort%n*currentCohort%md + nextc%n*nextc%md)/newn - - currentCohort%carbon_balance = (currentCohort%n*currentCohort%carbon_balance + & - nextc%n*nextc%carbon_balance)/newn - - currentCohort%storage_flux = (currentCohort%n*currentCohort%storage_flux + & - nextc%n*nextc%storage_flux)/newn - - currentCohort%b = (currentCohort%n*currentCohort%b + nextc%n*nextc%b)/newn - currentCohort%bsw = (currentCohort%n*currentCohort%bsw + nextc%n*nextc%bsw)/newn - currentCohort%bl = (currentCohort%n*currentCohort%bl + nextc%n*nextc%bl)/newn - - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn 569 ',currentCohort%br - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn 570 ',currentCohort%n - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn 571 ',nextc%br - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn 572 ',nextc%n - - currentCohort%br = (currentCohort%n*currentCohort%br + nextc%n*nextc%br)/newn - currentCohort%hite = (currentCohort%n*currentCohort%hite + nextc%n*nextc%hite)/newn - currentCohort%dbh = (currentCohort%n*currentCohort%dbh + nextc%n*nextc%dbh)/newn - - currentCohort%gpp_acc = (currentCohort%n*currentCohort%gpp_acc + nextc%n*nextc%gpp_acc)/newn - - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn III ',currentCohort%npp_acc - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn IV ',currentCohort%resp_acc - - currentCohort%npp_acc = (currentCohort%n*currentCohort%npp_acc + nextc%n*nextc%npp_acc)/newn - currentCohort%resp_acc = (currentCohort%n*currentCohort%resp_acc + nextc%n*nextc%resp_acc)/newn - - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn V ',currentCohort%npp_acc - if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn VI ',currentCohort%resp_acc - - currentCohort%resp_acc_hold = & - (currentCohort%n*currentCohort%resp_acc_hold + & - nextc%n*nextc%resp_acc_hold)/newn - currentCohort%npp_acc_hold = & - (currentCohort%n*currentCohort%npp_acc_hold + & - nextc%n*nextc%npp_acc_hold)/newn - currentCohort%gpp_acc_hold = & - (currentCohort%n*currentCohort%gpp_acc_hold + & - nextc%n*nextc%gpp_acc_hold)/newn - - currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim + nextc%n*nextc%canopy_trim)/newn - currentCohort%dmort = (currentCohort%n*currentCohort%dmort + nextc%n*nextc%dmort)/newn - currentCohort%fire_mort = (currentCohort%n*currentCohort%fire_mort + nextc%n*nextc%fire_mort)/newn - currentCohort%leaf_litter = (currentCohort%n*currentCohort%leaf_litter + nextc%n*nextc%leaf_litter)/newn - - ! mortality diagnostics - currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn - currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn - currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn - currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn - currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn - - ! npp diagnostics - currentCohort%npp_leaf = (currentCohort%n*currentCohort%npp_leaf + nextc%n*nextc%npp_leaf)/newn - currentCohort%npp_froot = (currentCohort%n*currentCohort%npp_froot + nextc%n*nextc%npp_froot)/newn - currentCohort%npp_bsw = (currentCohort%n*currentCohort%npp_bsw + nextc%n*nextc%npp_bsw)/newn - currentCohort%npp_bdead = (currentCohort%n*currentCohort%npp_bdead + nextc%n*nextc%npp_bdead)/newn - currentCohort%npp_bseed = (currentCohort%n*currentCohort%npp_bseed + nextc%n*nextc%npp_bseed)/newn - currentCohort%npp_store = (currentCohort%n*currentCohort%npp_store + nextc%n*nextc%npp_store)/newn - - ! recent canopy history - currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & - nextc%n*nextc%canopy_layer_yesterday)/newn - - do i=1, nlevleaf - if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then - currentCohort%year_net_uptake(i) = min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i)) - else - currentCohort%year_net_uptake(i) = (currentCohort%n*currentCohort%year_net_uptake(i) + & - nextc%n*nextc%year_net_uptake(i))/newn - endif - enddo - - if(use_fates_plant_hydro) call FuseCohortHydraulics(currentCohort,nextc,bc_in,newn) - - currentCohort%n = newn - !remove fused cohort from the list - nextc%taller%shorter => nextnextc - if (.not. associated(nextc%shorter)) then !this is the shortest cohort. - currentPatch%shortest => nextc%taller - else - nextnextc%taller => nextc%taller - endif - - if (associated(nextc)) then - if(use_fates_plant_hydro) call DeallocateHydrCohort(nextc) - deallocate(nextc) - endif - - endif ! Not a recruit - - endif !canopy layer - endif !pft - endif !index no. - endif !diff - - if (associated(nextc)) then - nextc => nextc%shorter - else - nextc => nextnextc !if we have removed next - endif - - enddo !end checking nextc cohort loop - - if (associated (currentCohort%shorter)) then - currentCohort => currentCohort%shorter - endif - enddo !end currentCohort cohort loop - - !---------------------------------------------------------------------! - ! Is the number of cohorts larger than the maximum? ! - !---------------------------------------------------------------------! - nocohorts = 0 - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - nocohorts = nocohorts + 1 - currentCohort => currentCohort%shorter - enddo - - if (nocohorts > maxcohorts) then - iterate = 1 - !---------------------------------------------------------------------! - ! Making profile tolerance larger means that more fusion will happen ! - !---------------------------------------------------------------------! - dynamic_fusion_tolerance = dynamic_fusion_tolerance * 1.1_r8 - - write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance + ! + ! !DESCRIPTION: + ! Join similar cohorts to reduce total number + ! + ! !USES: + use EDTypesMod , only : nlevleaf + use EDParamsMod , only : ED_val_cohort_fusion_tol + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) + ! + ! !ARGUMENTS + type (ed_patch_type), intent(inout), target :: patchptr + type (bc_in_type), intent(in) :: bc_in + ! + ! !LOCAL VARIABLES: + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort, nextc, nextnextc + integer :: i + integer :: fusion_took_place + integer :: maxcohorts !maximum total no of cohorts. Needs to be >numpft_edx2 + integer :: iterate !do we need to keep fusing to get below maxcohorts? + integer :: nocohorts + real(r8) :: newn + real(r8) :: diff + real(r8) :: dynamic_fusion_tolerance + + logical, parameter :: FUSE_DEBUG = .false. ! This debug is over-verbose + ! and gets its own flag + + !---------------------------------------------------------------------- + + !set initial fusion tolerance + dynamic_fusion_tolerance = ED_val_cohort_fusion_tol + + !This needs to be a function of the canopy layer, because otherwise, at canopy closure + !the number of cohorts doubles and very dissimilar cohorts are fused together + !because c_area and biomass are non-linear with dbh, this causes several mass inconsistancies + !in theory, all of this routine therefore causes minor losses of C and area, but these are below + !detection limit normally. + iterate = 1 + fusion_took_place = 0 + currentPatch => patchptr + maxcohorts = maxCohortsPerPatch + + !---------------------------------------------------------------------! + ! Keep doing this until nocohorts <= maxcohorts ! + !---------------------------------------------------------------------! + + if (associated(currentPatch%shortest)) then + do while(iterate == 1) + + currentCohort => currentPatch%tallest + + ! The following logic continues the loop while the current cohort is not the shortest cohort + ! if they point to the same target (ie equivalence), then the loop ends. + ! This loop is different than the simple "continue while associated" loop in that + ! it omits the last cohort (because it has already been compared by that point) + + do while ( .not.associated(currentCohort,currentPatch%shortest) ) + + nextc => currentPatch%tallest + + do while (associated(nextc)) + nextnextc => nextc%shorter + diff = abs((currentCohort%dbh - nextc%dbh)/(0.5*(currentCohort%dbh + nextc%dbh))) + + !Criteria used to divide up the height continuum into different cohorts. + + if (diff < dynamic_fusion_tolerance) then + + ! Don't fuse a cohort with itself! + if (.not.associated(currentCohort,nextc) ) then + + if (currentCohort%pft == nextc%pft) then + + ! check cohorts in same c. layer. before fusing + + if (currentCohort%canopy_layer == nextc%canopy_layer) then + + ! Note: because newly recruited cohorts that have not experienced + ! a day yet will have un-known flux quantities or change rates + ! we don't want them fusing with non-new cohorts. We allow them + ! to fuse with other new cohorts to keep the total number of cohorts + ! down. + + if( currentCohort%isnew.eqv.nextc%isnew ) then + + newn = currentCohort%n + nextc%n + fusion_took_place = 1 + + if ( FUSE_DEBUG .and. currentCohort%isnew ) then + write(fates_log(),*) 'Fusing Two Cohorts' + write(fates_log(),*) 'newn: ',newn + write(fates_log(),*) 'Cohort I, Cohort II' + write(fates_log(),*) 'n:',currentCohort%n,nextc%n + write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew + write(fates_log(),*) 'balive:',currentCohort%balive,nextc%balive + write(fates_log(),*) 'bdead:',currentCohort%bdead,nextc%bdead + write(fates_log(),*) 'bstore:',currentCohort%bstore,nextc%bstore + write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory + write(fates_log(),*) 'b:',currentCohort%b,nextc%b + write(fates_log(),*) 'bsw:',currentCohort%bsw,nextc%bsw + write(fates_log(),*) 'bl:',currentCohort%bl ,nextc%bl + write(fates_log(),*) 'br:',currentCohort%br,nextc%br + write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite + write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh + write(fates_log(),*) 'pft:',currentCohort%pft,nextc%pft + write(fates_log(),*) 'canopy_trim:',currentCohort%canopy_trim,nextc%canopy_trim + write(fates_log(),*) 'canopy_layer_yesterday:', & + currentCohort%canopy_layer_yesterday,nextc%canopy_layer_yesterday + do i=1, nlevleaf + write(fates_log(),*) 'leaf level: ',i,'year_net_uptake', & + currentCohort%year_net_uptake(i),nextc%year_net_uptake(i) + end do + end if + + currentCohort%balive = (currentCohort%n*currentCohort%balive & + + nextc%n*nextc%balive)/newn + currentCohort%bdead = (currentCohort%n*currentCohort%bdead & + + nextc%n*nextc%bdead)/newn + currentCohort%bstore = (currentCohort%n*currentCohort%bstore & + + nextc%n*nextc%bstore)/newn + currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & + + nextc%n*nextc%laimemory)/newn + currentCohort%b = (currentCohort%n*currentCohort%b & + + nextc%n*nextc%b)/newn + currentCohort%bsw = (currentCohort%n*currentCohort%bsw & + + nextc%n*nextc%bsw)/newn + currentCohort%bl = (currentCohort%n*currentCohort%bl & + + nextc%n*nextc%bl)/newn + currentCohort%br = (currentCohort%n*currentCohort%br & + + nextc%n*nextc%br)/newn + currentCohort%hite = (currentCohort%n*currentCohort%hite & + + nextc%n*nextc%hite)/newn + currentCohort%dbh = (currentCohort%n*currentCohort%dbh & + + nextc%n*nextc%dbh)/newn + currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + + nextc%n*nextc%canopy_trim)/newn + + call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & + currentCohort%size_class,currentCohort%size_by_pft_class) + + if(use_fates_plant_hydro) call FuseCohortHydraulics(currentCohort,nextc,bc_in,newn) + + ! recent canopy history + currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & + nextc%n*nextc%canopy_layer_yesterday)/newn + + ! Flux and biophysics variables have not been calculated for recruits we just default to + ! their initization values, which should be the same for eahc + + if ( .not.currentCohort%isnew) then + + currentCohort%md = (currentCohort%n*currentCohort%md + & + nextc%n*nextc%md)/newn + currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + & + nextc%n*nextc%seed_prod)/newn + currentCohort%root_md = (currentCohort%n*currentCohort%root_md + & + nextc%n*nextc%root_md)/newn + currentCohort%leaf_md = (currentCohort%n*currentCohort%leaf_md + & + nextc%n*nextc%leaf_md)/newn + currentCohort%carbon_balance = (currentCohort%n*currentCohort%carbon_balance + & + nextc%n*nextc%carbon_balance)/newn + currentCohort%storage_flux = (currentCohort%n*currentCohort%storage_flux + & + nextc%n*nextc%storage_flux)/newn + currentCohort%gpp_acc = (currentCohort%n*currentCohort%gpp_acc + & + nextc%n*nextc%gpp_acc)/newn + currentCohort%npp_acc = (currentCohort%n*currentCohort%npp_acc + & + nextc%n*nextc%npp_acc)/newn + currentCohort%resp_acc = (currentCohort%n*currentCohort%resp_acc + & + nextc%n*nextc%resp_acc)/newn + currentCohort%resp_acc_hold = & + (currentCohort%n*currentCohort%resp_acc_hold + & + nextc%n*nextc%resp_acc_hold)/newn + currentCohort%npp_acc_hold = & + (currentCohort%n*currentCohort%npp_acc_hold + & + nextc%n*nextc%npp_acc_hold)/newn + currentCohort%gpp_acc_hold = & + (currentCohort%n*currentCohort%gpp_acc_hold + & + nextc%n*nextc%gpp_acc_hold)/newn + + currentCohort%dmort = (currentCohort%n*currentCohort%dmort + & + nextc%n*nextc%dmort)/newn + currentCohort%fire_mort = (currentCohort%n*currentCohort%fire_mort + & + nextc%n*nextc%fire_mort)/newn + currentCohort%leaf_litter = (currentCohort%n*currentCohort%leaf_litter + & + nextc%n*nextc%leaf_litter)/newn + + ! mortality diagnostics + currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn + currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn + currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn + currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn + currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn + + ! npp diagnostics + currentCohort%npp_leaf = (currentCohort%n*currentCohort%npp_leaf + nextc%n*nextc%npp_leaf)/newn + currentCohort%npp_froot = (currentCohort%n*currentCohort%npp_froot + nextc%n*nextc%npp_froot)/newn + currentCohort%npp_bsw = (currentCohort%n*currentCohort%npp_bsw + nextc%n*nextc%npp_bsw)/newn + currentCohort%npp_bdead = (currentCohort%n*currentCohort%npp_bdead + nextc%n*nextc%npp_bdead)/newn + currentCohort%npp_bseed = (currentCohort%n*currentCohort%npp_bseed + nextc%n*nextc%npp_bseed)/newn + currentCohort%npp_store = (currentCohort%n*currentCohort%npp_store + nextc%n*nextc%npp_store)/newn + + ! biomass and dbh tendencies + currentCohort%ddbhdt = (currentCohort%n*currentCohort%ddbhdt + nextc%n*nextc%ddbhdt)/newn + currentCohort%dbalivedt = (currentCohort%n*currentCohort%dbalivedt + nextc%n*nextc%dbalivedt)/newn + currentCohort%dbdeaddt = (currentCohort%n*currentCohort%dbdeaddt + nextc%n*nextc%dbdeaddt)/newn + currentCohort%dbstoredt = (currentCohort%n*currentCohort%dbstoredt + nextc%n*nextc%dbstoredt)/newn + + do i=1, nlevleaf + if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then + currentCohort%year_net_uptake(i) = & + min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i)) + else + currentCohort%year_net_uptake(i) = (currentCohort%n*currentCohort%year_net_uptake(i) + & + nextc%n*nextc%year_net_uptake(i))/newn + endif + enddo + + end if !(currentCohort%isnew) + + currentCohort%n = newn + !remove fused cohort from the list + nextc%taller%shorter => nextnextc + if (.not. associated(nextc%shorter)) then !this is the shortest cohort. + currentPatch%shortest => nextc%taller + else + nextnextc%taller => nextc%taller + endif + + if (associated(nextc)) then + if(use_fates_plant_hydro) call DeallocateHydrCohort(nextc) + deallocate(nextc) + endif + + endif ! if( currentCohort%isnew.eqv.nextc%isnew ) then + + endif !canopy layer + endif !pft + endif !index no. + endif !diff + + if (associated(nextc)) then + nextc => nextc%shorter + else + nextc => nextnextc !if we have removed next + endif + + enddo !end checking nextc cohort loop + + if (associated (currentCohort%shorter)) then + currentCohort => currentCohort%shorter + endif + enddo !end currentCohort cohort loop + + !---------------------------------------------------------------------! + ! Is the number of cohorts larger than the maximum? ! + !---------------------------------------------------------------------! + nocohorts = 0 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + nocohorts = nocohorts + 1 + currentCohort => currentCohort%shorter + enddo + + if (nocohorts > maxcohorts) then + iterate = 1 + !---------------------------------------------------------------------! + ! Making profile tolerance larger means that more fusion will happen ! + !---------------------------------------------------------------------! + dynamic_fusion_tolerance = dynamic_fusion_tolerance * 1.1_r8 + + write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance + + else + iterate = 0 + endif - else - iterate = 0 - endif + if ( dynamic_fusion_tolerance .gt. 100._r8) then + ! something has gone terribly wrong and we need to report what + write(fates_log(),*) 'exceeded reasonable expectation of cohort fusion.' + currentCohort => currentPatch%tallest + nocohorts = 0 + do while(associated(currentCohort)) + write(fates_log(),*) 'cohort ', nocohorts, currentCohort%dbh, currentCohort%canopy_layer, currentCohort%n + nocohorts = nocohorts + 1 + currentCohort => currentCohort%shorter + enddo + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif - enddo !do while nocohorts>maxcohorts + enddo !do while nocohorts>maxcohorts - endif ! patch. + endif ! patch. - if (fusion_took_place == 1) then ! if fusion(s) occured sort cohorts - call sort_cohorts(currentPatch) - endif + if (fusion_took_place == 1) then ! if fusion(s) occured sort cohorts + call sort_cohorts(currentPatch) + endif end subroutine fuse_cohorts @@ -1099,6 +1167,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%status_coh = o%status_coh n%excl_weight = o%excl_weight n%prom_weight = o%prom_weight + n%size_class = o%size_class + n%size_by_pft_class = o%size_by_pft_class ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold @@ -1178,6 +1248,10 @@ subroutine copy_cohort( currentCohort,copyc ) if( use_fates_plant_hydro ) call CopyCohortHydraulics(n,o) + ! indices for binning + n%size_class = o%size_class + n%size_by_pft_class = o%size_by_pft_class + !Pointers n%taller => NULL() ! pointer to next tallest cohort n%shorter => NULL() ! pointer to next shorter cohort diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 old mode 100644 new mode 100755 index f0c081c8..e0abb086 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -225,11 +225,13 @@ real(r8) function c_area( cohort_in ) ! ============================================================================ use EDParamsMod , only : ED_val_grass_spread + use EDTypesMod , only : nclmax type(ed_cohort_type), intent(in) :: cohort_in real(r8) :: dbh ! Tree diameter at breat height. cm. real(r8) :: crown_area_to_dbh_exponent + integer :: can_layer_index ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth, ! but allowed to vary via the dbh2bl_dbh2carea_expnt_diff term (which has default value of zero) @@ -247,9 +249,22 @@ real(r8) function c_area( cohort_in ) end if dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft)) + + ! ---------------------------------------------------------------------------------- + ! The function c_area is called during the process of canopy position demotion + ! and promotion. As such, some cohorts are temporarily elevated to canopy positions + ! that are outside the number of alloted canopy spaces. Ie, a two story canopy + ! may have a third-story plant, if only for a moment. However, these plants + ! still need to generate a crown area to complete the promotion, demotion process. + ! So we allow layer index exceedence here and force it down to max. + ! (rgk/cdk 05/2017) + ! ---------------------------------------------------------------------------------- + + can_layer_index = min(cohort_in%canopy_layer,nclmax) + if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then c_area = 3.142_r8 * cohort_in%n * & - (cohort_in%patchptr%spread(cohort_in%canopy_layer)*dbh)**crown_area_to_dbh_exponent + (cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent else c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent end if diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 old mode 100644 new mode 100755 index f9896ba1..4bbf6194 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -23,6 +23,7 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : r8 => fates_r8 use FatesPlantHydraulicsMod, only : InitHydrCohort use FatesPlantHydraulicsMod, only : DeallocateHydrCohort + use EDParamsMod , only : fates_mortality_disturbance_fraction ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -86,7 +87,7 @@ subroutine disturbance_rates( site_in) currentCohort%patchptr => currentPatch call mortality_rates(currentCohort,cmort,hmort,bmort) - currentCohort%dmort = cmort+hmort+bmort + currentCohort%dmort = cmort+hmort+bmort currentCohort%c_area = c_area(currentCohort) ! Initialize diagnostic mortality rates @@ -99,6 +100,7 @@ subroutine disturbance_rates( site_in) if(currentCohort%canopy_layer == 1)then currentPatch%disturbance_rates(1) = currentPatch%disturbance_rates(1) + & + fates_mortality_disturbance_fraction * & min(1.0_r8,currentCohort%dmort)*hlm_freq_day*currentCohort%c_area/currentPatch%area endif @@ -289,7 +291,8 @@ subroutine spawn_patches( currentSite, bc_in) ! because this is the part of the original patch where no trees have actually fallen ! The diagnostic cmort,bmort and hmort rates have already been saved - currentCohort%n = currentCohort%n * (1.0_r8 - min(1.0_r8,currentCohort%dmort * hlm_freq_day)) + currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & + min(1.0_r8,currentCohort%dmort * hlm_freq_day)) nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. nc%cmort = nan ! The mortality diagnostics are set to nan because the cohort should dissappear @@ -401,10 +404,13 @@ subroutine spawn_patches( currentSite, bc_in) !update area of donor patch currentPatch%area = currentPatch%area - patch_site_areadis - !sort out the cohorts, since some of them may be so small as to need removing. - + ! sort out the cohorts, since some of them may be so small as to need removing. + ! the first call to terminate cohorts removes sparse number densities, + ! the second call removes for all other reasons (sparse culling must happen + ! before fusion) + call terminate_cohorts(currentSite, currentPatch, 1) call fuse_cohorts(currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch) + call terminate_cohorts(currentSite, currentPatch, 2) call sort_cohorts(currentPatch) currentPatch => currentPatch%younger @@ -420,8 +426,13 @@ subroutine spawn_patches( currentSite, bc_in) currentPatch%younger => new_patch currentSite%youngest_patch => new_patch + ! sort out the cohorts, since some of them may be so small as to need removing. + ! the first call to terminate cohorts removes sparse number densities, + ! the second call removes for all other reasons (sparse culling must happen + ! before fusion) + call terminate_cohorts(currentSite, new_patch, 1) call fuse_cohorts(new_patch, bc_in) - call terminate_cohorts(currentSite, new_patch) + call terminate_cohorts(currentSite, new_patch, 2) call sort_cohorts(new_patch) endif !end new_patch area @@ -514,6 +525,14 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis do p = 1,numpft_ed !move litter pool en mass into the new patch newPatch%root_litter(p) = newPatch%root_litter(p) + currentPatch%root_litter(p) * patch_site_areadis/newPatch%area newPatch%leaf_litter(p) = newPatch%leaf_litter(p) + currentPatch%leaf_litter(p) * patch_site_areadis/newPatch%area + + ! The fragmentation/decomposition flux from donor patches has already occured in existing patches. However + ! some of their area has been carved out for this new patches which is receiving donations. + ! Lets maintain conservation on that pre-existing mass flux in these newly disturbed patches + + newPatch%root_litter_out(p) = newPatch%root_litter_out(p) + currentPatch%root_litter_out(p) * patch_site_areadis/newPatch%area + newPatch%leaf_litter_out(p) = newPatch%leaf_litter_out(p) + currentPatch%leaf_litter_out(p) * patch_site_areadis/newPatch%area + enddo newPatch%spread = newPatch%spread + currentPatch%spread * patch_site_areadis/newPatch%area @@ -829,9 +848,10 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat new_patch%leaf_litter(p) = new_patch%leaf_litter(p) + canopy_mortality_leaf_litter(p) / litter_area * np_mult new_patch%root_litter(p) = new_patch%root_litter(p) + canopy_mortality_root_litter(p) / litter_area * np_mult + currentPatch%leaf_litter(p) = currentPatch%leaf_litter(p) + canopy_mortality_leaf_litter(p) / litter_area currentPatch%root_litter(p) = currentPatch%root_litter(p) + canopy_mortality_root_litter(p) / litter_area - + ! track as diagnostic fluxes currentSite%leaf_litter_diagnostic_input_carbonflux(p) = currentSite%leaf_litter_diagnostic_input_carbonflux(p) + & canopy_mortality_leaf_litter(p) * hlm_days_per_year / AREA @@ -924,9 +944,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ new_patch%leaf_litter_in(:) = 0._r8 new_patch%leaf_litter_out(:) = 0._r8 - new_patch%root_litter_in(:) = 0._r8 + new_patch%root_litter_in(:) = 0._r8 new_patch%root_litter_out(:) = 0._r8 - end subroutine create_patch @@ -1002,11 +1021,17 @@ subroutine zero_patch(cp_p) ! LITTER currentPatch%cwd_ag(:) = 0.0_r8 ! above ground coarse woody debris gc/m2. currentPatch%cwd_bg(:) = 0.0_r8 ! below ground coarse woody debris - currentPatch%root_litter(:) = 0.0_r8 - currentPatch%leaf_litter(:) = 0.0_r8 + currentPatch%root_litter(:) = 0.0_r8 ! In new disturbed patches, loops over donors to increment total, needs zero here + currentPatch%leaf_litter(:) = 0.0_r8 ! In new disturbed patches, loops over donors to increment total, needs zero here + + ! Cold-start initialized patches should have no litter flux in/out as they have not undergone any time. + ! Litter fluxes in/out also need to be initialized to zero for newly disturbed patches, as they + ! will incorporate the fluxes from donors over a loop, and need an initialization - currentPatch%leaf_litter_in(:) = 0.0_r8 - currentPatch%leaf_litter_out(:) = 0.0_r8 + currentPatch%leaf_litter_in(:) = 0.0_r8 ! As a newly created patch with no age, there is no flux in + currentPatch%leaf_litter_out(:) = 0.0_r8 ! As a newly created patch with no age, no frag or decomp has happened yet + currentPatch%root_litter_in(:) = 0.0_r8 ! As a newly created patch with no age, there is no flux in + currentPatch%root_litter_out(:) = 0.0_r8 ! As a newly created patch with no age, no frag or decomp has happened yet ! FIRE currentPatch%fuel_eff_moist = 0.0_r8 ! average fuel moisture content of the ground fuel @@ -1199,6 +1224,7 @@ subroutine fuse_patches( csite, bc_in ) end subroutine fuse_patches ! ============================================================================ + subroutine fuse_2_patches(dp, rp) ! ! !DESCRIPTION: @@ -1219,57 +1245,66 @@ subroutine fuse_2_patches(dp, rp) type (ed_cohort_type), pointer :: nextc ! Remembers next cohort in list type (ed_cohort_type), pointer :: storesmallcohort type (ed_cohort_type), pointer :: storebigcohort - integer :: c,p !counters for pft and litter size class. - integer :: tnull,snull ! are the tallest and shortest cohorts associated? - type(ed_patch_type), pointer :: youngerp ! pointer to the patch younger than donor - type(ed_patch_type), pointer :: olderp ! pointer to the patch older than donor - type(ed_site_type), pointer :: csite ! pointer to the donor patch's site - !--------------------------------------------------------------------- + integer :: c,p !counters for pft and litter size class. + integer :: tnull,snull ! are the tallest and shortest cohorts associated? + type(ed_patch_type), pointer :: youngerp ! pointer to the patch younger than donor + type(ed_patch_type), pointer :: olderp ! pointer to the patch older than donor + type(ed_site_type), pointer :: csite ! pointer to the donor patch's site + real(r8) :: inv_sum_area ! Inverse of the sum of the two patches areas + !----------------------------------------------------------------------------------------------- + + ! Generate a litany of area weighted averages + + inv_sum_area = 1.0_r8/(dp%area + rp%area) + + rp%age = (dp%age * dp%area + rp%age * rp%area) * inv_sum_area - !area weighted average of ages & litter - rp%age = (dp%age * dp%area + rp%age * rp%area)/(dp%area + rp%area) rp%age_class = get_age_class_index(rp%age) - do p = 1,numpft_ed - rp%seeds_in(p) = (rp%seeds_in(p)*rp%area + dp%seeds_in(p)*dp%area)/(rp%area + dp%area) - rp%seed_decay(p) = (rp%seed_decay(p)*rp%area + dp%seed_decay(p)*dp%area)/(rp%area + dp%area) - rp%seed_germination(p) = (rp%seed_germination(p)*rp%area + dp%seed_germination(p)*dp%area)/(rp%area + dp%area) - enddo do c = 1,ncwd - rp%cwd_ag(c) = (dp%cwd_ag(c)*dp%area + rp%cwd_ag(c)*rp%area)/(dp%area + rp%area) - rp%cwd_bg(c) = (dp%cwd_bg(c)*dp%area + rp%cwd_bg(c)*rp%area)/(dp%area + rp%area) + rp%cwd_ag(c) = (dp%cwd_ag(c)*dp%area + rp%cwd_ag(c)*rp%area) * inv_sum_area + rp%cwd_bg(c) = (dp%cwd_bg(c)*dp%area + rp%cwd_bg(c)*rp%area) * inv_sum_area enddo + + do p = 1,numpft_ed + rp%seeds_in(p) = (rp%seeds_in(p)*rp%area + dp%seeds_in(p)*dp%area) * inv_sum_area + rp%seed_decay(p) = (rp%seed_decay(p)*rp%area + dp%seed_decay(p)*dp%area) * inv_sum_area + rp%seed_germination(p) = (rp%seed_germination(p)*rp%area + dp%seed_germination(p)*dp%area) * inv_sum_area - do p = 1,numpft_ed - rp%leaf_litter(p) = (dp%leaf_litter(p)*dp%area + rp%leaf_litter(p)*rp%area)/(dp%area + rp%area) - rp%root_litter(p) = (dp%root_litter(p)*dp%area + rp%root_litter(p)*rp%area)/(dp%area + rp%area) + rp%leaf_litter(p) = (dp%leaf_litter(p)*dp%area + rp%leaf_litter(p)*rp%area) * inv_sum_area + rp%root_litter(p) = (dp%root_litter(p)*dp%area + rp%root_litter(p)*rp%area) * inv_sum_area + + rp%root_litter_out(p) = (dp%root_litter_out(p)*dp%area + rp%root_litter_out(p)*rp%area) * inv_sum_area + rp%leaf_litter_out(p) = (dp%leaf_litter_out(p)*dp%area + rp%leaf_litter_out(p)*rp%area) * inv_sum_area + + rp%root_litter_in(p) = (dp%root_litter_in(p)*dp%area + rp%root_litter_in(p)*rp%area) * inv_sum_area + rp%leaf_litter_in(p) = (dp%leaf_litter_in(p)*dp%area + rp%leaf_litter_in(p)*rp%area) * inv_sum_area + + rp%dleaf_litter_dt(p) = (dp%dleaf_litter_dt(p)*dp%area + rp%dleaf_litter_dt(p)*rp%area) * inv_sum_area + rp%droot_litter_dt(p) = (dp%droot_litter_dt(p)*dp%area + rp%droot_litter_dt(p)*rp%area) * inv_sum_area enddo - rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area)/(dp%area + rp%area) - rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area)/(dp%area + rp%area) - rp%sum_fuel = (dp%sum_fuel*dp%area + rp%sum_fuel*rp%area)/(dp%area + rp%area) - rp%fuel_bulkd = (dp%fuel_bulkd*dp%area + rp%fuel_bulkd*rp%area)/(dp%area + rp%area) - rp%fuel_sav = (dp%fuel_sav*dp%area + rp%fuel_sav*rp%area)/(dp%area + rp%area) - rp%fuel_mef = (dp%fuel_mef*dp%area + rp%fuel_mef*rp%area)/(dp%area + rp%area) - rp%ros_front = (dp%ros_front*dp%area + rp%ros_front*rp%area)/(dp%area + rp%area) - rp%effect_wspeed = (dp%effect_wspeed*dp%area + rp%effect_wspeed*rp%area)/(dp%area + rp%area) - rp%tau_l = (dp%tau_l*dp%area + rp%tau_l*rp%area)/(dp%area + rp%area) - rp%fuel_frac(:) = (dp%fuel_frac(:)*dp%area + rp%fuel_frac(:)*rp%area)/(dp%area + rp%area) - rp%tfc_ros = (dp%tfc_ros*dp%area + rp%tfc_ros*rp%area)/(dp%area + rp%area) - rp%fi = (dp%fi*dp%area + rp%fi*rp%area)/(dp%area + rp%area) - rp%fd = (dp%fd*dp%area + rp%fd*rp%area)/(dp%area + rp%area) - rp%ros_back = (dp%ros_back*dp%area + rp%ros_back*rp%area)/(dp%area + rp%area) - rp%ab = (dp%ab*dp%area + rp%ab*rp%area)/(dp%area + rp%area) - rp%nf = (dp%nf*dp%area + rp%nf*rp%area)/(dp%area + rp%area) - rp%sh = (dp%sh*dp%area + rp%sh*rp%area)/(dp%area + rp%area) - rp%frac_burnt = (dp%frac_burnt*dp%area + rp%frac_burnt*rp%area)/(dp%area + rp%area) - rp%burnt_frac_litter(:) = (dp%burnt_frac_litter(:)*dp%area + rp%burnt_frac_litter(:)*rp%area)/(dp%area + rp%area) - rp%btran_ft(:) = (dp%btran_ft(:)*dp%area + rp%btran_ft(:)*rp%area)/(dp%area + rp%area) - rp%dleaf_litter_dt(:) = (dp%dleaf_litter_dt(:)*dp%area + rp%dleaf_litter_dt(:)*rp%area)/(dp%area+rp%area) - rp%leaf_litter_in(:) = (dp%leaf_litter_in(:)*dp%area + rp%leaf_litter_in(:)*rp%area)/(dp%area+rp%area) - rp%leaf_litter_out(:) = (dp%leaf_litter_out(:)*dp%area + rp%leaf_litter_out(:)*rp%area)/(dp%area+rp%area) - + rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area + rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area + rp%sum_fuel = (dp%sum_fuel*dp%area + rp%sum_fuel*rp%area) * inv_sum_area + rp%fuel_bulkd = (dp%fuel_bulkd*dp%area + rp%fuel_bulkd*rp%area) * inv_sum_area + rp%fuel_sav = (dp%fuel_sav*dp%area + rp%fuel_sav*rp%area) * inv_sum_area + rp%fuel_mef = (dp%fuel_mef*dp%area + rp%fuel_mef*rp%area) * inv_sum_area + rp%ros_front = (dp%ros_front*dp%area + rp%ros_front*rp%area) * inv_sum_area + rp%effect_wspeed = (dp%effect_wspeed*dp%area + rp%effect_wspeed*rp%area) * inv_sum_area + rp%tau_l = (dp%tau_l*dp%area + rp%tau_l*rp%area) * inv_sum_area + rp%fuel_frac(:) = (dp%fuel_frac(:)*dp%area + rp%fuel_frac(:)*rp%area) * inv_sum_area + rp%tfc_ros = (dp%tfc_ros*dp%area + rp%tfc_ros*rp%area) * inv_sum_area + rp%fi = (dp%fi*dp%area + rp%fi*rp%area) * inv_sum_area + rp%fd = (dp%fd*dp%area + rp%fd*rp%area) * inv_sum_area + rp%ros_back = (dp%ros_back*dp%area + rp%ros_back*rp%area) * inv_sum_area + rp%ab = (dp%ab*dp%area + rp%ab*rp%area) * inv_sum_area + rp%nf = (dp%nf*dp%area + rp%nf*rp%area) * inv_sum_area + rp%sh = (dp%sh*dp%area + rp%sh*rp%area) * inv_sum_area + rp%frac_burnt = (dp%frac_burnt*dp%area + rp%frac_burnt*rp%area) * inv_sum_area + rp%burnt_frac_litter(:) = (dp%burnt_frac_litter(:)*dp%area + rp%burnt_frac_litter(:)*rp%area) * inv_sum_area + rp%btran_ft(:) = (dp%btran_ft(:)*dp%area + rp%btran_ft(:)*rp%area) * inv_sum_area rp%area = rp%area + dp%area !THIS MUST COME AT THE END! diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 old mode 100644 new mode 100755 index fefcc4cb..b1473722 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -31,7 +31,7 @@ module EDPhysiologyMod use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun - + use EDParamsMod , only : fates_mortality_disturbance_fraction implicit none private @@ -788,11 +788,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology + call mortality_rates(currentCohort,cmort,hmort,bmort) if (currentCohort%canopy_layer > 1)then - call mortality_rates(currentCohort,cmort,hmort,bmort) currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort) * currentCohort%n else - currentCohort%dndt = 0._r8 + currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) & + * (cmort+hmort+bmort) * currentCohort%n endif ! Height @@ -874,7 +875,6 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentCohort%npp_froot = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, & currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))) - if (Bleaf(currentCohort) > 0._r8)then if ( DEBUG ) write(fates_log(),*) 'EDphys A ',currentCohort%carbon_balance @@ -961,7 +961,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !FIX(RF,032414) - to fix high bl's. needed to prevent numerical errors without the ODEINT. if (currentCohort%balive > target_balive*1.1_r8)then va = 0.0_r8; vs = 1._r8 - write(fates_log(),*) 'using high bl cap',target_balive,currentCohort%balive + if (DEBUG) write(fates_log(),*) 'using high bl cap',target_balive,currentCohort%balive endif else @@ -1334,7 +1334,6 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use EDTypesMod, only : numpft_ed use FatesInterfaceMod, only : hlm_numlevdecomp_full use FatesInterfaceMod, only : hlm_numlevdecomp - use SoilBiogeochemVerticalProfileMod, only: surfprof_exp use EDPftvarcon, only : EDPftvarcon_inst use FatesConstantsMod, only : sec_per_day use EDParamsMod, only : ED_val_ag_biomass @@ -1345,11 +1344,8 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel - ! INTERF-TODO: remove the control parameters: exponential_rooting_profile, - ! pftspecific_rootingprofile, rootprof_exp, surfprof_exp - ! implicit none - ! + ! !ARGUMENTS integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) @@ -1383,6 +1379,10 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) real(r8), parameter :: rootprof_exp = 3. ! how steep profile is ! for root C inputs (1/ e-folding depth) (1/m) + ! NOTE(rgk, 201705) this parameter was brought over from SoilBiogeochemVerticalProfile + ! how steep profile is for surface components (1/ e_folding depth) (1/m) + real(r8), parameter :: surfprof_exp = 10. + ! NOTE(bja, 201608) as of clm4_5_10_r187 rootprof_beta is now a ! two dimensional array with the second dimension being water,1, ! or carbon,2,. These are currently hard coded, but may be diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 9ab4392c..6da2a280 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -830,11 +830,12 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) endif enddo enddo + if (lai_change(1,2,1).gt.0.0.and.lai_change(1,2,2).gt.0.0)then - ! write(fates_log(),*) 'lai_change(1,2,12)',lai_change(1,2,1:4) + ! write(fates_log(),*) 'lai_change(1,2,12)',lai_change(1,2,1:4) endif if (lai_change(1,2,2).gt.0.0.and.lai_change(1,2,3).gt.0.0)then - ! write(fates_log(),*) ' lai_change (1,2,23)',lai_change(1,2,1:4) + ! write(fates_log(),*) ' lai_change (1,2,23)',lai_change(1,2,1:4) endif if (lai_change(1,1,3).gt.0.0.and.lai_change(1,1,2).gt.0.0)then ! NO-OP @@ -846,7 +847,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) endif if (lai_change(1,1,4).gt.0.0.and.lai_change(1,1,5).gt.0.0)then ! NO-OP - ! write(fates_log(),*) 'first layer of lai_change 4 5',lai_change(1,1,1:5) + ! write(fates_log(),*) 'first layer of lai_change 4 5',lai_change(1,1,1:5) endif if (radtype == 1)then @@ -865,10 +866,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) write(fates_log(),*) 'Large Dir Radn consvn error',error ,ifp,ib write(fates_log(),*) 'diags', bc_out(s)%albd_parb(ifp,ib), bc_out(s)%ftdd_parb(ifp,ib), & bc_out(s)%ftid_parb(ifp,ib), bc_out(s)%fabd_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'ftweight',ftweight(1,1:2,1:4) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'ftweight',ftweight(1,1:numpft_ed,1:4) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'bc_in(s)%albgr_dir_rb(ib)',bc_in(s)%albgr_dir_rb(ib) @@ -884,16 +885,16 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) write(fates_log(),*) '>5% Dif Radn consvn error',error ,ifp,ib write(fates_log(),*) 'diags', bc_out(s)%albi_parb(ifp,ib), bc_out(s)%ftii_parb(ifp,ib), & bc_out(s)%fabi_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:2,1:4) - write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:2,1:4) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft_ed,1:4) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'bc_in(s)%albgr_dif_rb(ib)',bc_in(s)%albgr_dif_rb(ib) - write(fates_log(),*) 'rhol',rhol(1:2,:) - write(fates_log(),*) 'ftw',sum(ftweight(1,:,1)),ftweight(1,1:2,1) - write(fates_log(),*) 'present',currentPatch%present(1,1:2) - write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:2,1) + write(fates_log(),*) 'rhol',rhol(1:numpft_ed,:) + write(fates_log(),*) 'ftw',sum(ftweight(1,:,1)),ftweight(1,1:numpft_ed,1) + write(fates_log(),*) 'present',currentPatch%present(1,1:numpft_ed) + write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft_ed,1) bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error end if diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 old mode 100644 new mode 100755 index 5d0586a2..8a5ab941 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -1,4 +1,4 @@ -module SFMainMod + module SFMainMod ! ============================================================================ ! All subroutines realted to the SPITFIRE fire routine. @@ -12,6 +12,7 @@ module SFMainMod use FatesGlobals , only : fates_log use FatesInterfaceMod , only : bc_in_type + use EDPftvarcon , only : EDPftvarcon_inst use EDEcophysconType , only : EDecophyscon @@ -174,7 +175,7 @@ subroutine charecteristics_of_fuel ( currentSite ) ! There are SIX fuel classes ! 1) Leaf litter, 2:5) four CWD_AG pools (twig, s branch, l branch, trunk) and 6) live grass - ! NCWD =4 + ! NCWD =4 NFSC = 6 ! dl_sf = 1, tw_sf = 2, lb_sf = 4, tr_sf = 5, lg_sf = 6, ! zero fire arrays. @@ -304,7 +305,8 @@ subroutine wind_effect ( currentSite, bc_in) !*****************************************************************. ! Routine called daily from within ED within a site loop. - ! Calculates the effective windspeed based on vegetation charecteristics. + ! Calculates the effective windspeed based on vegetation charecteristics. + ! currentSite%wind is daily wind converted to m/min for Spitfire units use FatesConstantsMod, only : sec_per_min @@ -314,21 +316,21 @@ subroutine wind_effect ( currentSite, bc_in) type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - real(r8) :: wind ! daily wind in m/min - real(r8) :: total_grass_area ! per patch,in m2 - real(r8) :: tree_fraction ! site level. no units - real(r8) :: grass_fraction ! site level. no units - real(r8) :: bare_fraction ! site level. no units - integer :: iofp ! index of oldest fates patch + real(r8) :: total_grass_area ! per patch,in m2 + real(r8) :: tree_fraction ! site level. no units + real(r8) :: grass_fraction ! site level. no units + real(r8) :: bare_fraction ! site level. no units + integer :: iofp ! index of oldest fates patch + ! note - this is a patch level temperature, which probably won't have much inpact, ! unless we decide to ever calculated the NI for each patch. iofp = currentSite%oldest_patch%patchno - wind = bc_in%wind24_pa(iofp) * sec_per_min ! Convert to m/min for SPITFIRE units. + currentSite%wind = bc_in%wind24_pa(iofp) * sec_per_min !Convert to m/min for SPITFIRE if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'wind24', wind + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'wind24', currentSite%wind endif ! --- influence of wind speed, corrected for surface roughness---- ! --- averaged over the whole grid cell to prevent extreme divergence @@ -342,7 +344,7 @@ subroutine wind_effect ( currentSite, bc_in) currentCohort => currentPatch%tallest do while(associated(currentCohort)) - write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area + if (DEBUG) write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area else @@ -354,10 +356,10 @@ subroutine wind_effect ( currentSite, bc_in) grass_fraction = grass_fraction + min(currentPatch%area,total_grass_area)/AREA if(DEBUG)then - !write(fates_log(),*) 'SF currentPatch%area ',currentPatch%area - !write(fates_log(),*) 'SF currentPatch%total_area ',currentPatch%total_tree_area - !write(fates_log(),*) 'SF total_grass_area ',tree_fraction,grass_fraction - !write(fates_log(),*) 'SF AREA ',AREA + write(fates_log(),*) 'SF currentPatch%area ',currentPatch%area + write(fates_log(),*) 'SF currentPatch%total_area ',currentPatch%total_tree_area + write(fates_log(),*) 'SF total_grass_area ',tree_fraction,grass_fraction + write(fates_log(),*) 'SF AREA ',AREA endif currentPatch => currentPatch%younger @@ -376,7 +378,7 @@ subroutine wind_effect ( currentSite, bc_in) do while(associated(currentPatch)) currentPatch%total_tree_area = min(currentPatch%total_tree_area,currentPatch%area) ! effect_wspeed in units m/min - currentPatch%effect_wspeed = wind * (tree_fraction*0.4+(grass_fraction+bare_fraction)*0.6) + currentPatch%effect_wspeed = currentSite%wind * (tree_fraction*0.4+(grass_fraction+bare_fraction)*0.6) currentPatch => currentPatch%younger enddo !end patch loop @@ -391,32 +393,37 @@ subroutine rate_of_spread ( currentSite ) use SFParamsMod, only : SF_val_miner_total, SF_val_part_dens, & SF_val_miner_damp, SF_val_fuel_energy + use FatesInterfaceMod, only : hlm_current_day, hlm_current_month type(ed_site_type), intent(in), target :: currentSite type(ed_patch_type), pointer :: currentPatch - real(r8) dummy - ! Rothermal fire spread model parameters. - real(r8) beta,beta_op !weighted average of packing ratio (unitless) - real(r8) ir !reaction intensity (kJ/m2/min) - real(r8) xi,eps,phi_wind !all are unitless - real(r8) q_ig !heat of pre-ignition (kJ/kg) - real(r8) reaction_v_opt,reaction_v_max !reaction velocity (per min) - real(r8) moist_damp,mw_weight !moisture dampening coefficient and ratio fuel moisture to extinction - real(r8) bet !ratio of beta/beta_op - real(r8) a,b,c,e !function of fuel sav + real(r8) beta,beta_op ! weighted average of packing ratio (unitless) + real(r8) ir ! reaction intensity (kJ/m2/min) + real(r8) xi,eps,phi_wind ! all are unitless + real(r8) q_ig ! heat of pre-ignition (kJ/kg) + real(r8) reaction_v_opt,reaction_v_max !reaction velocity (per min)!optimum and maximum + real(r8) moist_damp,mw_weight ! moisture dampening coefficient and ratio fuel moisture to extinction + real(r8) beta_ratio ! ratio of beta/beta_op + real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation + real(r8) a,b,c,e ! function of fuel sav + real(r8),parameter::wind_max = 45.718_r8 !max wind speed (m/min)=150 ft/min per Lasslop etal 2014 + real(r8) wind_elev_fire !wind speed (m/min) at elevevation relevant for fire + + logical,parameter :: debug_windspeed = .false. !for debugging currentPatch=>currentSite%oldest_patch; do while(associated(currentPatch)) ! ---initialise parameters to zero.--- - bet = 0.0_r8; q_ig = 0.0_r8; eps = 0.0_r8; a = 0.0_r8; b = 0.0_r8; c = 0.0_r8; e = 0.0_r8 + beta_ratio = 0.0_r8; q_ig = 0.0_r8; eps = 0.0_r8; a = 0.0_r8; b = 0.0_r8; c = 0.0_r8; e = 0.0_r8 phi_wind = 0.0_r8; xi = 0.0_r8; reaction_v_max = 0.0_r8; reaction_v_opt = 0.0_r8; mw_weight = 0.0_r8 - moist_damp = 0.0_r8; ir = 0.0_r8; dummy = 0.0_r8; + moist_damp = 0.0_r8; ir = 0.0_r8; a_beta = 0.0_r8; currentPatch%ROS_front = 0.0_r8 + ! remove mineral content from net fuel load per Thonicke 2010 for ir calculation currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals @@ -429,7 +436,6 @@ subroutine rate_of_spread ( currentSite ) ! beta = packing ratio (unitless) ! fraction of fuel array volume occupied by fuel or compactness of fuel bed - beta = currentPatch%fuel_bulkd / SF_val_part_dens ! Equation A6 in Thonicke et al. 2010 @@ -438,7 +444,7 @@ subroutine rate_of_spread ( currentSite ) if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - beta ',beta if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - beta_op ',beta_op - bet = beta/beta_op !unitless + beta_ratio = beta/beta_op !unitless if(write_sf == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'esf ',currentPatch%fuel_eff_moist @@ -462,36 +468,47 @@ subroutine rate_of_spread ( currentSite ) if (DEBUG) then if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - c ',c - if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) & - 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed + if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - b ',b - if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - bet ',bet + if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - beta_ratio ',beta_ratio if ( hlm_masterproc == itrue .and.DEBUG) write(fates_log(),*) 'SF - e ',e endif ! Equation A5 in Thonicke et al. 2010 - ! convert effect_wspeed from m/min to ft/min for Rothermel ROS eqn ! phi_wind (unitless) - phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(bet**(-e)) + ! convert wind_elev_fire from m/min to ft/min for Rothermel ROS eqn + ! wind max per Lasslop et al 2014 to linearly reduce ROS for high wind speeds + !OLD! phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) + if (currentPatch%effect_wspeed .le. wind_max) then + wind_elev_fire = currentPatch%effect_wspeed + phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) + if (debug_windspeed) write(fates_log(),*) 'SF wind LESS max ', currentPatch%effect_wspeed + if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day + else + !max condition 225 ft/min (FIREMIP Rabin table A10 JSBACH-Spitfire) convert to 68.577 m/min + wind_elev_fire = max(0.0_r8,(68.577-0.5*currentPatch%effect_wspeed)) + phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) + if (debug_windspeed) write(fates_log(),*) 'SF wind GREATER max ', currentPatch%effect_wspeed + if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day + endif ! ---propagating flux---- - ! Equation A2 in Thonicke et al. - ! xi (unitless) - + ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972 + ! xi (unitless) xi = (exp((0.792_r8 + 3.7597_r8 * (currentPatch%fuel_sav**0.5_r8)) * (beta+0.1_r8))) / & (192_r8+7.9095_r8 * currentPatch%fuel_sav) ! ---reaction intensity---- ! Equation in table A1 Thonicke et al. 2010. a = 8.9033_r8 * (currentPatch%fuel_sav**(-0.7913_r8)) - dummy = exp(a*(1-bet)) + a_beta = exp(a*(1-beta_ratio)) !dummy variable for reaction_v_opt equation + ! Equation in table A1 Thonicke et al. 2010. ! reaction_v_max and reaction_v_opt = reaction velocity in units of per min - - ! Equation 36 in Rothermal 1972 12 + ! reaction_v_max = Equation 36 in Rothermal 1972 and Fig 12 reaction_v_max = 1.0_r8 / (0.0591_r8 + 2.926_r8* (currentPatch%fuel_sav**(-1.5_r8))) - ! Equation 38 in Rothermal 1972 and Fig 11 - reaction_v_opt = reaction_v_max*(bet**a)*dummy + ! reaction_v_opt = Equation 38 in Rothermal 1972 and Fig 11 + reaction_v_opt = reaction_v_max*(beta_ratio**a)*a_beta ! mw_weight = relative fuel moisture/fuel moisture of extinction ! average values for litter pools (dead leaves, twigs, small and large branches) plus grass @@ -508,7 +525,7 @@ subroutine rate_of_spread ( currentSite ) ! endif ! ir = reaction intenisty in kJ/m2/min - ! currentPatch%sum_fuel needs to be converted from kgC/m2 to kgBiomass/m2 for ir calculation + ! currentPatch%sum_fuel converted from kgC/m2 to kgBiomass/m2 for ir calculation ir = reaction_v_opt*(currentPatch%sum_fuel/0.45_r8)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp ! write(fates_log(),*) 'ir',gamma_aptr,moist_damp,SF_val_fuel_energy,SF_val_miner_damp @@ -523,7 +540,8 @@ subroutine rate_of_spread ( currentSite ) endif ! Equation 10 in Thonicke et al. 2010 ! backward ROS from Can FBP System (1992) in m/min - currentPatch%ROS_back = currentPatch%ROS_front*exp(-0.012_r8*currentPatch%effect_wspeed) + ! backward ROS wind not changed by vegetation + currentPatch%ROS_back = currentPatch%ROS_front*exp(-0.012_r8*currentSite%wind) currentPatch => currentPatch%younger @@ -620,6 +638,7 @@ subroutine fire_intensity ( currentSite ) !returns the updated currentPatch%FI value for each patch. !currentPatch%FI average fire intensity of flaming front during day. Backward ROS plays no role here. kJ/m/s or kW/m. + !currentSite%FDI probability that an ignition will start a fire !currentPatch%ROS_front forward ROS (m/min) !currentPatch%TFC_ROS total fuel consumed by flaming front (kgC/m2) @@ -627,13 +646,12 @@ subroutine fire_intensity ( currentSite ) use SFParamsMod, only : SF_val_fdi_alpha,SF_val_fuel_energy, & SF_val_max_durat, SF_val_durat_slope - type(ed_site_type), intent(in), target :: currentSite + type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), pointer :: currentPatch real(r8) ROS !m/s real(r8) W !kgBiomass/m2 - real(r8) :: d_fdi !change in the NI on this day to give fire duration. currentPatch => currentSite%oldest_patch; @@ -648,11 +666,15 @@ subroutine fire_intensity ( currentSite ) if (currentPatch%FI >= fire_threshold) then ! 50kW/m is the threshold for a self-sustaining fire currentPatch%fire = 1 ! Fire... :D - ! This is like but not identical to equation 7 in Thonicke et al. 2010. WHY? - d_FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI) !follows Venevsky et al GCB 2002 + ! Equation 7 from Venevsky et al GCB 2002 (modification of equation 8 in Thonicke et al. 2010) + ! FDI 0.1 = low, 0.3 moderate, 0.75 high, and 1 = extreme ignition potential for alpha 0.000337 + currentSite%FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI) ! Equation 14 in Thonicke et al. 2010 ! fire duration in minutes - currentPatch%FD = SF_val_max_durat / (1.0_r8 + SF_val_max_durat * exp(SF_val_durat_slope*d_FDI)) + + currentPatch%FD = (SF_val_max_durat+1) / (1.0_r8 + SF_val_max_durat * & + exp(SF_val_durat_slope*currentSite%FDI)) + if(write_SF == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fire duration minutes',currentPatch%fd endif @@ -730,40 +752,44 @@ subroutine area_burnt ( currentSite ) ! THIS SHOULD HAVE THE COLUMN AND LU AREA WEIGHT ALSO, NO? gridarea = km2_to_m2 ! 1M m2 in a km2 - currentPatch%NF = ED_val_nignitions * currentPatch%area/area /365 + !NF = number of lighting strikes per day per km2 + currentPatch%NF = ED_val_nignitions * currentPatch%area/area /365 ! If there are 15 lightening strickes per year, per km2. (approx from NASA product) ! then there are 15/365 s/km2 each day. ! Equation 1 in Thonicke et al. 2010 ! To Do: Connect here with the Li & Levis GDP fire suppression algorithm. - ! Equation 16 in arora and boer model. + ! Equation 16 in arora and boer model JGR 2005 !currentPatch%AB = currentPatch%AB *3.0_r8 + + !size of fire = equation 14 Arora and Boer JGR 2005 size_of_fire = ((3.1416_r8/(4.0_r8*lb))*((df+db)**2.0_r8)) - !AB is daily area burnt = size of fires in m2 * number of ignitions - currentPatch%AB = size_of_fire * currentPatch%NF + !AB = daily area burnt = size fires in m2 * num ignitions * prob ignition starts fire + currentPatch%AB = size_of_fire * currentPatch%NF * currentSite%FDI patch_area_in_m2 = gridarea*currentPatch%area/area - if (currentPatch%AB > patch_area_in_m2 ) then !all of patch burnt. - - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno, & - currentPatch%area/area,currentPatch%ab,patch_area_in_m2 - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, & - currentPatch%NF,currentPatch%FI,size_of_fire - - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'litter', & - currentPatch%sum_fuel,currentPatch%CWD_AG,currentPatch%leaf_litter - ! turn km2 into m2. work out total area burnt. - currentPatch%AB = patch_area_in_m2 - endif + currentPatch%frac_burnt = currentPatch%AB / patch_area_in_m2 if(write_SF == itrue)then if ( hlm_masterproc == itrue ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt endif + + if (currentPatch%frac_burnt > 1 ) then !all of patch burnt. + + currentPatch%frac_burnt = 1.0 ! capping at 1 same as %AB/patch_area_in_m2 + + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, & + currentPatch%NF,currentPatch%FI,size_of_fire + + endif + + endif endif! fire - currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt + currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt * currentPatch%area/area currentPatch => currentPatch%younger @@ -864,8 +890,8 @@ subroutine crown_damage ( currentSite ) if ((currentCohort%hite > 0.0_r8).and.(currentPatch%SH >= & (currentCohort%hite-currentCohort%hite*EDecophyscon%crown(currentCohort%pft)))) then - currentCohort%cfa = (currentPatch%SH-currentCohort%hite* & - EDecophyscon%crown(currentCohort%pft))/(currentCohort%hite-currentCohort%hite* & + currentCohort%cfa = (currentPatch%SH-currentCohort%hite*(1- & + EDecophyscon%crown(currentCohort%pft)))/(currentCohort%hite* & EDecophyscon%crown(currentCohort%pft)) else diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 old mode 100644 new mode 100755 index f7ada2f3..709c99c5 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -6,6 +6,7 @@ module EDInitMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : ifalse + use FatesConstantsMod , only : itrue use FatesGlobals , only : endrun => fates_endrun use EDTypesMod , only : nclmax use FatesGlobals , only : fates_log @@ -21,12 +22,20 @@ module EDInitMod use EDTypesMod , only : numpft_ed use FatesInterfaceMod , only : bc_in_type use EDTypesMod , only : use_fates_plant_hydro - + + ! CIME GLOBALS + use shr_log_mod , only : errMsg => shr_log_errMsg + implicit none private logical :: DEBUG = .false. + integer, parameter :: do_inv_init = ifalse + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + public :: zero_site public :: init_patches public :: set_site_properties @@ -184,65 +193,92 @@ end subroutine set_site_properties ! ============================================================================ subroutine init_patches( nsites, sites, bc_in) - ! - ! !DESCRIPTION: - !initialize patches on new ground - ! - ! !USES: - use EDParamsMod , only : ED_val_maxspread - use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps - ! - ! !ARGUMENTS - integer, intent(in) :: nsites - type(ed_site_type) , intent(inout), target :: sites(nsites) - type(bc_in_type), intent(in) :: bc_in(nsites) - ! - ! !LOCAL VARIABLES: - integer :: s - real(r8) :: cwd_ag_local(ncwd) - real(r8) :: cwd_bg_local(ncwd) - real(r8) :: spread_local(nclmax) - real(r8) :: leaf_litter_local(numpft_ed) - real(r8) :: root_litter_local(numpft_ed) - real(r8) :: age !notional age of this patch - type(ed_patch_type), pointer :: newp - !---------------------------------------------------------------------- - - cwd_ag_local(:) = 0.0_r8 !ED_val_init_litter -- arbitrary value for litter pools. kgC m-2 - cwd_bg_local(:) = 0.0_r8 !ED_val_init_litter - leaf_litter_local(:) = 0.0_r8 - root_litter_local(:) = 0.0_r8 - spread_local(:) = ED_val_maxspread - age = 0.0_r8 - - !FIX(SPM,032414) clean this up...inits out of this loop - do s = 1, nsites - - allocate(newp) - - newp%patchno = 1 - newp%younger => null() - newp%older => null() - - sites(s)%youngest_patch => newp - sites(s)%youngest_patch => newp - sites(s)%oldest_patch => newp - - ! make new patch... - call create_patch(sites(s), newp, age, AREA, & - spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & - root_litter_local) - - call init_cohorts(newp, bc_in(s)) - - ! This sets the rhizosphere shells based on the plant initialization - ! The initialization of the plant-relevant hydraulics variables - ! were set from a call inside of the init_cohorts()->create_cohort() subroutine - if (use_fates_plant_hydro) then - call updateSizeDepRhizHydProps(sites(s), bc_in(s)) - end if - - enddo + ! + ! !DESCRIPTION: + ! initialize patches + ! This may be call a near bare ground initialization, or it may + ! load patches from an inventory. + + ! + + + use EDParamsMod , only : ED_val_maxspread + use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps + use FatesInventoryInitMod, only : initialize_sites_by_inventory + + ! + ! !ARGUMENTS + integer, intent(in) :: nsites + type(ed_site_type) , intent(inout), target :: sites(nsites) + type(bc_in_type), intent(in) :: bc_in(nsites) + ! + ! !LOCAL VARIABLES: + integer :: s + real(r8) :: cwd_ag_local(ncwd) + real(r8) :: cwd_bg_local(ncwd) + real(r8) :: spread_local(nclmax) + real(r8) :: leaf_litter_local(numpft_ed) + real(r8) :: root_litter_local(numpft_ed) + real(r8) :: age !notional age of this patch + type(ed_patch_type), pointer :: newp + + ! List out some nominal patch values that are used for Near Bear Ground initializations + ! as well as initializing inventory + ! --------------------------------------------------------------------------------------------- + cwd_ag_local(:) = 0.0_r8 !ED_val_init_litter -- arbitrary value for litter pools. kgC m-2 + cwd_bg_local(:) = 0.0_r8 !ED_val_init_litter + leaf_litter_local(:) = 0.0_r8 + root_litter_local(:) = 0.0_r8 + spread_local(:) = ED_val_maxspread + age = 0.0_r8 + ! --------------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------------- + ! Two primary options, either a Near Bear Ground (NBG) or Inventory based cold-start + ! --------------------------------------------------------------------------------------------- + + if (do_inv_init .eq. itrue) then + + call initialize_sites_by_inventory(nsites,sites,bc_in) + + do s = 1, nsites + if (use_fates_plant_hydro) then + call updateSizeDepRhizHydProps(sites(s), bc_in(s)) + end if + enddo + + else + + !FIX(SPM,032414) clean this up...inits out of this loop + do s = 1, nsites + + allocate(newp) + + newp%patchno = 1 + newp%younger => null() + newp%older => null() + + sites(s)%youngest_patch => newp + sites(s)%youngest_patch => newp + sites(s)%oldest_patch => newp + + ! make new patch... + call create_patch(sites(s), newp, age, AREA, & + spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & + root_litter_local) + + call init_cohorts(newp, bc_in(s)) + + ! This sets the rhizosphere shells based on the plant initialization + ! The initialization of the plant-relevant hydraulics variables + ! were set from a call inside of the init_cohorts()->create_cohort() subroutine + if (use_fates_plant_hydro) then + call updateSizeDepRhizHydProps(sites(s), bc_in(s)) + end if + + enddo + + end if end subroutine init_patches diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 old mode 100644 new mode 100755 index c9bc00f3..4ab507ef --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -147,12 +147,15 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) ! puts cohorts in right order call sort_cohorts(currentPatch) - + + ! kills cohorts that are too few + call terminate_cohorts(currentSite, currentPatch, 1) + ! fuses similar cohorts call fuse_cohorts(currentPatch, bc_in ) - ! kills cohorts that are too small - call terminate_cohorts(currentSite, currentPatch) + ! kills cohorts for various other reasons + call terminate_cohorts(currentSite, currentPatch, 2) currentPatch => currentPatch%younger @@ -403,7 +406,9 @@ subroutine ed_update_site( currentSite, bc_in ) currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) - call terminate_cohorts(currentSite, currentPatch) + ! Is termination really needed here? canopy_structure just called it several times! (rgk) + call terminate_cohorts(currentSite, currentPatch, 1) + call terminate_cohorts(currentSite, currentPatch, 2) ! FIX(SPM,040314) why is this needed for BFB restarts? Look into this at some point cohort_number = count_cohorts(currentPatch) @@ -445,8 +450,6 @@ subroutine ed_total_balance_check (currentSite, call_index ) ! Fluxes in are NPP. Fluxes out are decay of CWD and litter into SOM pools. ! ed_allsites_inst%flux_out and ed_allsites_inst%flux_in are set where they occur ! in the code. - use EDTypesMod , only : AREA - ! ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: currentSite @@ -475,7 +478,7 @@ subroutine ed_total_balance_check (currentSite, call_index ) biomass_stock = 0.0_r8 litter_stock = 0.0_r8 - seed_stock = sum(currentSite%seed_bank)*AREA + seed_stock = sum(currentSite%seed_bank) currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 8bc4e9f9..eb994766 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -72,6 +72,8 @@ module EDParamsMod public :: FatesParamsInit public :: FatesRegisterParams public :: FatesReceiveParams + + real(r8), protected :: fates_mortality_disturbance_fraction = 1.0_r8 ! the fraction of canopy mortality that results in disturbance (i.e. transfer of area from new to old patch) contains diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 old mode 100644 new mode 100755 index 785e6cbd..349ffc3a --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -11,7 +11,6 @@ module EDTypesMod save integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site - integer, parameter :: maxCohortsPerPatch = 160 ! maximum number of cohorts to live on a patch integer, parameter :: nclmax = 2 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for understory in two-canopy system @@ -24,6 +23,8 @@ module EDTypesMod integer, parameter :: numpft_ed = 2 ! number of PFTs used in ED. + integer, parameter :: maxCohortsPerPatch = nclmax * numpft_ed * nlevleaf ! maximum number of cohorts to live on a patch + ! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of ! land-ice abledo for vis and nir. This should be a parameter, which would ! get us on track to start using multi-spectral or hyper-spectral (RGK 02-2017) @@ -76,14 +77,14 @@ module EDTypesMod ! SPITFIRE - integer, parameter :: NCWD = 4 ! number of coarse woody debris pools - integer , parameter :: NFSC = NCWD+2 ! number fuel size classes (really this is a mix of cwd size classes, leaf litter, and grass types) + integer, parameter :: NCWD = 4 ! number of coarse woody debris pools (twig,s branch,l branch, trunk) + integer , parameter :: NFSC = NCWD+2 ! number fuel size classes (4 cwd size classes, leaf litter, and grass) integer, parameter :: lg_sf = 6 ! array index of live grass pool for spitfire integer, parameter :: dl_sf = 1 ! array index of dead leaf pool for spitfire (dead grass and dead leaves) integer, parameter :: tw_sf = 2 ! array index of twig pool for spitfire integer, parameter :: tr_sf = 5 ! array index of dead trunk pool for spitfire integer, parameter :: lb_sf = 4 ! array index of large branch pool for spitfire - real(r8), parameter :: fire_threshold = 35.0_r8 ! threshold for fires that spread or go out. KWm-2 + 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.05_r8 ! min plant density for hgt bin to be used in height profile comparisons @@ -97,8 +98,8 @@ module EDTypesMod 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) real(r8), parameter :: min_n_safemath = 1.0E-15_r8 ! in some cases, we want to immediately remove super small - ! number densities of cohorts to prevent FPEs, this is usually - ! just relevant in the first day after recruitment + ! number densities of cohorts to prevent FPEs + character*4 yearchar ! special mode to cause PFTs to create seed mass of all currently-existing PFTs @@ -540,9 +541,10 @@ module EDTypesMod real(r8) :: dseed_dt(numpft_ed) real(r8) :: seed_rain_flux(numpft_ed) ! flux of seeds from exterior KgC/m2/year (needed for C balance purposes) - ! FIRE + ! FIRE + real(r8) :: wind ! daily wind in m/min for Spitfire units real(r8) :: acc_ni ! daily nesterov index accumulating over time. - real(r8) :: ab ! daily burnt area: m2 + real(r8) :: fdi ! daily probability an ignition event will start a fire real(r8) :: frac_burnt ! fraction of soil burnt in this day. real(r8) :: total_burn_flux_to_atm ! total carbon burnt to the atmosphere in this day. KgC/site real(r8) :: cwd_ag_burned(ncwd) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 4d0a6a04..740955d7 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -82,6 +82,12 @@ module FatesHistoryInterfaceMod ! Indices to site by size-class by pft variables integer, private :: ih_nplant_si_scag + integer, private :: ih_nplant_canopy_si_scag + integer, private :: ih_nplant_understory_si_scag + integer, private :: ih_ddbh_canopy_si_scag + integer, private :: ih_ddbh_understory_si_scag + integer, private :: ih_mortality_canopy_si_scag + integer, private :: ih_mortality_understory_si_scag ! Indices to (site) variables integer, private :: ih_nep_si @@ -165,6 +171,12 @@ module FatesHistoryInterfaceMod integer, private :: ih_mortality_understory_si_scls integer, private :: ih_demotion_rate_si_scls integer, private :: ih_promotion_rate_si_scls + integer, private :: ih_trimming_canopy_si_scls + integer, private :: ih_trimming_understory_si_scls + integer, private :: ih_crown_area_canopy_si_scls + integer, private :: ih_crown_area_understory_si_scls + integer, private :: ih_ddbh_canopy_si_scls + integer, private :: ih_ddbh_understory_si_scls ! lots of non-default diagnostics for understanding canopy versus understory carbon balances integer, private :: ih_rdark_canopy_si_scls @@ -279,6 +291,10 @@ module FatesHistoryInterfaceMod integer, private :: ih_fabd_sha_si_cnlf integer, private :: ih_fabi_sun_si_cnlf integer, private :: ih_fabi_sha_si_cnlf + integer, private :: ih_ts_net_uptake_si_cnlf + integer, private :: ih_year_net_uptake_si_cnlf + integer, private :: ih_crownarea_si_cnlf + ! indices to (site x [canopy layer x leaf layer x pft]) variables integer, private :: ih_parsun_z_si_cnlfpft @@ -299,6 +315,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_fabd_sha_top_si_can integer, private :: ih_fabi_sun_top_si_can integer, private :: ih_fabi_sha_top_si_can + integer, private :: ih_crownarea_si_can ! The number of variable dim/kind types we have defined (static) integer, parameter :: fates_history_num_dimensions = 13 @@ -1086,6 +1103,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) use EDParamsMod , only : ED_val_ag_biomass use EDTypesMod , only : get_sizeage_class_index + use EDTypesMod , only : nlevleaf ! Arguments class(fates_history_interface_type) :: this @@ -1106,6 +1124,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) integer :: i_scpf,i_pft,i_scls ! iterators for scpf, pft, and scls dims integer :: i_cwd,i_fuel ! iterators for cwd and fuel dims integer :: iscag ! size-class x age index + integer :: ican, ileaf, cnlf_indx ! iterators for leaf and canopy level real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column @@ -1175,6 +1194,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nplant_understory_si_scpf => this%hvars(ih_nplant_understory_si_scpf)%r82d, & hio_ddbh_canopy_si_scpf => this%hvars(ih_ddbh_canopy_si_scpf)%r82d, & hio_ddbh_understory_si_scpf => this%hvars(ih_ddbh_understory_si_scpf)%r82d, & + hio_ddbh_canopy_si_scls => this%hvars(ih_ddbh_canopy_si_scls)%r82d, & + hio_ddbh_understory_si_scls => this%hvars(ih_ddbh_understory_si_scls)%r82d, & hio_gpp_canopy_si_scpf => this%hvars(ih_gpp_canopy_si_scpf)%r82d, & hio_gpp_understory_si_scpf => this%hvars(ih_gpp_understory_si_scpf)%r82d, & hio_ar_canopy_si_scpf => this%hvars(ih_ar_canopy_si_scpf)%r82d, & @@ -1196,6 +1217,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_demotion_rate_si_scls => this%hvars(ih_demotion_rate_si_scls)%r82d, & hio_demotion_carbonflux_si => this%hvars(ih_demotion_carbonflux_si)%r81d, & hio_promotion_rate_si_scls => this%hvars(ih_promotion_rate_si_scls)%r82d, & + hio_trimming_canopy_si_scls => this%hvars(ih_trimming_canopy_si_scls)%r82d, & + hio_trimming_understory_si_scls => this%hvars(ih_trimming_understory_si_scls)%r82d, & + hio_crown_area_canopy_si_scls => this%hvars(ih_crown_area_canopy_si_scls)%r82d, & + hio_crown_area_understory_si_scls => this%hvars(ih_crown_area_understory_si_scls)%r82d, & hio_promotion_carbonflux_si => this%hvars(ih_promotion_carbonflux_si)%r81d, & hio_canopy_mortality_carbonflux_si => this%hvars(ih_canopy_mortality_carbonflux_si)%r81d, & hio_understory_mortality_carbonflux_si => this%hvars(ih_understory_mortality_carbonflux_si)%r81d, & @@ -1241,7 +1266,15 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_cwd_bg_in_si_cwdsc => this%hvars(ih_cwd_bg_in_si_cwdsc)%r82d, & hio_cwd_ag_out_si_cwdsc => this%hvars(ih_cwd_ag_out_si_cwdsc)%r82d, & hio_cwd_bg_out_si_cwdsc => this%hvars(ih_cwd_bg_out_si_cwdsc)%r82d, & - hio_nplant_si_scag => this%hvars(ih_nplant_si_scag)%r82d) + hio_crownarea_si_cnlf => this%hvars(ih_crownarea_si_cnlf)%r82d, & + hio_crownarea_si_can => this%hvars(ih_crownarea_si_can)%r82d, & + hio_nplant_si_scag => this%hvars(ih_nplant_si_scag)%r82d, & + hio_nplant_canopy_si_scag => this%hvars(ih_nplant_canopy_si_scag)%r82d, & + hio_nplant_understory_si_scag => this%hvars(ih_nplant_understory_si_scag)%r82d, & + hio_ddbh_canopy_si_scag => this%hvars(ih_ddbh_canopy_si_scag)%r82d, & + hio_ddbh_understory_si_scag => this%hvars(ih_ddbh_understory_si_scag)%r82d, & + hio_mortality_canopy_si_scag => this%hvars(ih_mortality_canopy_si_scag)%r82d, & + hio_mortality_understory_si_scag => this%hvars(ih_mortality_understory_si_scag)%r82d) ! --------------------------------------------------------------------------------- @@ -1433,15 +1466,24 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! update SCPF/SCLS- and canopy/subcanopy- partitioned quantities if (ccohort%canopy_layer .eq. 1) then + hio_nplant_canopy_si_scag(io_si,iscag) = hio_nplant_canopy_si_scag(io_si,iscag) + ccohort%n + hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + hio_ddbh_canopy_si_scag(io_si,iscag) = hio_ddbh_canopy_si_scag(io_si,iscag) + & + ccohort%ddbhdt*ccohort%n hio_bstor_canopy_si_scpf(io_si,scpf) = hio_bstor_canopy_si_scpf(io_si,scpf) + & ccohort%bstore * ccohort%n hio_bleaf_canopy_si_scpf(io_si,scpf) = hio_bleaf_canopy_si_scpf(io_si,scpf) + & ccohort%bl * ccohort%n hio_canopy_biomass_pa(io_pa) = hio_canopy_biomass_pa(io_pa) + n_density * ccohort%b * g_per_kg hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n hio_nplant_canopy_si_scpf(io_si,scpf) = hio_nplant_canopy_si_scpf(io_si,scpf) + ccohort%n hio_nplant_canopy_si_scls(io_si,scls) = hio_nplant_canopy_si_scls(io_si,scls) + ccohort%n + hio_trimming_canopy_si_scls(io_si,scls) = hio_trimming_canopy_si_scls(io_si,scls) + & + ccohort%n * ccohort%canopy_trim + hio_crown_area_canopy_si_scls(io_si,scls) = hio_crown_area_canopy_si_scls(io_si,scls) + & + ccohort%c_area hio_gpp_canopy_si_scpf(io_si,scpf) = hio_gpp_canopy_si_scpf(io_si,scpf) + & n_perm2*ccohort%gpp_acc_hold hio_ar_canopy_si_scpf(io_si,scpf) = hio_ar_canopy_si_scpf(io_si,scpf) + & @@ -1449,11 +1491,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! growth increment hio_ddbh_canopy_si_scpf(io_si,scpf) = hio_ddbh_canopy_si_scpf(io_si,scpf) + & ccohort%ddbhdt*ccohort%n + hio_ddbh_canopy_si_scls(io_si,scls) = hio_ddbh_canopy_si_scls(io_si,scls) + & + ccohort%ddbhdt*ccohort%n ! sum of all mortality hio_mortality_canopy_si_scls(io_si,scls) = hio_mortality_canopy_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * & ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 ! hio_leaf_md_canopy_si_scls(io_si,scls) = hio_leaf_md_canopy_si_scls(io_si,scls) + & @@ -1488,15 +1532,24 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_yesterdaycanopylevel_canopy_si_scls(io_si,scls) + & ccohort%canopy_layer_yesterday * ccohort%n else + hio_nplant_understory_si_scag(io_si,iscag) = hio_nplant_understory_si_scag(io_si,iscag) + ccohort%n + hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n + hio_ddbh_understory_si_scag(io_si,iscag) = hio_ddbh_understory_si_scag(io_si,iscag) + & + ccohort%ddbhdt*ccohort%n hio_bstor_understory_si_scpf(io_si,scpf) = hio_bstor_understory_si_scpf(io_si,scpf) + & ccohort%bstore * ccohort%n hio_bleaf_understory_si_scpf(io_si,scpf) = hio_bleaf_understory_si_scpf(io_si,scpf) + & ccohort%bl * ccohort%n hio_understory_biomass_pa(io_pa) = hio_understory_biomass_pa(io_pa) + n_density * ccohort%b * g_per_kg hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n hio_nplant_understory_si_scpf(io_si,scpf) = hio_nplant_understory_si_scpf(io_si,scpf) + ccohort%n hio_nplant_understory_si_scls(io_si,scls) = hio_nplant_understory_si_scls(io_si,scls) + ccohort%n + hio_trimming_understory_si_scls(io_si,scls) = hio_trimming_understory_si_scls(io_si,scls) + & + ccohort%n * ccohort%canopy_trim + hio_crown_area_understory_si_scls(io_si,scls) = hio_crown_area_understory_si_scls(io_si,scls) + & + ccohort%c_area hio_gpp_understory_si_scpf(io_si,scpf) = hio_gpp_understory_si_scpf(io_si,scpf) + & n_perm2*ccohort%gpp_acc_hold hio_ar_understory_si_scpf(io_si,scpf) = hio_ar_understory_si_scpf(io_si,scpf) + & @@ -1504,11 +1557,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! growth increment hio_ddbh_understory_si_scpf(io_si,scpf) = hio_ddbh_understory_si_scpf(io_si,scpf) + & ccohort%ddbhdt*ccohort%n + hio_ddbh_understory_si_scls(io_si,scls) = hio_ddbh_understory_si_scls(io_si,scls) + & + ccohort%ddbhdt*ccohort%n ! sum of all mortality hio_mortality_understory_si_scls(io_si,scls) = hio_mortality_understory_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * ccohort%n hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%imort + ccohort%fmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%fmort) * & ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 ! hio_leaf_md_understory_si_scls(io_si,scls) = hio_leaf_md_understory_si_scls(io_si,scls) + & @@ -1544,10 +1599,32 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%canopy_layer_yesterday * ccohort%n endif ! + ! consider imort as understory mortality even if it happens in cohorts that may have been promoted as part of the patch creation... + hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & + (ccohort%imort) * ccohort%n + hio_mortality_understory_si_scls(io_si,scls) = hio_mortality_understory_si_scls(io_si,scls) + & + (ccohort%imort) * ccohort%n + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & + (ccohort%imort) * & + ccohort%b * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & + (ccohort%imort) * ccohort%n + ! ccohort%canopy_layer_yesterday = real(ccohort%canopy_layer, r8) end associate end if + + ! resolve some canopy area profiles, both total and of occupied leaves + ican = ccohort%canopy_layer + ! + hio_crownarea_si_can(io_si, ican) = hio_crownarea_si_can(io_si, ican) + ccohort%c_area / AREA + ! + do ileaf=1,ccohort%nv + cnlf_indx = ileaf + (ican-1) * nlevleaf + hio_crownarea_si_cnlf(io_si, cnlf_indx) = hio_crownarea_si_cnlf(io_si, cnlf_indx) + & + ccohort%c_area / AREA + end do ccohort => ccohort%taller enddo ! cohort loop @@ -1791,6 +1868,8 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) hio_npp_si_age => this%hvars(ih_npp_si_age)%r82d, & hio_parsun_z_si_cnlf => this%hvars(ih_parsun_z_si_cnlf)%r82d, & hio_parsha_z_si_cnlf => this%hvars(ih_parsha_z_si_cnlf)%r82d, & + hio_ts_net_uptake_si_cnlf => this%hvars(ih_ts_net_uptake_si_cnlf)%r82d, & + hio_year_net_uptake_si_cnlf => this%hvars(ih_year_net_uptake_si_cnlf)%r82d, & hio_parsun_z_si_cnlfpft => this%hvars(ih_parsun_z_si_cnlfpft)%r82d, & hio_parsha_z_si_cnlfpft => this%hvars(ih_parsha_z_si_cnlfpft)%r82d, & hio_laisun_z_si_cnlf => this%hvars(ih_laisun_z_si_cnlf)%r82d, & @@ -1909,45 +1988,61 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! accumulate fluxes on canopy- and understory- separated fluxes if (ccohort%canopy_layer .eq. 1) then + ! + ! bulk fluxes are in gC / m2 / s hio_gpp_canopy_pa(io_pa) = hio_gpp_canopy_pa(io_pa) + & ccohort%gpp_tstep * g_per_kg * n_density * per_dt_tstep hio_ar_canopy_pa(io_pa) = hio_ar_canopy_pa(io_pa) + & ccohort%resp_tstep * g_per_kg * n_density * per_dt_tstep ! + ! size-resolved respiration fluxes are in kg C / ha / yr hio_rdark_canopy_si_scls(io_si,scls) = hio_rdark_canopy_si_scls(io_si,scls) + & - ccohort%rdark * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%rdark * ccohort%n * sec_per_day * days_per_year hio_livestem_mr_canopy_si_scls(io_si,scls) = hio_livestem_mr_canopy_si_scls(io_si,scls) + & - ccohort%livestem_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%livestem_mr * ccohort%n * sec_per_day * days_per_year hio_livecroot_mr_canopy_si_scls(io_si,scls) = hio_livecroot_mr_canopy_si_scls(io_si,scls) + & - ccohort%livecroot_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%livecroot_mr * ccohort%n * sec_per_day * days_per_year hio_froot_mr_canopy_si_scls(io_si,scls) = hio_froot_mr_canopy_si_scls(io_si,scls) + & - ccohort%froot_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%froot_mr * ccohort%n * sec_per_day * days_per_year hio_resp_g_canopy_si_scls(io_si,scls) = hio_resp_g_canopy_si_scls(io_si,scls) + & - ccohort%resp_g * g_per_kg * ccohort%n * sec_per_day * days_per_year * per_dt_tstep + ccohort%resp_g * ccohort%n * sec_per_day * days_per_year * per_dt_tstep hio_resp_m_canopy_si_scls(io_si,scls) = hio_resp_m_canopy_si_scls(io_si,scls) + & - ccohort%resp_m * g_per_kg * ccohort%n * sec_per_day * days_per_year * per_dt_tstep + ccohort%resp_m * ccohort%n * sec_per_day * days_per_year * per_dt_tstep else + ! + ! bulk fluxes are in gC / m2 / s hio_gpp_understory_pa(io_pa) = hio_gpp_understory_pa(io_pa) + & ccohort%gpp_tstep * g_per_kg * n_density * per_dt_tstep hio_ar_understory_pa(io_pa) = hio_ar_understory_pa(io_pa) + & ccohort%resp_tstep * g_per_kg * n_density * per_dt_tstep ! + ! size-resolved respiration fluxes are in kg C / ha / yr hio_rdark_understory_si_scls(io_si,scls) = hio_rdark_understory_si_scls(io_si,scls) + & - ccohort%rdark * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%rdark * ccohort%n * sec_per_day * days_per_year hio_livestem_mr_understory_si_scls(io_si,scls) = hio_livestem_mr_understory_si_scls(io_si,scls) + & - ccohort%livestem_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%livestem_mr * ccohort%n * sec_per_day * days_per_year hio_livecroot_mr_understory_si_scls(io_si,scls) = hio_livecroot_mr_understory_si_scls(io_si,scls) + & - ccohort%livecroot_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%livecroot_mr * ccohort%n * sec_per_day * days_per_year hio_froot_mr_understory_si_scls(io_si,scls) = hio_froot_mr_understory_si_scls(io_si,scls) + & - ccohort%froot_mr * g_per_kg * ccohort%n * sec_per_day * days_per_year + ccohort%froot_mr * ccohort%n * sec_per_day * days_per_year hio_resp_g_understory_si_scls(io_si,scls) = hio_resp_g_understory_si_scls(io_si,scls) + & - ccohort%resp_g * g_per_kg * ccohort%n * sec_per_day * days_per_year * per_dt_tstep + ccohort%resp_g * ccohort%n * sec_per_day * days_per_year * per_dt_tstep hio_resp_m_understory_si_scls(io_si,scls) = hio_resp_m_understory_si_scls(io_si,scls) + & - ccohort%resp_m * g_per_kg * ccohort%n * sec_per_day * days_per_year * per_dt_tstep + ccohort%resp_m * ccohort%n * sec_per_day * days_per_year * per_dt_tstep endif end associate endif + !!! resolve some canopy profile terms that are also on the cohort indices + ican = ccohort%canopy_layer + do ileaf=1,ccohort%nv + cnlf_indx = ileaf + (ican-1) * nlevleaf + hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & + ccohort%ts_net_uptake(ileaf) * ccohort%c_area / AREA + hio_year_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_year_net_uptake_si_cnlf(io_si, cnlf_indx) + & + ccohort%year_net_uptake(ileaf) * ccohort%c_area / AREA + end do + ccohort => ccohort%taller enddo ! cohort loop @@ -2651,7 +2746,7 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_ar_canopy_pa ) call this%set_history_var(vname='GPP_UNDERSTORY', units='gC/m^2/s', & - long='gross primary production of understory plants', use_default='inactive', & + long='gross primary production of understory plants', use_default='active', & avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & ivar=ivar, initialize=initialize_variables, index = ih_gpp_understory_pa ) @@ -2806,6 +2901,31 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_can_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & ivar=ivar, initialize=initialize_variables, index = ih_fabi_sha_top_si_can ) + !!! canopy-resolved fluxes and structure + call this%set_history_var(vname='TS_NET_UPTAKE_CNLF', units='kgC/m2/s', & + long='net carbon uptake by each canopy and leaf layer er unit ground area (i.e. divide by CROWNAREA_CNLF)', & + use_default='inactive', & + avgflag='A', vtype=site_cnlf_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_ts_net_uptake_si_cnlf ) + + call this%set_history_var(vname='YEAR_NET_UPTAKE_CNLF', units='kgC/m2/y', & + long='yearly net carbon uptake by each canopy and leaf layer per unit ground area (i.e. divide by CROWNAREA_CNLF)', & + use_default='inactive', & + avgflag='A', vtype=site_cnlf_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_year_net_uptake_si_cnlf ) + + call this%set_history_var(vname='CROWNAREA_CNLF', units='m2/m2', & + long='total crown area that is occupied by leaves in each canopy and leaf layer', & + use_default='inactive', & + avgflag='A', vtype=site_cnlf_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_crownarea_si_cnlf ) + + call this%set_history_var(vname='CROWNAREA_CAN', units='m2/m2', & + long='total crown area in each canopy layer', & + use_default='active', & + avgflag='A', vtype=site_can_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_crownarea_si_can ) + ! slow carbon fluxes associated with mortality from or transfer betweeen canopy and understory call this%set_history_var(vname='DEMOTION_CARBONFLUX', units = 'gC/m2/s', & long='demotion-associated biomass carbon flux from canopy to understory', use_default='active', & @@ -2818,12 +2938,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_promotion_carbonflux_si ) call this%set_history_var(vname='MORTALITY_CARBONFLUX_CANOPY', units = 'gC/m2/s', & - long='flux of biomass carbon from live to dead pools from mortality of canopy plants', use_default='active', & + long='flux of biomass carbon from live to dead pools from mortality of canopy plants', use_default='inactive', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_canopy_mortality_carbonflux_si ) call this%set_history_var(vname='MORTALITY_CARBONFLUX_UNDERSTORY', units = 'gC/m2/s', & - long='flux of biomass carbon from live to dead pools from mortality of understory plants', use_default='active', & + long='flux of biomass carbon from live to dead pools from mortality of understory plants',use_default='inactive',& avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_understory_mortality_carbonflux_si ) @@ -2833,6 +2953,36 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_scag ) + call this%set_history_var(vname='NPLANT_CANOPY_SCAG',units = 'plants/ha', & + long='number of plants per hectare in canopy in each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_canopy_si_scag ) + + call this%set_history_var(vname='NPLANT_UNDERSTORY_SCAG',units = 'plants/ha', & + long='number of plants per hectare in understory in each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_understory_si_scag ) + + call this%set_history_var(vname='DDBH_CANOPY_SCAG',units = 'cm/yr/ha', & + long='growth rate of canopy plantsnumber of plants per hectare in canopy in each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_canopy_si_scag ) + + call this%set_history_var(vname='DDBH_UNDERSTORY_SCAG',units = 'cm/yr/ha', & + long='growth rate of understory plants in each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_understory_si_scag ) + + call this%set_history_var(vname='MORTALITY_CANOPY_SCAG',units = 'plants/ha/yr', & + long='mortality rate of canopy plants in each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_mortality_canopy_si_scag ) + + call this%set_history_var(vname='MORTALITY_UNDERSTORY_SCAG',units = 'plants/ha/yr', & + long='mortality rate of understory plantsin each size x age class', use_default='inactive', & + avgflag='A', vtype=site_scag_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_mortality_understory_si_scag ) + ! Carbon Flux (grid dimension x scpf) (THESE ARE DEFAULT INACTIVE!!! ! (BECAUSE THEY TAKE UP SPACE!!! @@ -3073,6 +3223,16 @@ subroutine define_history_vars(this, initialize_variables) ! size-class only variables + call this%set_history_var(vname='DDBH_CANOPY_SCLS', units = 'cm/yr/ha', & + long='diameter growth increment by pft/size',use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_canopy_si_scls ) + + call this%set_history_var(vname='DDBH_UNDERSTORY_SCLS', units = 'cm/yr/ha', & + long='diameter growth increment by pft/size',use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_understory_si_scls ) + call this%set_history_var(vname='YESTERDAYCANLEV_CANOPY_SCLS', units = 'indiv/ha', & long='Yesterdays canopy level for canopy plants by size class', use_default='inactive', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -3118,6 +3278,26 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_mortality_understory_si_scls ) + call this%set_history_var(vname='TRIMMING_CANOPY_SCLS', units = 'indiv/ha', & + long='trimming term of canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_trimming_canopy_si_scls ) + + call this%set_history_var(vname='TRIMMING_UNDERSTORY_SCLS', units = 'indiv/ha', & + long='trimming term of understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_trimming_understory_si_scls ) + + call this%set_history_var(vname='CROWN_AREA_CANOPY_SCLS', units = 'm2/ha', & + long='total crown area of canopy plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_crown_area_canopy_si_scls ) + + call this%set_history_var(vname='CROWN_AREA_UNDERSTORY_SCLS', units = 'm2/ha', & + long='total crown area of understory plants by size class', use_default='inactive', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_crown_area_understory_si_scls ) + call this%set_history_var(vname='LEAF_MD_CANOPY_SCLS', units = 'kg C / ha / yr', & long='LEAF_MD for canopy plants by size class', use_default='inactive', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 new file mode 100644 index 00000000..28ab8dd2 --- /dev/null +++ b/main/FatesInventoryInitMod.F90 @@ -0,0 +1,916 @@ +module FatesInventoryInitMod + + !----------------------------------------------------------------------------------------------- + ! This module contains the majority of the code used to read forest inventory data and + ! initialize a simulation from that data. Some important points: + ! - This procedure is called through the host model's "cold-start" initialization and not a + ! restart type of simulation. + ! - This procedure, if called from a massive grid, is both probably inappropriate, and probably + ! time consuming. + ! - This procedure is assumed to be called over a small subset of sites, for instance a single + ! site, or a small collection of sparse/irregularly spaced group of sites + ! + ! Created: Ryan Knox June 2017 + ! This code borrows heavily in concept from what is done in ED2. We will also do our best to + ! maintain compatibility with the PSS/CSS file formats that were used in ED2. + ! See: https://github.com/EDmodel/ED2/blob/master/ED/src/io/ed_read_ed10_20_history.f90 + ! At the time of writing this ED2 is unlicensed, and only concepts were borrowed with no direct + ! code copied. + !----------------------------------------------------------------------------------------------- + + ! CIME GLOBALS + + use shr_log_mod , only : errMsg => shr_log_errMsg + + ! FATES GLOBALS + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : pi_const + use FatesGlobals , only : endrun => fates_endrun + use FatesGlobals , only : fates_log + use FatesInterfaceMod, only : bc_in_type + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : ed_cohort_type + use EDTypesMod , only : area + use EDPftvarcon , only : EDPftvarcon_inst + use EDEcophysConType , only : EDecophyscon + + implicit none + private + + ! This derived type is to allow an array of pointers to the LL patch structure + ! This is different than allocating a vector of patches. This is needed for + ! quickly matching cohort string identifiers, the indices that match thos identifiers + ! with a patch. BY having a vector of patch pointers that lines up with the string + ! identifier array, this can be done quickly. + type pp_array + type(ed_patch_type), pointer :: cpatch + end type pp_array + + ! For now we will use a hard-coded file name for the inventory file list + character(len=*), parameter :: inv_file_list = 'inventory_file_list.txt' + + character(len=*), parameter, private :: sourcefile = __FILE__ + + logical, parameter :: debug_inv = .false. ! Debug flag for devs + + ! String length specifiers + integer, parameter :: patchname_strlen = 64 + integer, parameter :: cohortname_strlen = 64 + integer, parameter :: line_strlen = 512 + integer, parameter :: path_strlen = 256 + + + real(r8), parameter :: max_site_adjacency_deg = 0.05_r8 ! The maximum distance in degrees + ! allowed between a site's coordinate + ! defined in model memory and a physical + ! site listed in the file + + public :: initialize_sites_by_inventory + +contains + + ! ============================================================================================== + + subroutine initialize_sites_by_inventory(nsites,sites,bc_in) + + ! !USES: + use shr_file_mod, only : shr_file_getUnit + use shr_file_mod, only : shr_file_freeUnit + use EDTypesMod, only : nclmax + use EDTypesMod, only : numpft_ed + use EDTypesMod, only : maxpft + use EDTypesMod, only : ncwd + use EDParamsMod, only : ED_val_maxspread + use EDPatchDynamicsMod, only : create_patch + use EDPatchDynamicsMod, only : fuse_patches + use EDCohortDynamicsMod, only : fuse_cohorts + use EDCohortDynamicsMod, only : sort_cohorts + use EDcohortDynamicsMod, only : count_cohorts + + ! Arguments + integer, intent(in) :: nsites + type(ed_site_type), intent(inout), target :: sites(nsites) + type(bc_in_type), intent(in) :: bc_in(nsites) + + ! Locals + type(ed_patch_type), pointer :: currentpatch + type(ed_cohort_type), pointer :: currentcohort + type(ed_patch_type), pointer :: newpatch + type(ed_patch_type), pointer :: olderpatch + integer :: sitelist_file_unit ! fortran file unit for site list + integer :: pss_file_unit ! fortran file unit for the pss file + integer :: css_file_unit ! fortran file unit for the css file + integer :: nfilesites ! number of sites in file list + logical :: lopen ! logical, file is open + logical :: lexist ! logical, file exists + integer :: ios ! integer, "IO" status + character(len=line_strlen) :: header_str ! large string for whole lines + real(r8) :: age_init ! dummy value for creating a patch + real(r8) :: area_init ! dummy value for creating a patch + real(r8) :: spread_init(nclmax) ! dummy value for creating a patch + real(r8) :: cwd_ag_init(ncwd) ! dummy value for creating a patch + real(r8) :: cwd_bg_init(ncwd) ! dummy value for creating a patch + real(r8) :: leaf_litter_init(maxpft) ! dummy value for creating a patch + real(r8) :: root_litter_init(maxpft) ! dummy value for creating a patch + integer :: s ! site index + integer :: ipa ! patch index + integer :: total_cohorts ! cohort counter for error checking + integer, allocatable :: inv_format_list(:) ! list of format specs + character(len=path_strlen), allocatable :: inv_css_list(:) ! list of css file names + character(len=path_strlen), allocatable :: inv_pss_list(:) ! list of pss file names + real(r8), allocatable :: inv_lat_list(:) ! list of lat coords + real(r8), allocatable :: inv_lon_list(:) ! list of lon coords + integer :: invsite ! index of inventory site + ! closest to actual site + character(len=patchname_strlen) :: patch_name ! patch ID string in the file + integer :: npatches ! number of patches found in PSS + type(pp_array), allocatable :: patch_pointer_vec(:) ! vector of pointers to patch LL + character(len=patchname_strlen), allocatable :: patch_name_vec(:) ! vector of patch ID strings + real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) + real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) + + real(r8), parameter :: max_ba_diff = 1.0e-2 ! 1% is the maximum allowable + ! change in BA due to fusion + + ! I. Load the inventory list file, do some file handle checks + ! ------------------------------------------------------------------------------------------ + + sitelist_file_unit = shr_file_getUnit() + inquire(file=trim(inv_file_list),exist=lexist,opened=lopen) + if( .not.lexist ) then ! The inventory file list DNE + write(fates_log(), *) 'An inventory Initialization was requested.' + write(fates_log(), *) 'However the inventory file: ',trim(inv_file_list),' DNE' + write(fates_log(), *) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if( lopen ) then ! The inventory file should not be open + write(fates_log(), *) 'The inventory list file is open but should not be.' + write(fates_log(), *) 'Aborting.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + open(unit=sitelist_file_unit,file=trim(inv_file_list),status='OLD',action='READ',form='FORMATTED') + rewind(sitelist_file_unit) + + ! There should be at least 1 line + read(sitelist_file_unit,fmt='(A)',iostat=ios) header_str + read(sitelist_file_unit,fmt='(A)',iostat=ios) header_str + if( ios /= 0 ) then + write(fates_log(), *) 'The inventory file does not contain at least two lines' + write(fates_log(), *) 'of data, ie a header and 1 site. Aborting.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + rewind(unit=sitelist_file_unit) + + + ! Count the number of sites that are listed in this file, and allocate storage arrays + ! ------------------------------------------------------------------------------------------ + + nfilesites = count_inventory_sites(sitelist_file_unit) + + allocate(inv_format_list(nfilesites)) + allocate(inv_pss_list(nfilesites)) + allocate(inv_css_list(nfilesites)) + allocate(inv_lat_list(nfilesites)) + allocate(inv_lon_list(nfilesites)) + + + ! Check through the sites that are listed and do some sanity checks + ! ------------------------------------------------------------------------------------------ + call assess_inventory_sites(sitelist_file_unit, nfilesites, inv_format_list, & + inv_pss_list, inv_css_list, & + inv_lat_list, inv_lon_list) + + ! We can close the list file now + close(sitelist_file_unit, iostat = ios) + if( ios /= 0 ) then + write(fates_log(), *) 'The inventory file needed to be closed, but was still open' + write(fates_log(), *) 'aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + call shr_file_freeUnit(sitelist_file_unit) + + + ! For each site, identify the most proximal PSS/CSS couplet, read-in the data + ! allocate linked lists and assign to memory + do s = 1, nsites + invsite = & + minloc( (sites(s)%lat-inv_lat_list(:))**2.0_r8 + & + (sites(s)%lon-inv_lon_list(:))**2.0_r8 , dim=1) + + ! Do a sanity check on the distance separation between physical site and model site + if ( sqrt( (sites(s)%lat-inv_lat_list(invsite))**2.0_r8 + & + (sites(s)%lon-inv_lon_list(invsite))**2.0_r8 ) > max_site_adjacency_deg ) then + write(fates_log(), *) 'Model site at lat:',sites(s)%lat,' lon:',sites(s)%lon + write(fates_log(), *) 'has no reasonably proximal site in the inventory site list.' + write(fates_log(), *) 'Closest is at lat:',inv_lat_list(invsite),' lon:',inv_lon_list(invsite) + write(fates_log(), *) 'Separation must be less than ',max_site_adjacency_deg,' degrees' + write(fates_log(), *) 'Exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + ! Open the PSS/CSS couplet and initialize the ED data structures. + ! Lets start withe the PSS + ! --------------------------------------------------------------------------------------- + + pss_file_unit = shr_file_getUnit() + open(unit=pss_file_unit,file=trim(inv_pss_list(invsite)), & + status='OLD',action='READ',form='FORMATTED') + rewind(pss_file_unit) + read(pss_file_unit,fmt=*) header_str + + ! Do one quick pass through just to count lines + ipa = 0 + countpatchloop: do + read(pss_file_unit,fmt=*,iostat=ios) header_str + if(ios/=0) exit + ipa = ipa + 1 + end do countpatchloop + rewind(pss_file_unit) + read(pss_file_unit,fmt=*) header_str + + npatches = ipa + allocate(patch_name_vec(npatches)) + allocate(patch_pointer_vec(npatches)) + + + do ipa=1,npatches + + allocate(newpatch) + + newpatch%patchno = ipa + newpatch%younger => null() + newpatch%older => null() + + ! This call doesn't do much asside from initializing the patch with + ! nominal values, NaNs, zero's and allocating some vectors. We should + ! be able to get the following values from the patch files. But on + ! the patch creation step, we don't have that information. + + age_init = 0.0_r8 + area_init = 0.0_r8 + spread_init(:) = ED_val_maxspread + cwd_ag_init(:) = 0.0_r8 + cwd_bg_init(:) = 0.0_r8 + leaf_litter_init(1:numpft_ed) = 0.0_r8 + root_litter_init(1:numpft_ed) = 0.0_r8 + + call create_patch(sites(s), newpatch, age_init, area_init, spread_init, & + cwd_ag_init, cwd_bg_init, & + leaf_litter_init(1:numpft_ed), root_litter_init(1:numpft_ed) ) + + if( inv_format_list(invsite) == 1 ) then + call set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) + end if + + ! Add it to the site's patch list + ! ------------------------------------------------------------------------------------ + + patch_name_vec(ipa) = trim(patch_name) + patch_pointer_vec(ipa)%cpatch => newpatch + + if(ipa == 1) then + ! This is the first patch to be added + ! It starts off as the oldest and youngest patch in the list + sites(s)%youngest_patch => newpatch + sites(s)%oldest_patch => newpatch + else + + ! Insert this patch into the patch LL + ! First check the two end-cases + + ! Youngest Patch + if(newpatch%age<=sites(s)%youngest_patch%age)then + newpatch%older => sites(s)%youngest_patch + newpatch%younger => null() + sites(s)%youngest_patch%younger => newpatch + sites(s)%youngest_patch => newpatch + + ! Oldest Patch + else if(newpatch%age>sites(s)%oldest_patch%age)then + newpatch%older => null() + newpatch%younger => sites(s)%oldest_patch + sites(s)%oldest_patch%older => newpatch + sites(s)%oldest_patch => newpatch + + ! Somewhere in the middle + else + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + olderpatch => currentpatch%older + if(associated(currentpatch%older)) then + if(newpatch%age >= currentpatch%age .and. & + newpatch%age < olderpatch%age) then + ! Set the new patches pointers + newpatch%older => currentpatch%older + newpatch%younger => currentpatch + ! Fix the patch's older pointer + currentpatch%older => newpatch + ! Fix the older patch's younger pointer + olderpatch%younger => newpatch + end if + end if + currentPatch => olderpatch + enddo + end if + end if + end do + + close(pss_file_unit,iostat=ios) + if( ios /= 0 ) then + write(fates_log(), *) 'The pss file: ',inv_pss_list(invsite),' could not be closed' + write(fates_log(), *) 'aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + call shr_file_freeUnit(pss_file_unit) + + if(debug_inv) then + write(fates_log(),*) 'Raw List of Inventory Patches, Age Sorted:' + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + write(fates_log(),*) ' AGE: ',currentpatch%age,' AREA: ',currentpatch%area + currentPatch => currentpatch%older + enddo + end if + + ! OPEN THE CSS FILE + ! --------------------------------------------------------------------------------------- + css_file_unit = shr_file_getUnit() + open(unit=css_file_unit,file=trim(inv_css_list(invsite)), & + status='OLD',action='READ',form='FORMATTED') + rewind(css_file_unit) + read(css_file_unit,fmt=*) header_str + + ! Read in each cohort line. Each line is associated with a patch from the PSS + ! file via a patch name identification string. We pass the whole site pointer + ! to this routine, because inside the routine we identify the patch by making + ! comparisons with patch_name_vec and identifying the patch pointer + ! from patch_pointer_vec + + invcohortloop: do + if ( inv_format_list(invsite) == 1 ) then + call set_inventory_edcohort_type1(sites(s),bc_in(s),css_file_unit, & + npatches, patch_pointer_vec,patch_name_vec, ios) + end if + if ( ios/=0 ) exit + end do invcohortloop + + close(css_file_unit,iostat=ios) + if( ios/=0 ) then + write(fates_log(), *) 'The css file: ',inv_css_list(invsite),' could not be closed' + write(fates_log(), *) 'aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + call shr_file_freeUnit(css_file_unit) + + deallocate(patch_pointer_vec,patch_name_vec) + + ! Report Basal Area (as a check on if things were read in) + ! ------------------------------------------------------------------------------ + basal_area_pref = 0.0_r8 + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + currentcohort => currentpatch%tallest + do while(associated(currentcohort)) + basal_area_pref = basal_area_pref + & + currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const + currentcohort => currentcohort%shorter + end do + currentPatch => currentpatch%older + enddo + + write(fates_log(),*) '-------------------------------------------------------' + write(fates_log(),*) 'Basal Area from inventory, BEFORE fusion' + write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon + write(fates_log(),*) basal_area_pref,' [m2/ha]' + write(fates_log(),*) '-------------------------------------------------------' + + ! Update the patch index numbers and fuse the cohorts in the patches + ! ---------------------------------------------------------------------------------------- + ipa=1 + total_cohorts = 0 + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + currentpatch%patchno = ipa + ipa=ipa+1 + + ! Perform Cohort Fusion + call fuse_cohorts(currentpatch,bc_in(s)) + call sort_cohorts(currentpatch) + total_cohorts = total_cohorts + count_cohorts(currentpatch) + + currentPatch => currentpatch%older + enddo + + if(total_cohorts .eq. 0)then + write(fates_log(), *) 'The inventory initialization produced no cohorts.' + write(fates_log(), *) 'aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Fuse patches + ! ---------------------------------------------------------------------------------------- + call fuse_patches(sites(s), bc_in(s) ) + + ! Report Basal Area (as a check on if things were read in) + ! ---------------------------------------------------------------------------------------- + basal_area_postf = 0.0_r8 + currentpatch => sites(s)%youngest_patch + do while(associated(currentpatch)) + currentcohort => currentpatch%tallest + do while(associated(currentcohort)) + basal_area_postf = basal_area_postf + & + currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const + currentcohort => currentcohort%shorter + end do + currentPatch => currentpatch%older + enddo + + write(fates_log(),*) '-------------------------------------------------------' + write(fates_log(),*) 'Basal Area from inventory, AFTER fusion' + write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon + write(fates_log(),*) basal_area_postf,' [m2/ha]' + write(fates_log(),*) '-------------------------------------------------------' + + ! Check to see if the fusion process has changed too much + ! We are sensitive to fusion in inventories because we may be asking for a massive amount + ! of fusion. For instance some init files are directly from inventory, where a cohort + ! is synomomous with a single plant. + + if( abs(basal_area_postf-basal_area_pref)/basal_area_pref > max_ba_diff ) then + write(fates_log(),*) 'Inventory Fusion Changed total biomass beyond reasonable limit' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end do + + deallocate(inv_format_list, inv_pss_list, inv_css_list, inv_lat_list, inv_lon_list) + + return + end subroutine initialize_sites_by_inventory + + ! ============================================================================================== + + function count_inventory_sites(sitelist_file_unit) result(nsites) + + ! Simple function that counts the number of lines in the inventory descriptor file + + ! Arguments + integer, intent(in) :: sitelist_file_unit + + ! Locals + character(len=line_strlen) :: header_str + character(len=line_strlen) :: site_str + integer :: ios + integer :: nsites + + + ! Set the file position to the top of the file + ! Read in the header line + ! Read through sites and check coordinates and file existence + rewind(unit=sitelist_file_unit) + read(sitelist_file_unit,fmt='(A)') header_str + nsites = 0 + do + read(sitelist_file_unit,fmt='(A)',iostat=ios) site_str + if (ios/=0) exit + nsites = nsites + 1 + end do + + return + end function count_inventory_sites + + ! ============================================================================================== + + subroutine assess_inventory_sites(sitelist_file_unit,nsites, inv_format_list, & + inv_pss_list,inv_css_list, & + inv_lat_list,inv_lon_list) + + ! ------------------------------------------------------------------------------------------- + ! This subroutine looks through the inventory descriptor file + ! and line by line reads information about the available inventory + ! sites, and saves their information (such as location and file path) + ! to arrays. This routine also does some simple checks to make + ! sure it is not reading nonsense + ! + ! File Format for the inventory site file: + ! 1 line header + ! 1 line listing each available inventory site with the following fields: + ! type latitude longitude pss-name css-name + ! + ! The fields for each site are described as follows: + ! + ! + ! + ! type integer We will accomodate different file format with different + ! field values as the need arises. format 1 will read in + ! datasets via "set_inventory_edpatch_type1()", + ! "set_inventory_edcohort_type1()" + ! + ! latitude float The geographic latitude coordinate of the site + ! longitude float The geogarphic longitude coordinate of the site + ! pss-name string The full path to the patch descriptor file (PSS) + ! css-name string The full path to the cohort descriptor file (CSS) + ! ------------------------------------------------------------------------------------------- + + + ! Arguments + integer, intent(in) :: sitelist_file_unit ! file unit for sitelist + integer, intent(in) :: nsites ! number of inventory sites + integer, intent(inout) :: inv_format_list(nsites) ! array of formats for each inventory site + character(len=path_strlen),intent(inout) :: inv_pss_list(nsites) ! array of pss file paths for each site + character(len=path_strlen),intent(inout) :: inv_css_list(nsites) ! array of css file paths for each site + real(r8),intent(inout) :: inv_lat_list(nsites) ! array of latitudes for each site + real(r8),intent(inout) :: inv_lon_list(nsites) ! array of longitudes for each site + + ! Locals + character(len=line_strlen) :: header_str ! a string to hold the header information + character(len=line_strlen) :: site_str ! a string to hold each site-line in the file + integer :: isite ! site index + integer :: ios ! fortran read status flag + character(len=path_strlen) :: pss_file + character(len=path_strlen) :: css_file + real(r8) :: site_lat ! inventory site latitude + real(r8) :: site_lon ! site longitude + integer :: iblnk ! Index used for string parsing + integer :: file_format ! format type (1=legacy ED pss/css) + logical :: lexist ! file existence flag + + rewind(unit=sitelist_file_unit) + read(sitelist_file_unit,fmt='(4A)') header_str + + do isite=1,nsites + + ! Read in the whole line + read(sitelist_file_unit,fmt='(a)',iostat=ios) site_str + + ! Parse the format identifier + read(site_str,*) file_format + + ! Parse the latitude + site_str = adjustl(site_str) + iblnk = index(site_str,' ') + site_str = adjustl(site_str(iblnk:)) + read(site_str,*) site_lat + + ! Parse the longitude + site_str = adjustl(site_str) + iblnk = index(site_str,' ') + site_str = adjustl(site_str(iblnk:)) + read(site_str,*) site_lon + + ! Parse the pss file name + site_str = adjustl(site_str) + iblnk = index(site_str,' ') + site_str = adjustl(site_str(iblnk:)) + iblnk = index(site_str,' ') + read(site_str(:iblnk),fmt='(1A)') pss_file + + ! Parse the css file name + site_str = adjustl(site_str) + iblnk = index(site_str,' ') + site_str = adjustl(site_str(iblnk:)) + read(site_str,fmt='(1A)') css_file + + + if ( site_lat < -90.0_r8 .or. site_lat > 90.0_r8 ) then + write(fates_log(), *) 'read invalid latitude: ',site_lat,' from inventory site list' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Longitude should be converted to positive coordinate + + if( site_lon<0.0_r8 ) site_lon = 360.0_r8 + site_lon + + if ( site_lon < 0.0_r8 .or. site_lon > 360.0_r8 ) then + write(fates_log(), *) 'read invalid longitude: ',site_lon,' from inventory site list' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + inquire(file=trim(pss_file),exist=lexist) + if( .not.lexist ) then + write(fates_log(), *) 'the following pss file could not be found:' + write(fates_log(), *) trim(pss_file) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + inquire(file=trim(css_file),exist=lexist) + if( .not.lexist ) then + write(fates_log(), *) 'the following css file could not be found:' + write(fates_log(), *) trim(css_file) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! If we have made it to this point, then in all likelihood, the PSS/CSS + ! File has probably been correctly identified + + inv_format_list(isite) = file_format + inv_pss_list(isite) = pss_file + inv_css_list(isite) = css_file + inv_lat_list(isite) = site_lat + inv_lon_list(isite) = site_lon + + end do + + return + end subroutine assess_inventory_sites + + ! ============================================================================================== + + subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) + + ! -------------------------------------------------------------------------------------------- + ! This subroutine reads in a line of an inventory patch file (pss) + ! And populates a new patch with that information. + ! This routine specifically reads PSS files that are "Type 1" formatted + ! + ! The file is formatted text, which contains 1 header line to label columns + ! and then 1 line for each patch containing the following fields: + ! + ! time (year) year of measurement + ! patch (string) patch id string + ! trk (integer) LU type index (0 non-forest, 1 secondary, 2 primary + ! age (years) Time since this patch was disturbed (created) + ! area (fraction) Fraction of the site occupied by this patch + ! water (NA) Water content of soil (NOT USED) + ! fsc (kg/m2) Fast Soil Carbon + ! stsc (kg/m2) Structural Soil Carbon + ! stsl (kg/m2) Structural Soil Lignan + ! ssc (kg/m2) Slow Soil Carbon + ! psc (NA) Passive Soil Carbon (NOT USED) + ! msn (kg/m2) Mineralized Soil Nitrogen + ! fsn (kg/m2) Fast Soil Nitrogen + ! -------------------------------------------------------------------------------------------- + + use EDTypesMod, only: get_age_class_index + use EDtypesMod, only: AREA + use EDTypesMod, only: numpft_ed + use EDTypesMod, only: ncwd + use SFParamsMod , only : SF_val_CWD_frac + use EDParamsMod , only : ED_val_ag_biomass + + ! Arguments + type(ed_patch_type),intent(inout), target :: newpatch ! Patch structure + integer,intent(in) :: pss_file_unit ! Self explanatory + integer,intent(in) :: ipa ! Patch index (line number) + integer,intent(out) :: ios ! Return flag + character(len=patchname_strlen),intent(out) :: patch_name ! unique string identifier of patch + + ! Locals + real(r8) :: p_time ! Time patch was recorded + real(r8) :: p_trk ! Land Use index (see above descriptions) + character(len=patchname_strlen) :: p_name ! unique string identifier of patch + real(r8) :: p_age ! Patch age [years] + real(r8) :: p_area ! Patch area [fraction] + real(r8) :: p_water ! Patch water (unused) + real(r8) :: p_fsc ! Patch fast soil carbon + real(r8) :: p_stsc ! Patch structural soil carbon + real(r8) :: p_stsl ! Patch structural soil lignans + real(r8) :: p_ssc ! Patch slow soil carbon + real(r8) :: p_psc ! Patch P soil carbon + real(r8) :: p_msn ! Patch mean soil nitrogen + real(r8) :: p_fsn ! Patch fast soil nitrogen + integer :: icwd ! index for counting CWD pools + integer :: ipft ! index for counting PFTs + real(r8) :: pftfrac ! the inverse of the total number of PFTs + + character(len=128),parameter :: wr_fmt = & + '(F5.2,2X,A4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' + + + read(pss_file_unit,fmt=*,iostat=ios) p_time, p_name, p_trk, p_age, p_area, & + p_water,p_fsc, p_stsc, p_stsl, p_ssc, & + p_psc, p_msn, p_fsn + + if (ios/=0) return + + patch_name = trim(p_name) + + if( debug_inv) then + write(*,fmt=wr_fmt) & + p_time, p_name, p_trk, p_age, p_area, & + p_water,p_fsc, p_stsc, p_stsl, p_ssc, & + p_psc, p_msn, p_fsn + end if + + ! Fill in the patch's memory structures + + newpatch%age = p_age + newpatch%age_class = get_age_class_index(newpatch%age) + newpatch%area = p_area*AREA + + ! The non-litter patch soil variables need to be sent to the HLM + ! This is not being sent as of this message (RGK 06-2017) + + ! --------------------------------------------------------------------- + ! The litter and CWD could be estimated from at least two methods + ! 1) after reading in the cohort data, assuming a SS turnover rate + ! 2) again assuming SS, back out the CWD and Litter flux rates into + ! the SSC, STSC and FSC pools that balance with their decomp rates + ! and then use those flux rates, to calculate the CWD and litter + ! pool sizes based on another SS model where flux out balances with + ! mortality and litter fluxes into the non-decomposing pool + ! + ! This is significant science modeling and does not have a simple + ! first hack solution. (RGK 06-2017) + ! ---------------------------------------------------------------------- + + do icwd = 1, ncwd + newpatch%cwd_ag(icwd) = 0.0_r8 + newpatch%cwd_bg(icwd) = 0.0_r8 + end do + + do ipft = 1, numpft_ed + newpatch%leaf_litter(ipft) = 0.0_r8 + newpatch%root_litter(ipft) = 0.0_r8 + end do + + return + end subroutine set_inventory_edpatch_type1 + + + ! ============================================================================================== + + subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & + patch_pointer_vec,patch_name_vec,ios) + + ! -------------------------------------------------------------------------------------------- + ! This subroutine reads in a line of an inventory cohort file (css) + ! And populates a new cohort with that information. + ! This routine specifically reads CSS files that are "Type 1" formatted + ! + ! The file formatted text, which contains 1 header line to label columns + ! and then 1 line for each cohort containing the following fields: + ! + ! FILE FORMAT: + ! time (year) year of measurement + ! patch (string) patch id string associated with this cohort + ! index (integer) cohort index + ! dbh (cm) diameter at breast height + ! height (m) height of the tree + ! pft (integer) the plant functional type index (must be consistent with param file) + ! n (/m2) The plant number density + ! bdead (kgC/plant)The dead biomass per indiv of this cohort (NOT USED) + ! balive (kgC/plant)The live biomass per indiv of this cohort (NOT USED) + ! avgRG (cm/yr?) Average Radial Growth (NOT USED) + ! -------------------------------------------------------------------------------------------- + + use EDTypesMod , only : numpft_ed + use EDGrowthFunctionsMod, only : hite + use EDGrowthFunctionsMod, only : bleaf + use EDGrowthFunctionsMod, only : bdead + use EDCohortDynamicsMod , only : create_cohort + + ! Arguments + type(ed_site_type),intent(inout), target :: csite ! current site + type(bc_in_type),intent(in) :: bc_in ! boundary conditions + integer, intent(in) :: css_file_unit ! Self explanatory + integer, intent(in) :: npatches ! number of patches + type(pp_array), intent(in) :: patch_pointer_vec(npatches) + character(len=patchname_strlen), intent(in) :: patch_name_vec(npatches) + integer,intent(out) :: ios ! Return flag + + ! Locals + real(r8) :: c_time ! Time patch was recorded + character(len=patchname_strlen) :: p_name ! The patch associated with this cohort + character(len=cohortname_strlen) :: c_name ! cohort index + real(r8) :: c_dbh ! diameter at breast height (cm) + real(r8) :: c_height ! tree height (m) + integer :: c_pft ! plant functional type index + real(r8) :: c_nplant ! plant density (/m2) + real(r8) :: c_bdead ! dead biomass (kg) + real(r8) :: c_balive ! live biomass (kg) + real(r8) :: c_avgRG ! avg radial growth (NOT USED) + integer :: cstatus ! + type(ed_patch_type), pointer :: cpatch ! current patch pointer + type(ed_cohort_type), pointer :: temp_cohort ! temporary patch (needed for allom funcs) + integer :: ipa ! patch idex + logical :: matched_patch ! check if cohort was matched w/ patch + + character(len=128),parameter :: wr_fmt = & + '(F7.1,2X,A20,2X,A20,2X,F5.2,2X,F5.2,2X,I4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' + + real(r8), parameter :: abnormal_large_nplant = 1000.0_r8 ! Used to catch bad values + real(r8), parameter :: abnormal_large_dbh = 500.0_r8 ! I've never heard of a tree > 3m + + read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_name, c_dbh, c_height, & + c_pft, c_nplant, c_bdead, c_balive, c_avgRG + + if( debug_inv) then + write(*,fmt=wr_fmt) & + c_time, p_name, c_name, c_dbh, c_height, & + c_pft, c_nplant, c_bdead, c_balive, c_avgRG + end if + + if (ios/=0) return + + ! Identify the patch based on the patch_name + matched_patch = .false. + do ipa=1,npatches + if( trim(p_name) == trim(patch_name_vec(ipa))) then + cpatch => patch_pointer_vec(ipa)%cpatch + matched_patch = .true. + end if + end do + + if(.not.matched_patch)then + write(fates_log(), *) 'could not match a cohort with a patch' + write(fates_log(),fmt=wr_fmt) & + c_time, p_name, c_name, c_dbh, c_height, & + c_pft, c_nplant, c_bdead, c_balive, c_avgRG + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Run some sanity checks on the input data + ! pft, nplant and dbh are the critical ones in this format specification + ! ------------------------------------------------------------------------------------------- + + if (c_pft > numpft_ed ) then + write(fates_log(), *) 'inventory pft: ',c_pft + write(fates_log(), *) 'An inventory cohort file specified a pft index' + write(fates_log(), *) 'greater than the maximum specified pfts ed_numpft' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (c_pft <= 0 ) then + write(fates_log(), *) 'inventory pft: ',c_pft + write(fates_log(), *) 'The inventory produced a cohort with <=0 pft index' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (c_dbh <=0 ) then + write(fates_log(), *) 'inventory dbh: ', c_dbh + write(fates_log(), *) 'The inventory produced a cohort with <= 0 dbh' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (c_dbh > abnormal_large_dbh ) then + write(fates_log(), *) 'inventory dbh: ', c_nplant + write(fates_log(), *) 'The inventory produced a cohort with very large diameter [cm]' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (c_nplant <=0 ) then + write(fates_log(), *) 'inventory nplant: ', c_nplant + write(fates_log(), *) 'The inventory produced a cohort with <= 0 density /m2' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (c_nplant > abnormal_large_nplant ) then + write(fates_log(), *) 'inventory nplant: ', c_nplant + write(fates_log(), *) 'The inventory produced a cohort with very large density /m2' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + allocate(temp_cohort) ! A temporary cohort is needed because we want to make + ! use of the allometry functions + + + + temp_cohort%pft = c_pft + temp_cohort%n = c_nplant * cpatch%area + temp_cohort%hite = Hite(temp_cohort) + temp_cohort%dbh = c_dbh + temp_cohort%canopy_trim = 1.0_r8 + temp_cohort%bdead = Bdead(temp_cohort) + temp_cohort%balive = Bleaf(temp_cohort)*(1.0_r8 + EDPftvarcon_inst%froot_leaf(c_pft) & + + EDecophyscon%sapwood_ratio(c_pft)*temp_cohort%hite) + temp_cohort%b = temp_cohort%balive + temp_cohort%bdead + + if( EDPftvarcon_inst%evergreen(c_pft) == 1) then + temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) + temp_cohort%laimemory = 0._r8 + cstatus = 2 + endif + + if( EDPftvarcon_inst%season_decid(c_pft) == 1 ) then !for dorment places + temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) !stored carbon in new seedlings. + if(csite%status == 2)then + temp_cohort%laimemory = 0.0_r8 + else + temp_cohort%laimemory = Bleaf(temp_cohort) + endif + ! reduce biomass according to size of store, this will be recovered when elaves com on. + temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory + cstatus = csite%status + endif + + if ( EDPftvarcon_inst%stress_decid(c_pft) == 1 ) then + temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) + temp_cohort%laimemory = Bleaf(temp_cohort) + temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory + cstatus = csite%dstatus + endif + + call create_cohort(cpatch, c_pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & + temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore, & + temp_cohort%laimemory, cstatus, temp_cohort%canopy_trim, 1, bc_in) + + deallocate(temp_cohort) ! get rid of temporary cohort + + return + end subroutine set_inventory_edcohort_type1 + +end module FatesInventoryInitMod