Merge pull request #1448 from glemieux/nbp-gpp-restart-fix

`gpp_site` restart fix

Author: glemieux

biogeochem/EDPatchDynamicsMod.F90

@@ -1966,9 +1966,13 @@ subroutine TransLitterNewPatch(currentSite,        &
        curr_litt  => currentPatch%litter(el)
        new_litt  => newPatch%litter(el)
 
-       ! Distribute the fragmentation litter flux rates. This is only used for diagnostics
-       ! at this point.  Litter fragmentation has already been passed to the output
-       ! boundary flux arrays.
+       ! Distribute the fragmentation litter flux rates. The mean site-level
+       ! flux rate must be preserved, so when we create new patches
+       ! from disturbance, we must area weight the contributions of the
+       ! donor patches.  This is because the host model will call
+       ! FatesSoilBGCFluxMod:FluxIntoLitterPools() which uses these
+       ! litt%<>_frac() arrays to fill site level output fluxes, and
+       ! this is called over the next day on the model timestep.       
 
        do c = 1,ncwd 
           new_litt%ag_cwd_frag(c) = new_litt%ag_cwd_frag(c) + &

main/EDMainMod.F90

@@ -426,9 +426,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
 
     current_fates_landuse_state_vector = currentSite%get_current_landuse_statevector()
 
-    ! Clear site GPP and AR passing to HLM
-    bc_out%gpp_site = 0._r8
-    bc_out%ar_site = 0._r8
+    
 
     ! Patch level biomass are required for C-based harvest
     call get_harvestable_carbon(currentSite, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c)
@@ -645,20 +643,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
 
           currentCohort%npp_acc_hold  = currentCohort%npp_acc_hold - &
                currentCohort%resp_excess_hold*real( hlm_days_per_year,r8)
-
-          ! Passing gpp_acc_hold to HLM 
-          bc_out%gpp_site = bc_out%gpp_site + currentCohort%gpp_acc_hold * &
-               AREA_INV * currentCohort%n / real( hlm_days_per_year,r8) / sec_per_day
-          bc_out%ar_site = bc_out%ar_site + (currentCohort%resp_m_acc_hold + &
-               currentCohort%resp_g_acc_hold + currentCohort%resp_excess_hold*real(hlm_days_per_year,r8) ) * & 
-               AREA_INV * currentCohort%n / real( hlm_days_per_year,r8) / sec_per_day
 
           ! Update the mass balance tracking for the daily nutrient uptake flux
           ! Then zero out the daily uptakes, they have been used
 
-          ! -----------------------------------------------------------------------------
-
-
 
           call EffluxIntoLitterPools(currentSite, currentPatch, currentCohort, bc_in )
 
@@ -693,7 +681,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
                currentCohort%resp_m_acc*currentCohort%n + &
                currentCohort%resp_excess_hold*currentCohort%n + &
                currentCohort%resp_g_acc_hold*currentCohort%n/real( hlm_days_per_year,r8)
-
+          
+          
+          
           call currentCohort%prt%CheckMassConservation(ft,5)
 
           ! Update the leaf biophysical rates based on proportion of leaf
@@ -851,13 +841,23 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting )
     real(r8) :: total_stock  ! dummy variable for receiving from sitemassstock
     !-----------------------------------------------------------------------
 
+    site_cmass => currentSite%mass_balance(element_pos(carbon12_element))
+    
     ! check patch order (set second argument to true)
     if (debug) then
        call set_patchno(currentSite,.true.,1)
     end if
+
+    ! Pass site-level mass fluxes to output boundary conditions
+    ! [kg/site/day] * [site/m2 day/sec] = [kgC/m2/s]
+    bc_out%gpp_site = site_cmass%gpp_acc * area_inv / sec_per_day
+    bc_out%ar_site  = site_cmass%aresp_acc * area_inv / sec_per_day
 
     if(hlm_use_sp.eq.ifalse .and. (.not.is_restarting))then
-      call canopy_spread(currentSite)
+       call canopy_spread(currentSite)
+    else
+       site_cmass%gpp_acc = 0._r8
+       site_cmass%aresp_acc = 0._r8
     end if
 
     call TotalBalanceCheck(currentSite,6)
@@ -925,11 +925,13 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting )
     bc_out%seed_c_si = bc_out%seed_c_si * g_per_kg * AREA_INV
 
     ! Set boundary condition to HLM for carbon loss to atm from fires and grazing
-    ! [kgC/ha/day]*[m2/ha]*[day/s] = [kg/m2/s] 
-    site_cmass => currentSite%mass_balance(element_pos(carbon12_element))
+    ! [kgC/ha/day]*[ha/m2]*[day/s] = [kg/m2/s] 
+    
     bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * ha_per_m2 * days_per_sec
     bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * ha_per_m2 * days_per_sec
 
+    
+
   end subroutine ed_update_site
 
   !-------------------------------------------------------------------------------!
@@ -1178,14 +1180,6 @@ subroutine bypass_dynamics(currentSite, bc_out)
           ! Shouldn't need to zero any nutrient fluxes
           ! as they should just be zero, no uptake
           ! in ST3 mode.
-
-          ! Passing 
-          bc_out%gpp_site = bc_out%gpp_site + currentCohort%gpp_acc_hold * &
-               AREA_INV * currentCohort%n / real( hlm_days_per_year,r8) / sec_per_day
-          bc_out%ar_site = bc_out%ar_site + (currentCohort%resp_m_acc_hold + &
-               currentCohort%resp_g_acc_hold +  &
-               currentCohort%resp_excess_hold*real( hlm_days_per_year,r8)) * & 
-               AREA_INV * currentCohort%n / real( hlm_days_per_year,r8) / sec_per_day
 
           currentCohort => currentCohort%taller
        enddo

main/FatesInterfaceMod.F90

@@ -159,7 +159,7 @@ module FatesInterfaceMod
       ! instance is fine.
 
       type(bc_pconst_type) :: bc_pconst
-      
+
 
    end type fates_interface_type
 

main/FatesInterfaceTypesMod.F90

@@ -861,6 +861,4 @@ subroutine ZeroBCOutCarbonFluxes(bc_out)
 
   end subroutine ZeroBCOutCarbonFluxes
 
-   
-       
 end module FatesInterfaceTypesMod

main/FatesRestartInterfaceMod.F90

@@ -49,8 +49,9 @@ module FatesRestartInterfaceMod
   use EDTypesMod,              only : area
   use EDTypesMod,              only : set_patchno
   use EDParamsMod,             only : nlevleaf
-  use PRTGenericMod,           only : prt_global
+  use PRTGenericMod,           only : carbon12_element
   use PRTGenericMod,           only : num_elements
+  use PRTGenericMod,           only : element_pos
   use FatesRunningMeanMod,     only : rmean_type
   use FatesRunningMeanMod,     only : ema_lpa
   use FatesRadiationMemMod,    only : num_swb,norman_solver,twostr_solver
@@ -111,6 +112,9 @@ module FatesRestartInterfaceMod
   integer :: ir_area_bareground_si
   integer :: ir_snow_depth_si
   integer :: ir_landuse_config_si
+  integer :: ir_gpp_acc_si
+  integer :: ir_aresp_acc_si
+  
   integer :: ir_ncohort_pa
   integer :: ir_canopy_layer_co
   integer :: ir_canopy_layer_yesterday_co
@@ -341,7 +345,8 @@ module FatesRestartInterfaceMod
 
   ! The number of variable dim/kind types we have defined (static)
   integer, parameter, public :: fates_restart_num_dimensions = 2   !(cohort,column)
-  integer, parameter, public :: fates_restart_num_dim_kinds = 4    !(cohort-int,cohort-r8,site-int,site-r8)
+  integer, parameter, public :: fates_restart_num_dim_kinds = 4    !(cohort-int,cohort-r8,
+                                                                   ! site-int,site-r8)
 
   ! integer constants for storing logical data
   integer, parameter, public :: old_cohort = 0
@@ -769,6 +774,16 @@ subroutine define_restart_vars(this, initialize_variables)
          units='kgC/m2', flushval = flushzero, &
          hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_landuse_config_si )
 
+    call this%set_restart_var(vname='fates_massbal_gpp', vtype=site_r8, &
+         long_name='accumulated gpp over previous day cycle', &
+         units='kgC/m2/s', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gpp_acc_si )
+
+    call this%set_restart_var(vname='fates_massbal_ar', vtype=site_r8, &
+         long_name='accumulated autotrophic respiration over previous day cycle', &
+         units='kgC/m2/s', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_aresp_acc_si )
+    
     ! -----------------------------------------------------------------------------------
     ! Variables stored within cohort vectors
     ! Note: Some of these are multi-dimensional variables in the patch/site dimension
@@ -2229,6 +2244,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
 
     integer  :: ft               ! functional type index
     integer  :: el               ! element loop index
+    integer  :: c_el             ! element loop index for carbon12
     integer  :: ilyr             ! soil layer index
     integer  :: nlevsoil         ! total soil layers in patch of interest
     integer  :: k,j,i            ! indices to the radiation matrix
@@ -2566,7 +2582,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
              end do
           end if
 
-
+          c_el = element_pos(carbon12_element)
+          this%rvars(ir_gpp_acc_si)%r81d(io_idx_si) = sites(s)%mass_balance(c_el)%gpp_acc
+          this%rvars(ir_aresp_acc_si)%r81d(io_idx_si) = sites(s)%mass_balance(c_el)%aresp_acc
 
           ! canopy spread term
           rio_spread_si(io_idx_si)   = sites(s)%spread
@@ -3281,6 +3299,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
      integer  :: patchespersite   ! number of patches per site
      integer  :: cohortsperpatch  ! number of cohorts per patch
      integer  :: el               ! loop counter for elements
+     integer  :: c_el             ! loop counter for carbon12
      integer  :: nlevsoil         ! number of soil layers
      integer  :: ilyr             ! soil layer loop counter
      integer  :: iscpf            ! multiplex loop counter for size x pft
@@ -3599,6 +3618,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
 
              end do
           end if
+
+          c_el = element_pos(carbon12_element)
+          sites(s)%mass_balance(c_el)%gpp_acc   = this%rvars(ir_gpp_acc_si)%r81d(io_idx_si)
+          sites(s)%mass_balance(c_el)%aresp_acc = this%rvars(ir_aresp_acc_si)%r81d(io_idx_si)
+          
 
           sites(s)%spread = rio_spread_si(io_idx_si)