Merge pull request #821 from glemieux/satphen-ers-fix

Exact restart fix for satellite phenology mode
NGEET · Jan 14, 2022 · 858d0fd · 858d0fd
2 parents b271492 + 1188a50
commit 858d0fd
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 9 deletions.
diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90
@@ -1289,7 +1289,6 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     integer  :: ft               ! plant functional type
     integer  :: ifp              ! the number of the vegetated patch (1,2,3). In SP mode bareground patch is 0
     integer  :: patchn           ! identification number for each patch.
-    real(r8) :: canopy_leaf_area ! total amount of leaf area in the vegetated area. m2.
     real(r8) :: leaf_c           ! leaf carbon [kg]
     real(r8) :: fnrt_c           ! fineroot carbon [kg]
     real(r8) :: sapw_c           ! sapwood carbon [kg]
@@ -1318,7 +1317,6 @@ subroutine canopy_summarization( nsites, sites, bc_in )
           !zero cohort-summed variables.
           currentPatch%total_canopy_area = 0.0_r8
           currentPatch%total_tree_area = 0.0_r8
-          canopy_leaf_area = 0.0_r8
 
           !update cohort quantitie s
           currentCohort => currentPatch%shortest
@@ -1347,11 +1345,6 @@ subroutine canopy_summarization( nsites, sites, bc_in )
                 call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,&
                      currentCohort%pft,currentCohort%c_area)
              endif
-             currentCohort%treelai = tree_lai(leaf_c,             &
-                  currentCohort%pft, currentCohort%c_area, currentCohort%n, &
-                  currentCohort%canopy_layer, currentPatch%canopy_layer_tlai,currentCohort%vcmax25top )
-
-             canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
 
              if(currentCohort%canopy_layer==1)then
                 currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area
@@ -1938,6 +1931,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        bc_out(s)%dleaf_pa(:) = 0._r8
        bc_out(s)%z0m_pa(:) = 0._r8
        bc_out(s)%displa_pa(:) = 0._r8
+
        currentPatch => sites(s)%oldest_patch
        c = fcolumn(s)
        do while(associated(currentPatch))
@@ -1963,6 +1957,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
              ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic
              ! Roughness length and displacement height used by the HLM
              ! use total LAI + SAI to weight the leaft characteristic dimension
+             ! Avoid this if running in satellite phenology mode
              ! ----------------------------------------------------------------------------
 
              if (currentPatch%total_canopy_area > nearzero) then
@@ -1983,7 +1978,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
                 currentCohort => currentPatch%shortest
                 do while(associated(currentCohort))
 
-                   ! mkae sure that allometries are correct
+                   if (hlm_use_sp.eq.ifalse) then
+                   ! make sure that allometries are correct
                    call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,&
                         currentCohort%pft,currentCohort%c_area)
 
@@ -1995,6 +1991,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
                         currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
                         currentPatch%canopy_layer_tlai, currentCohort%treelai , &
                         currentCohort%vcmax25top,4)
+                   endif
 
                    total_patch_leaf_stem_area = total_patch_leaf_stem_area + &
                         (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area

diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90
@@ -11,6 +11,7 @@ module EDCohortDynamicsMod
   use FatesInterfaceTypesMod     , only : hlm_use_planthydro
   use FatesInterfaceTypesMod     , only : hlm_use_sp
   use FatesInterfaceTypesMod     , only : hlm_use_cohort_age_tracking
+  use FatesInterfaceTypesMod     , only : hlm_is_restart
   use FatesConstantsMod     , only : r8 => fates_r8
   use FatesConstantsMod     , only : fates_unset_int
   use FatesConstantsMod     , only : itrue,ifalse
@@ -1992,9 +1993,10 @@ subroutine UpdateCohortBioPhysRates(currentCohort)
        ! We assume that leaf age does not effect the specific leaf area, so the mass
        ! fractions are applicable to these rates
 
+       ipft = currentCohort%pft
+
        if(sum(frac_leaf_aclass(1:nleafage))>nearzero) then
 
-          ipft = currentCohort%pft
 
           frac_leaf_aclass(1:nleafage) =  frac_leaf_aclass(1:nleafage) / &
                 sum(frac_leaf_aclass(1:nleafage))
@@ -2011,6 +2013,13 @@ subroutine UpdateCohortBioPhysRates(currentCohort)
           currentCohort%kp25top    = sum(param_derived%kp25top(ipft,1:nleafage) * &
                 frac_leaf_aclass(1:nleafage))
 
+       elseif (hlm_use_sp .eq. itrue .and. hlm_is_restart .eq. itrue) then
+
+          currentCohort%vcmax25top = EDPftvarcon_inst%vcmax25top(ipft,1)
+          currentCohort%jmax25top  = param_derived%jmax25top(ipft,1)
+          currentCohort%tpu25top   = param_derived%tpu25top(ipft,1)
+          currentCohort%kp25top    = param_derived%kp25top(ipft,1)
+
        else
 
           currentCohort%vcmax25top = 0._r8