Initial carbon linkages between ED and soil BGC

Merge andre-ed-v0.1.0-buggy-science into master initial implementation of the carbon linkages from ED to the CLM soil biogeochemistry cycle. I originally did this in the CLM SVN, which Ben pulled into the NGT github on the "buggy science" branch. The purpose of this code is to take all the litter that flows out of ED, aggregated onto column-level variables and disaggregated vertically and into litter decomposability pools, and then send this to the BGC code that is shared with the rest of the model. The code behavior should be identical to the master for the variables that are output in master, but there is now a much larger set of carbon cycle variables that were not output before. New output variables include soil and litter carbon pools, combined ED + BGC pools (TOTECOSYSC), and a bunch of fluxes (heterotrophic respiration, NEP, NBP, etc). I've also included some extra diagnostics, including high-frequency NPP on the column level (current standard NPP is daily-averaged on patch level, which I've left as-is). For carbon balance, I am now tracking all fluxes into the ED pools and the BGC pools, separately and together, over each day (since ED fluxes are resolved on the daily timestep and CLM fluxes are resolved on the half-hourly timestep). Carbon balance errors are not trivially small, so instead of implementing a tolerance-checking crasher I am just writing the errors to history for now. When we get these smaller we can replace them with logic that crashes when a tolerance is exceeded. I find some of the behavior hard to check in the history files because of issues #24 and #26, but otherwise things seem reasonable. I have avoided crashing bugs associated with issue #30 but not actually tried to solve that. Code Review: C Koven, R Knox Closes: #35 Test suite: ed - lawrencium pass; ed - yellowstone intel, pgi, gnu Test baseline: d8a9ee5 Test namelist changes: yes Test answer changes: yes Test summary: expected fails in f19 and f09 restart, answer changing and nl changes as expected.
NGEET · Mar 23, 2016 · 3da8a77 · 3da8a77
2 parents 2c3e52d + 4230198
commit 3da8a77
Show file tree

Hide file tree

Showing 6 changed files with 1,155 additions and 136 deletions.
diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90
@@ -733,7 +733,7 @@ subroutine fuse_cohorts(patchptr)
                                   nextc%n*nextc%year_net_uptake(i))/newn                
                             endif
                          enddo
-
+                         
                          currentCohort%n = newn     
                          !remove fused cohort from the list
                          nextc%taller%shorter => nextnextc        

diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -200,8 +200,11 @@ subroutine spawn_patches( currentSite )
 
     ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. 
     currentPatch => currentSite%youngest_patch
+
+    ! zero site-level fire fluxes
     currentSite%cwd_ag_burned       = 0.0_r8
     currentSite%leaf_litter_burned  = 0.0_r8
+    currentSite%total_burn_flux_to_atm = 0.0_r8    
 
     site_areadis = 0.0_r8
     do while(associated(currentPatch))
@@ -548,12 +551,14 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
           burned_litter = new_patch%cwd_ag(c) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(c+1) !kG/m2/day
           new_patch%cwd_ag(c) = new_patch%cwd_ag(c) - burned_litter
           currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day
+          currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day
        enddo
 
        do p = 1,numpft_ed
           burned_litter = new_patch%leaf_litter(p) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(dg_sf)
           new_patch%leaf_litter(p) = new_patch%leaf_litter(p) - burned_litter
-          currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/dat
+          currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day
+          currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day
       enddo
 
        !************************************/     
@@ -615,6 +620,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
                      SF_val_CWD_frac(c) * bstem * currentCohort%cfa
                 currentSite%flux_out  = currentSite%flux_out + dead_tree_density * &
                      AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa
+                currentSite%total_burn_flux_to_atm  = currentSite%total_burn_flux_to_atm + dead_tree_density * &
+                     AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa
 
              enddo
 
@@ -625,6 +632,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
                      dead_tree_density * currentCohort%bl * currentCohort%cfa
                 currentSite%flux_out  = currentSite%flux_out + &
                      dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa
+                currentSite%total_burn_flux_to_atm  = currentSite%total_burn_flux_to_atm + &
+                     dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa
 
              enddo
 
@@ -655,6 +664,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
              !KgC/gridcell/day
              currentSite%flux_out = currentSite%flux_out + burned_leaves * currentCohort%n * &
                   patch_site_areadis/currentPatch%area * AREA 
+             currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm+ burned_leaves * currentCohort%n * &
+                  patch_site_areadis/currentPatch%area * AREA 
 
           endif
           currentCohort%cfa = 0.0_r8        
@@ -1113,7 +1124,7 @@ subroutine fuse_patches( csite )
        if(nopatches > maxpatch)then
           iterate = 1
           profiletol = profiletol * 1.1_r8
-          write(iulog,*) 'maxpatch exceeded, triggering patch fusion iteration.',profiletol,nopatches
+
           !---------------------------------------------------------------------!
           ! Making profile tolerance larger means that more fusion will happen  !
           !---------------------------------------------------------------------!        

diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90
@@ -125,15 +125,14 @@ subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, wa
        currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - currentPatch%root_litter_out(p) 
     enddo
 
-    currentPatch%leaf_litter_in(:)  = 0.0_r8
-    currentPatch%root_litter_in(:)  = 0.0_r8
-    currentPatch%leaf_litter_out(:) = 0.0_r8
-    currentPatch%root_litter_out(:) = 0.0_r8
-
-    currentPatch%CWD_AG_in(:)       = 0.0_r8
-    currentPatch%cwd_bg_in(:)       = 0.0_r8
-    currentPatch%CWD_AG_out(:)      = 0.0_r8
-    currentPatch%cwd_bg_out(:)      = 0.0_r8
+    ! currentPatch%leaf_litter_in(:)  = 0.0_r8
+    ! currentPatch%root_litter_in(:)  = 0.0_r8
+    ! currentPatch%leaf_litter_out(:) = 0.0_r8
+    ! currentPatch%root_litter_out(:) = 0.0_r8
+    ! currentPatch%CWD_AG_in(:)       = 0.0_r8
+    ! currentPatch%cwd_bg_in(:)       = 0.0_r8
+    ! currentPatch%CWD_AG_out(:)      = 0.0_r8
+    ! currentPatch%cwd_bg_out(:)      = 0.0_r8
 
   end subroutine non_canopy_derivs
 
@@ -636,6 +635,8 @@ subroutine seeds_in( cp_pnt )
     currentSite  => currentPatch%siteptr
 
     currentPatch%seeds_in(:) = 0.0_r8
+    currentPatch%seed_rain_flux(:) = 0.0_r8
+
     currentCohort => currentPatch%tallest
     do while (associated(currentCohort))
        p = currentCohort%pft
@@ -649,6 +650,7 @@ subroutine seeds_in( cp_pnt )
        if (EXTERNAL_RECRUITMENT == 1) then !external seed rain - needed to prevent extinction  
           do p = 1,numpft_ed
            currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + EDecophyscon%seed_rain(p) !KgC/m2/year
+           currentPatch%seed_rain_flux(p) = currentPatch%seed_rain_flux(p) + EDecophyscon%seed_rain(p) !KgC/m2/year
           enddo
        endif
        currentPatch => currentPatch%younger
@@ -1251,4 +1253,6 @@ subroutine cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_i
 
   end subroutine cwd_out
 
+
+
 end module EDPhysiologyMod
diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90
@@ -857,10 +857,10 @@ subroutine ED_Norman_Radiation (bounds, &
                             enddo
                          enddo
                          if (lai_change(1,2,1).gt.0.0.and.lai_change(1,2,2).gt.0.0)then
-                             write(iulog,*) 'lai_change(1,2,12)',lai_change(1,2,1:4)
+!                             write(iulog,*) '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(iulog,*) ' lai_change (1,2,23)',lai_change(1,2,1:4)
+!                             write(iulog,*) ' 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