Merge pull request NGEET#458 from rgknox/rgknox-spawnpatch-cfix

Checks on Patch Spawning

biogeochem/EDPatchDynamicsMod.F90

@@ -329,7 +329,7 @@ subroutine spawn_patches( currentSite, bc_in)
     real(r8) :: store_c              ! storage carbon [kg]
     real(r8) :: struct_c             ! structure carbon [kg]
     real(r8) :: total_c              ! total carbon of plant [kg]
-    !---------------------------------------------------------------------
+
 
     storesmallcohort => null() ! storage of the smallest cohort for insertion routine
     storebigcohort   => null() ! storage of the largest cohort for insertion routine 
@@ -351,12 +351,17 @@ subroutine spawn_patches( currentSite, bc_in)
           call endrun(msg=errMsg(sourcefile, __LINE__))          
        end if
 
-       site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate
+       ! Only create new patches that have non-negligible amount of land
+       if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then
+          site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate
+       end if
+
        currentPatch => currentPatch%older     
 
     enddo ! end loop over patches. sum area disturbed for all patches. 
 
-    if (site_areadis > 0.0_r8) then  
+    if (site_areadis > nearzero) then  
+
        cwd_ag_local = 0.0_r8
        cwd_bg_local = 0.0_r8
        leaf_litter_local = 0.0_r8
@@ -378,6 +383,8 @@ subroutine spawn_patches( currentSite, bc_in)
           ! This is the amount of patch area that is disturbed, and donated by the donor
           patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate
 
+          if (patch_site_areadis > nearzero) then
+
           call average_patch_properties(currentPatch, new_patch, patch_site_areadis)
 
           if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. &
@@ -682,10 +689,6 @@ subroutine spawn_patches( currentSite, bc_in)
           enddo ! currentCohort 
           call sort_cohorts(currentPatch)
 
-          !zero disturbance accumulators
-          currentPatch%disturbance_rate  = 0._r8
-          currentPatch%disturbance_rates = 0._r8
-
           !update area of donor patch
           currentPatch%area = currentPatch%area - patch_site_areadis
 
@@ -698,6 +701,12 @@ subroutine spawn_patches( currentSite, bc_in)
           call terminate_cohorts(currentSite, currentPatch, 2)
           call sort_cohorts(currentPatch)
 
+          end if    ! if (patch_site_areadis > nearzero) then
+
+          !zero disturbance rate trackers
+          currentPatch%disturbance_rate  = 0._r8
+          currentPatch%disturbance_rates = 0._r8
+
           currentPatch => currentPatch%younger
 
        enddo ! currentPatch patch loop. 
@@ -723,9 +732,11 @@ subroutine spawn_patches( currentSite, bc_in)
 
     endif !end new_patch area 
 
+
     call check_patch_area(currentSite)
     call set_patchno(currentSite)
 
+    return
   end subroutine spawn_patches
 
   ! ============================================================================

main/EDMainMod.F90

@@ -538,6 +538,12 @@ subroutine ed_total_balance_check (currentSite, call_index )
     real(r8) :: error           ! How much carbon did we gain or lose (should be zero!) 
     real(r8) :: error_frac      ! Error as a fraction of total biomass
     real(r8) :: net_flux        ! Difference between recorded fluxes in and out. KgC/site
+    real(r8) :: leaf_c
+    real(r8) :: fnrt_c
+    real(r8) :: sapw_c
+    real(r8) :: store_c
+    real(r8) :: struct_c
+    real(r8) :: repro_c
 
     ! nb. There is no time associated with these variables 
     ! because this routine can be called between any two 
@@ -579,7 +585,7 @@ subroutine ed_total_balance_check (currentSite, call_index )
     !               burned_litter * new_patch%area !kG/site/day
     ! -----------------------------------------------------------------------------------
 
-    if ( error_frac > 10e-6 ) then
+    if ( error_frac > 10e-6_r8 ) then
        write(fates_log(),*) 'carbon balance error detected'
        write(fates_log(),*) 'error fraction relative to biomass stock:',error_frac
        write(fates_log(),*) 'call index: ',call_index
@@ -593,33 +599,40 @@ subroutine ed_total_balance_check (currentSite, call_index )
        write(fates_log(),*) 'seeds',seed_stock
        write(fates_log(),*) 'previous total',currentSite%old_stock  
 
-       if(debug)then
+       write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon
+
+       ! 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
+
           change_in_stock = 0.0_r8
           biomass_stock   = 0.0_r8
           litter_stock    = 0.0_r8
 
           seed_stock   =  sum(currentSite%seed_bank)*AREA
-          currentPatch => currentSite%oldest_patch 
+          currentPatch => currentSite%oldest_patch
           do while(associated(currentPatch))
              write(fates_log(),*) '---------------------------------------'
              write(fates_log(),*) currentPatch%area , sum(currentPatch%cwd_ag), sum(currentPatch%cwd_bg)
              write(fates_log(),*) sum(currentPatch%leaf_litter),sum(currentPatch%root_litter)
              write(fates_log(),*)'---'
              currentCohort => currentPatch%tallest
              do while(associated(currentCohort))
-                write(fates_log(),*) 'structure: ',currentCohort%prt%GetState(struct_organ,all_carbon_elements)
-                write(fates_log(),*) 'storage: ',currentCohort%prt%GetState(store_organ,all_carbon_elements)
+                write(fates_log(),*) 'pft: ',currentCohort%pft
+                write(fates_log(),*) 'dbh: ',currentCohort%dbh
+                leaf_c   = currentCohort%prt%GetState(leaf_organ,all_carbon_elements)
+                struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements)
+                store_c  = currentCohort%prt%GetState(store_organ,all_carbon_elements)
+                fnrt_c   = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements)
+                repro_c  = currentCohort%prt%GetState(repro_organ,all_carbon_elements)
+                sapw_c   = currentCohort%prt%GetState(sapw_organ,all_carbon_elements)
+
+                write(fates_log(),*) 'lc: ',leaf_c,' dc: ',struct_c,' stc: ',store_c
+                write(fates_log(),*) 'fc: ',fnrt_c,' rc: ',repro_c,' sac: ',sapw_c
                 write(fates_log(),*) 'N plant: ',currentCohort%n
-                currentCohort => currentCohort%shorter;
-             enddo !end cohort loop 
+                currentCohort => currentCohort%shorter
+             enddo !end cohort loop                                                                                                   
              currentPatch => currentPatch%younger
-          enddo !end patch loop
-       end if
-
-       write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon
-
-       ! 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
+          enddo !end patch loop                                                                                                      
           write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day
           call endrun(msg=errMsg(sourcefile, __LINE__))
        end if

main/EDTypesMod.F90

@@ -145,11 +145,11 @@ module EDTypesMod
   integer , parameter :: N_HITE_BINS          = 60         ! no. of hite bins used to distribute LAI
 
   ! COHORT TERMINATION
-  real(r8), parameter :: min_npm2       = 1.0E-8_r8  ! minimum cohort number density per m2 before termination
-  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
+  real(r8), parameter :: min_npm2       = 1.0E-7_r8               ! minimum cohort number density per m2 before termination
+  real(r8), parameter :: min_patch_area = 0.01_r8                 ! smallest allowable patch area before termination
+  real(r8), parameter :: min_nppatch    = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area)
+  real(r8), parameter :: min_n_safemath = 1.0E-12_r8              ! in some cases, we want to immediately remove super small
+                                                                  ! number densities of cohorts to prevent FPEs
 
   character*4 yearchar