diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90
index c68e6a1efb..e9167c2513 100644
--- a/biogeochem/EDCanopyStructureMod.F90
+++ b/biogeochem/EDCanopyStructureMod.F90
@@ -2266,6 +2266,13 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area)
    
    ! Number of actual vegetation layers in this cohort's crown
    currentCohort%nv =  count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
+
+   if( currentCohort%nv .ne. minloc(dlower_vai, DIM=1, MASK=(dlower_vai>(currentCohort%treelai+currentCohort%treesai))) ) then
+      write(fates_log(),*) 'We use two methods of finding maximum leaf layers, and they are not equivalent'
+      write(fates_log(),*) 'count method:',currentCohort%nv
+      write(fates_log(),*) 'minloc method:',minloc(dlower_vai, DIM=1, MASK=(dlower_vai>(currentCohort%treelai+currentCohort%treesai)))
+      call endrun(msg=errMsg(sourcefile, __LINE__))
+   end if
    
   end subroutine UpdateCohortLAI
   
diff --git a/radiation/FatesTwoStreamUtilsMod.F90 b/radiation/FatesTwoStreamUtilsMod.F90
index 3abe0ce849..0d7ef85419 100644
--- a/radiation/FatesTwoStreamUtilsMod.F90
+++ b/radiation/FatesTwoStreamUtilsMod.F90
@@ -78,6 +78,20 @@ subroutine FatesConstructRadElements(site,fcansno_pa,coszen_pa)
     integer :: max_elements     ! Maximum number of scattering elements on the site
     integer :: n_scr            ! The size of the scratch arrays
     logical :: allocate_scratch ! Whether to re-allocate the scratch arrays
+    integer  :: icolmax         ! Column index for each layer with largest area footprint
+    real(r8) :: areamax         ! The area footprint of the largest column
+    
+    ! its possible that there is more horizontal area taken up by the cohorts
+    ! than there is ground, which is simply a result numerical and algorithmic
+    ! imprecision in the rest of FATES. If this is true, then we need
+    ! to somehow force the total area of our scattering elements to be exactly
+    ! 1 and not slightly more. One way is to just chop off some area of the
+    ! largest scattering element (simple way), the other is to chop off that
+    ! area but also increase the LAI+SAI. The latter is conservative, but
+    ! could create indexing problems when transfering fluxes back into FATES arrays
+    
+    logical, parameter :: do_simple_area_correct = .true.
+    
     
     ! These parameters are not used yet
     !real(r8) :: max_vai_diff_per_elem ! The maximum vai difference in any element
@@ -119,6 +133,13 @@ subroutine FatesConstructRadElements(site,fcansno_pa,coszen_pa)
          ! an air element is needed for all the non
          ! occupied space, even if the canopy_open_frac
          ! is zero.
+         ! If the area of the elements does not match
+         ! the area of the canopy space within 1.e-7_r8
+         ! then we either add the space in the form of air
+         ! or we compress the space by literally squeezing
+         ! the elements (which consequently increases their
+         ! LAI and SAI to conserve area)
+        
             
          if(patch%total_canopy_area>nearzero)then
             canopy_frac(:) = 0._r8
@@ -131,7 +152,8 @@ subroutine FatesConstructRadElements(site,fcansno_pa,coszen_pa)
          else
             canopy_frac(:) = 0._r8
          end if
-         
+
+         ! Add the air element if the canopy area is not filled
          do ican = 1,patch%ncl_p
             if( (1._r8-canopy_frac(ican))>area_err_thresh ) then
                n_col(ican) = n_col(ican) + 1
@@ -229,25 +251,57 @@ subroutine FatesConstructRadElements(site,fcansno_pa,coszen_pa)
                twostr%scelg(ican,n_col(ican))%sai  = 0._r8
             end if
 
+            ! Check to see if any of these layers are extremely over-full and stop the
+            ! model if so. This should not happen and would most likely be an issue
+            ! with the canopy ppa promotion/demotion logic
+            ! This is more of a sanity check, so we use a 1% threshold
+            if(debug)then
+               if_very_overfull: if( (canopy_frac(ican)-1._r8)>0.01_r8 ) then
+                  write(fates_log(),*) 'One of the fates canopy layers takes up'
+                  write(fates_log(),*) 'more than 100% of the area footprint, exceeding a'
+                  write(fates_log(),*) 'precision threshold of 0.01'
+                  write(fates_log(),*) 'Aborting'
+                  write(fates_log(),*) 'canopy layer: ',ican,' canopy_frac:',canopy_frac(ican)
+                  call endrun(msg=errMsg(sourcefile, __LINE__))
+               end if if_very_overfull
+            end if
             ! If the layer is overfull, remove some from area from
-            ! the first element that is 10x larger than the threshold
+            ! the element with the largest footprint
 
             if_overfull: if( (canopy_frac(ican)-1._r8)>area_err_thresh ) then
+
+               ! First find the element with the largest footprint
+               icolmax = -1
+               areamax = 0
                do icol = 1,n_col(ican)
-                  if(twostr%scelg(ican,icol)%area > 10._r8*(canopy_frac(ican)-1._r8))then
-                      area_ratio = (twostr%scelg(ican,icol)%area + (1._r8-canopy_frac(ican)))/twostr%scelg(ican,icol)%area
-                     twostr%scelg(ican,icol)%area = twostr%scelg(ican,icol)%area * area_ratio
-                     twostr%scelg(ican,icol)%lai  = twostr%scelg(ican,icol)%lai / area_ratio
-                     twostr%scelg(ican,icol)%sai  = twostr%scelg(ican,icol)%sai / area_ratio
-                     canopy_frac(ican) = 1.0_r8
-                     exit if_overfull
+                  if(twostr%scelg(ican,icol)%area>areamax) then
+                     icolmax = icol
+                     areamax = twostr%scelg(ican,icol)%area
                   end if
                end do
                
-               !write(fates_log(),*) 'overfull areas'
-               !twostr%cosz = coszen_pa(ifp)
-               !     call twostr%Dump(1,lat=site%lat,lon=site%lon)
-               !     call endrun(msg=errMsg(sourcefile, __LINE__))
+               ! Test out a simpler way to correct area errors
+               if(do_simple_area_correct) then
+                  if(debug) then
+                     if((canopy_frac(ican)-1._r8)>0.5_r8*twostr%scelg(ican,icolmax)%area)then
+                        write(fates_log(),*) 'An area correction is being applied where '
+                        write(fates_log(),*) 'the correction is greater than 50% of the area of the largest donor.'
+                        write(fates_log(),*) 'This will have to large of an impact on the donor and is not representative'
+                        write(fates_log(),*) 'of the actual composition of the canopy'
+                        write(fates_log(),*) 'Aborting'
+                        write(fates_log(),*) 'canopy layer: ',ican,' canopy_frac:',canopy_frac(ican)
+                        write(fates_log(),*) 'existing donor area',twostr%scelg(ican,icolmax)%area
+                        call endrun(msg=errMsg(sourcefile, __LINE__))
+                     end if
+                  end if
+                  twostr%scelg(ican,icolmax)%area = twostr%scelg(ican,icolmax)%area - (canopy_frac(ican)-1._r8)
+               else
+                  area_ratio = (twostr%scelg(ican,icolmax)%area + (1._r8-canopy_frac(ican)))/twostr%scelg(ican,icolmax)%area
+                  twostr%scelg(ican,icolmax)%area = twostr%scelg(ican,icolmax)%area * area_ratio
+                  twostr%scelg(ican,icolmax)%lai  = twostr%scelg(ican,icolmax)%lai / area_ratio
+                  twostr%scelg(ican,icolmax)%sai  = twostr%scelg(ican,icolmax)%sai / area_ratio
+               end if
+
             end if if_overfull
 
          end do