Merge pull request #282 from rgknox/rgknox-cbalance-check

carbon balance fix

biogeochem/EDPhysiologyMod.F90

@@ -1006,6 +1006,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
     if ( DEBUG ) write(fates_log(),*) 'EDPhys dbstoredt I ',currentCohort%dbstoredt
 
     currentCohort%seed_prod = (1.0_r8 - gr_fract) * currentCohort%carbon_balance
+
     if (abs(currentCohort%npp_acc_hold-(currentCohort%dbalivedt+currentCohort%dbdeaddt+currentCohort%dbstoredt+ &
          currentCohort%seed_prod+currentCohort%md)) > 0.0000000001_r8)then
        write(fates_log(),*) 'error in carbon check growth derivs',currentCohort%npp_acc_hold- &
@@ -1041,7 +1042,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
   end subroutine Growth_Derivatives
 
   ! ============================================================================
-  subroutine recruitment( t, currentSite, currentPatch, bc_in )
+  subroutine recruitment( currentSite, currentPatch, bc_in )
     !
     ! !DESCRIPTION:
     ! spawn new cohorts of juveniles of each PFT             
@@ -1051,7 +1052,6 @@ subroutine recruitment( t, currentSite, currentPatch, bc_in )
     use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys
     !
     ! !ARGUMENTS    
-    integer, intent(in) :: t
     type(ed_site_type), intent(inout), target  :: currentSite
     type(ed_patch_type), intent(inout), pointer :: currentPatch
     type(bc_in_type), intent(in)                :: bc_in
@@ -1076,20 +1076,14 @@ subroutine recruitment( t, currentSite, currentPatch, bc_in )
             + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite)
        temp_cohort%bstore      = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) &
             + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))
+
        if (hlm_use_ed_prescribed_phys .eq. ifalse) then
           temp_cohort%n           = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day &
                / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore)
        else
           ! prescribed recruitment rates. number per sq. meter per year
           temp_cohort%n        = currentPatch%area * EDPftvarcon_inst%prescribed_recruitment(ft) * hlm_freq_day
        endif
-
-       if (t == 1)then
-          write(fates_log(),*) 'filling in cohorts where there are none left; this will break carbon balance', &
-               currentPatch%patchno,currentPatch%area
-          temp_cohort%n = 0.1_r8*currentPatch%area
-          write(fates_log(),*) 'cohort n',ft,temp_cohort%n
-       endif
 
        temp_cohort%laimemory = 0.0_r8     
        if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == 1)then
@@ -1198,10 +1192,11 @@ subroutine CWD_Input( currentSite, currentPatch)
 
           dead_n_natural = dead_n - dead_n_dlogging - dead_n_ilogging
 
+
           currentPatch%leaf_litter_in(pft) = currentPatch%leaf_litter_in(pft) + &
-               (currentCohort%bl+currentCohort%leaf_litter/hlm_freq_day)* dead_n          
+               currentCohort%bl * dead_n          
           currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + &
-                (currentCohort%br+currentCohort%bstore)     * dead_n
+               (currentCohort%br+currentCohort%bstore)     * dead_n
 
           ! Update diagnostics that track resource management
           currentSite%resources_management%delta_litter_stock  = &
@@ -1223,12 +1218,17 @@ subroutine CWD_Input( currentSite, currentPatch)
 
              ! Send AGB component of boles from non direct-logging activities to AGB litter pool
              if (c==ncwd) then
+
                 currentPatch%cwd_AG_in(c) = currentPatch%cwd_AG_in(c) + (currentCohort%bdead+currentCohort%bsw) * &
-                      SF_val_CWD_frac(c) * (dead_n_natural+dead_n_ilogging)  * &
-                      EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)
+                     SF_val_CWD_frac(c) * (dead_n_natural+dead_n_ilogging)  * &
+                     EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)
 
-
              else
+
+                currentPatch%cwd_AG_in(c) = currentPatch%cwd_AG_in(c) + (currentCohort%bdead+currentCohort%bsw) * &
+                     SF_val_CWD_frac(c) * dead_n  * &
+                     EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)
+
                 ! Send AGB component of boles from direct-logging activities to export/harvest pool
                 ! Generate trunk product (kgC/day/site)
                 trunk_product = (currentCohort%bdead+currentCohort%bsw) * &

main/EDMainMod.F90

@@ -14,6 +14,7 @@ module EDMainMod
   use FatesInterfaceMod        , only : hlm_current_month
   use FatesInterfaceMod        , only : hlm_current_day 
   use FatesInterfaceMod        , only : hlm_use_planthydro 
+  use FatesInterfaceMod        , only : hlm_reference_date
   use FatesInterfaceMod        , only : hlm_use_ed_st3 
   use FatesInterfaceMod        , only : bc_in_type
   use FatesInterfaceMod        , only : hlm_masterproc
@@ -145,7 +146,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in)
        do while (associated(currentPatch))                 
 
           ! adds small cohort of each PFT
-          call recruitment(0, currentSite, currentPatch, bc_in)
+          call recruitment(currentSite, currentPatch, bc_in)
 
           currentPatch => currentPatch%younger
        enddo
@@ -514,7 +515,7 @@ subroutine ed_total_balance_check (currentSite, call_index )
     !            burn fractions:  
     !               currentSite%flux_out = currentSite%flux_out + &
     !               burned_litter * new_patch%area !kG/site/day
-    !            
+    ! -----------------------------------------------------------------------------------
 
     if ( error_frac > 10e-6 ) then
        write(fates_log(),*) 'carbon balance error detected'
@@ -552,11 +553,13 @@ subroutine ed_total_balance_check (currentSite, call_index )
        end if
 
        write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon
-
-       ! Only abort if this is not the first call
-       if(abs(total_stock - currentSite%old_stock)>1.e-15_r8)then
-   !       call endrun(msg=errMsg(sourcefile, __LINE__))
+
+       ! If this is the first day of simulation, carbon balance reports but does not end the run
+       if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then
+          write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day
+          call endrun(msg=errMsg(sourcefile, __LINE__))
        end if
+
     endif
 
     currentSite%flux_in   = 0.0_r8