Merge pull request #2 from rgknox/fire_mortality

conflict resolutions and fixes to canopy layer diagnostics

biogeochem/EDCanopyStructureMod.F90

@@ -504,10 +504,10 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr)
 
                   allocate(copyc)
                   call InitPRTCohort(copyc)
-                  call copy_cohort(currentCohort, copyc)
                   if( hlm_use_planthydro.eq.itrue ) then
                      call InitHydrCohort(currentSite,copyc)
                   endif
+		  call copy_cohort(currentCohort, copyc)
 
                   newarea = currentCohort%c_area - cc_loss
                   copyc%n = currentCohort%n*newarea/currentCohort%c_area 
@@ -851,11 +851,11 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr)
 
                      allocate(copyc)
                      call InitPRTCohort(copyc)
-                     call copy_cohort(currentCohort, copyc) !makes an identical copy...
                      if( hlm_use_planthydro.eq.itrue ) then
                         call InitHydrCohort(CurrentSite,copyc)
                      endif
-
+                     call copy_cohort(currentCohort, copyc) !makes an identical copy...
+
                      newarea = currentCohort%c_area - cc_gain !new area of existing cohort
 
                      call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
@@ -1422,11 +1422,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
                    ! is obscured by snow.
 
                    layer_top_hite = currentCohort%hite - &
-                         ( dble(iv-1.0)/currentCohort%NV * currentCohort%hite *  &
+                         ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite *  &
                          EDPftvarcon_inst%crown(currentCohort%pft) )
 
                    layer_bottom_hite = currentCohort%hite - &
-                         ( dble(iv)/currentCohort%NV * currentCohort%hite * &
+                         ( real(iv,r8)/currentCohort%NV * currentCohort%hite * &
                          EDPftvarcon_inst%crown(currentCohort%pft) )
 
                    fraction_exposed = 1.0_r8
@@ -1449,7 +1449,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
 
                    if(iv==currentCohort%NV) then
                       remainder = (currentCohort%treelai + currentCohort%treesai) - &
-                            (dinc_ed*dble(currentCohort%nv-1.0_r8))
+                            (dinc_ed*real(currentCohort%nv-1,r8))
                       if(remainder > dinc_ed )then
                          write(fates_log(), *)'ED: issue with remainder', &
                                currentCohort%treelai,currentCohort%treesai,dinc_ed, & 

biogeochem/EDCohortDynamicsMod.F90

@@ -372,6 +372,7 @@ subroutine nan_cohort(cc_p)
     currentCohort%NV                 = fates_unset_int  ! Number of leaf layers: -
     currentCohort%status_coh         = fates_unset_int  ! growth status of plant  (2 = leaves on , 1 = leaves off)
     currentCohort%size_class         = fates_unset_int  ! size class index
+    currentCohort%size_class_lasttimestep = fates_unset_int  ! size class index
     currentCohort%size_by_pft_class  = fates_unset_int  ! size by pft classification index
 
     currentCohort%n                  = nan ! number of individuals in cohort per 'area' (10000m2 default)     
@@ -476,6 +477,8 @@ subroutine zero_cohort(cc_p)
     currentcohort%ts_net_uptake(:)   = 0._r8
     currentcohort%seed_prod          = 0._r8
     currentcohort%fraction_crown_burned = 0._r8 
+    currentCohort%size_class            = 1
+    currentCohort%size_class_lasttimestep = 0
     currentcohort%npp_acc_hold       = 0._r8 
     currentcohort%gpp_acc_hold       = 0._r8  
     currentcohort%dmort              = 0._r8 
@@ -489,7 +492,6 @@ subroutine zero_cohort(cc_p)
     currentcohort%prom_weight        = 0._r8
     currentcohort%crownfire_mort     = 0._r8
     currentcohort%cambial_mort       = 0._r8
-
 
   end subroutine zero_cohort
 
@@ -728,6 +730,9 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
      real(r8) :: diff
      real(r8) :: dynamic_fusion_tolerance
 
+     integer  :: largersc, smallersc, sc_i        ! indices for tracking the growth flux caused by fusion
+     real(r8) :: larger_n, smaller_n
+
      logical, parameter :: fuse_debug = .false.   ! This debug is over-verbose
                                                  ! and gets its own flag
 
@@ -847,6 +852,42 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                 currentCohort%canopy_layer_yesterday  = (currentCohort%n*currentCohort%canopy_layer_yesterday  + &
                                       nextc%n*nextc%canopy_layer_yesterday)/newn
 
+                                ! keep track of the size class bins so that we can monitor growth fluxes
+                                ! compare the values.  if they are the same, then nothing needs to be done. if not, track the diagnostic flux
+                                if (currentCohort%size_class_lasttimestep .ne. nextc%size_class_lasttimestep ) then
+                                   !
+                                   ! keep track of which was which, irresespective of which cohort they were in
+                                   if (currentCohort%size_class_lasttimestep .gt. nextc%size_class_lasttimestep) then
+                                      largersc = currentCohort%size_class_lasttimestep
+                                      smallersc = nextc%size_class_lasttimestep
+                                      larger_n = currentCohort%n
+                                      smaller_n = nextc%n
+                                   else
+                                      largersc = nextc%size_class_lasttimestep
+                                      smallersc = currentCohort%size_class_lasttimestep
+                                      larger_n = nextc%n
+                                      smaller_n = currentCohort%n
+                                   endif
+                                   !
+                                   ! it is possible that fusion has caused cohorts separated by at least two size bin deltas to join.  
+                                   ! so slightly complicated to keep track of because the resulting cohort could be in one of the old bins or in between
+                                   ! structure as a loop to handle the general case
+                                   !
+                                   ! first the positive growth case
+                                   do sc_i = smallersc + 1, currentCohort%size_class
+                                      currentSite%growthflux_fusion(sc_i, currentCohort%pft) = &
+                                           currentSite%growthflux_fusion(sc_i, currentCohort%pft) + smaller_n
+                                   end do
+                                   !
+                                   ! next the negative growth case
+                                   do sc_i = currentCohort%size_class + 1, largersc
+                                      currentSite%growthflux_fusion(sc_i, currentCohort%pft) = &
+                                           currentSite%growthflux_fusion(sc_i, currentCohort%pft) - larger_n
+                                   end do
+                                   ! now that we've tracked the change flux.  reset the memory of the prior timestep
+                                   currentCohort%size_class_lasttimestep = currentCohort%size_class
+                                endif
+
                                 ! Flux and biophysics variables have not been calculated for recruits we just default to 
                                 ! their initization values, which should be the same for eahc
 
@@ -1215,6 +1256,7 @@ subroutine copy_cohort( currentCohort,copyc )
     n%excl_weight     = o%excl_weight               
     n%prom_weight     = o%prom_weight               
     n%size_class      = o%size_class
+    n%size_class_lasttimestep      = o%size_class_lasttimestep
     n%size_by_pft_class = o%size_by_pft_class
 
     ! This transfers the PRT objects over.
@@ -1292,6 +1334,7 @@ subroutine copy_cohort( currentCohort,copyc )
 
     ! indices for binning
     n%size_class      = o%size_class
+    n%size_class_lasttimestep      = o%size_class_lasttimestep
     n%size_by_pft_class   = o%size_by_pft_class
 
     !Pointers

biogeochem/EDPatchDynamicsMod.F90

@@ -350,12 +350,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
@@ -377,6 +382,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. &
@@ -543,7 +550,7 @@ subroutine spawn_patches( currentSite, bc_in)
                      nc%n * currentCohort%crownfire_mort / hlm_freq_day
                 currentSite%fmort_carbonflux(levcan) = currentSite%fmort_carbonflux(levcan) + &
                      (nc%n * currentCohort%fire_mort) * &
-                     currentCohort%b_total() * g_per_kg * days_per_sec * ha_per_m2
+                     total_c * g_per_kg * days_per_sec * ha_per_m2
 
                 ! loss of individual from fire in new patch.
                 nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) 
@@ -690,10 +697,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
 
@@ -706,6 +709,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. 
@@ -731,9 +740,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
 
   ! ============================================================================
@@ -1153,14 +1164,20 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
              !not right to recalcualte dmort here.
              canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day * fates_mortality_disturbance_fraction)
 
-
-
              canopy_mortality_woody_litter(p)= canopy_mortality_woody_litter(p) + &
                   canopy_dead*(struct_c + sapw_c)
              canopy_mortality_leaf_litter(p) = canopy_mortality_leaf_litter(p) + &
                   canopy_dead*leaf_c
+
+             ! Some plants upon death will transfer storage carbon to seed production
+             ! Storage carbon that is not transferred to seeds goes to root litter flux
+
              canopy_mortality_root_litter(p) = canopy_mortality_root_litter(p) + &
-                  canopy_dead*(fnrt_c + store_c)
+                  canopy_dead*(fnrt_c + store_c*(1.0_r8-EDPftvarcon_inst%allom_frbstor_repro(p)) )
+
+             currentSite%seed_bank(p) = currentSite%seed_bank(p) + &
+                   canopy_dead * store_c * EDPftvarcon_inst%allom_frbstor_repro(p)/AREA
+
 
              if( hlm_use_planthydro == itrue ) then
                 call AccumulateMortalityWaterStorage(currentSite,currentCohort, canopy_dead)
@@ -1324,6 +1341,8 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l
     new_patch%frac_burnt         = 0._r8  
     new_patch%total_tree_area    = 0.0_r8  
     new_patch%NCL_p              = 1
+
+
 
   end subroutine create_patch
 
@@ -1448,6 +1467,12 @@ subroutine zero_patch(cp_p)
     currentPatch%c_stomata                  = 0.0_r8 ! This is calculated immediately before use
     currentPatch%c_lblayer                  = 0.0_r8
 
+    currentPatch%solar_zenith_flag          = .false.
+    currentPatch%solar_zenith_angle         = nan
+
+    currentPatch%gnd_alb_dir(:)             = nan
+    currentPatch%gnd_alb_dif(:)             = nan
+
   end subroutine zero_patch
 
   ! ============================================================================

biogeochem/EDPhysiologyMod.F90

@@ -172,7 +172,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in )
     ! update fragmenting pool fluxes
     call cwd_input( currentSite, currentPatch)
     call cwd_out( currentSite, currentPatch, bc_in)
-
+  
     do p = 1,numpft
        currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + &
             (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - &
@@ -589,7 +589,7 @@ subroutine phenology( currentSite, bc_in )
     if ((t >= currentSite%dleafondate - 30.and.t <= currentSite%dleafondate + 30).or.(t > 360 - 15.and. &
          currentSite%dleafondate < 15))then ! are we in the window?
        ! TODO: CHANGE THIS MATH, MOVE THE DENOMENATOR OUTSIDE OF THE SUM (rgk 01-2017)
-       if (sum(currentSite%water_memory(1:numWaterMem)/dble(numWaterMem)) &
+       if (sum(currentSite%water_memory(1:numWaterMem)/real(numWaterMem,r8)) &
             >= ED_val_phen_drought_threshold.and.currentSite%dstatus == 1.and.t >= 10)then 
           ! leave some minimum time between leaf off and leaf on to prevent 'flickering'.  
           if (timesincedleafoff > ED_val_phen_doff_time)then  
@@ -792,6 +792,7 @@ subroutine seeds_in( currentSite, cp_pnt )
     type(ed_cohort_type), pointer :: currentCohort
     integer :: p
     logical :: pft_present(maxpft)
+    real(r8) :: store_c_to_repro   ! carbon sent from storage to reproduction upon death [kg/plant]
     real(r8) :: npfts_present
     !----------------------------------------------------------------------
 
@@ -822,10 +823,18 @@ subroutine seeds_in( currentSite, cp_pnt )
        currentPatch => cp_pnt
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort))
+
+          ! a certain fraction of bstore goes to clonal reproduction when plants die
+          store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * &
+                EDPftvarcon_inst%allom_frbstor_repro(currentCohort%pft)
+
           do p = 1, numpft
              if (pft_present(p)) then
-                currentPatch%seeds_in(p) = currentPatch%seeds_in(p) +  currentCohort%seed_prod * currentCohort%n / &
-                     (currentPatch%area * npfts_present)
+
+                  currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + &
+                        (currentCohort%seed_prod * currentCohort%n - &
+                        currentCohort%dndt*store_c_to_repro) &
+                        /(currentPatch%area * npfts_present)		  
              endif
           end do
           currentCohort => currentCohort%shorter
@@ -836,8 +845,15 @@ subroutine seeds_in( currentSite, cp_pnt )
     currentCohort => currentPatch%tallest
     do while (associated(currentCohort))
        p = currentCohort%pft
+
+       ! a certain fraction of bstore goes to clonal reproduction when plants die
+       store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * &
+             EDPftvarcon_inst%allom_frbstor_repro(p)
+
        currentPatch%seeds_in(p) = currentPatch%seeds_in(p) +  &
-             currentCohort%seed_prod * currentCohort%n/currentPatch%area
+           (currentCohort%seed_prod * currentCohort%n - &
+	   currentCohort%dndt*store_c_to_repro)/currentPatch%area
+
        currentCohort => currentCohort%shorter
     enddo !cohort loop
 
@@ -913,13 +929,18 @@ subroutine seed_germination( currentSite, currentPatch )
 
     do p = 1,numpft
        currentPatch%seed_germination(p) =  min(currentSite%seed_bank(p) * &
-             EDPftvarcon_inst%germination_timescale(p),max_germination)
+             EDPftvarcon_inst%germination_timescale(p),max_germination)     
+       !set the germination only under the growing season...c.xu
+       if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == 1)then 
+             currentPatch%seed_germination(p) = 0.0_r8
+       endif
+       if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == 1)then
+             currentPatch%seed_germination(p) = 0.0_r8
+       endif
     enddo
 
   end subroutine seed_germination
 
-  ! ============================================================================
-
   subroutine recruitment( currentSite, currentPatch, bc_in )
     !
     ! !DESCRIPTION:
@@ -1134,9 +1155,8 @@ subroutine CWD_Input( currentSite, currentPatch)
                 ! the litter flux has already been counted since it captured
                 ! the losses of live trees and those flagged for death
 
-
           currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + &
-               (fnrt_c + store_c ) * dead_n
+               (fnrt_c + store_c*(1._r8-EDPftvarcon_inst%allom_frbstor_repro(pft)) ) * dead_n
 
           ! Update diagnostics that track resource management
           currentSite%resources_management%delta_litter_stock  = &

biogeochem/FatesAllometryMod.F90

@@ -2156,7 +2156,7 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h )
      ! of the diameter increment
      counter = 0
      step_frac = step_frac0
-     do while( (bdead-bt_dead) > calloc_abs_error )
+     do while( (bdead-bt_dead) > calloc_abs_error .and. dbt_dead_dd>0.0_r8)
 
         ! vulnerable to div0
         dd    = step_frac*(bdead-bt_dead)/dbt_dead_dd

biogeophys/EDBtranMod.F90

@@ -25,6 +25,7 @@ module EDBtranMod
 
    public :: btran_ed
    public :: get_active_suction_layers
+   public :: check_layer_water
 
 contains 
 
@@ -197,7 +198,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out)
                             cpatch%rootr_ft(ft,j) * pftgs(ft)/sum_pftgs
                     else
                        bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + &
-                            cpatch%rootr_ft(ft,j) * 1._r8/dble(numpft)
+                            cpatch%rootr_ft(ft,j) * 1._r8/real(numpft,r8)
                     end if
                  enddo
               enddo