nlevleaf too low is causing TotalBalanceCheck fails

[JessicaNeedham]

I’m getting mysterious failures in global fixed biogeography no comp mode runs.

I’m on @glemieux 's E3SM branch https://github.com/glemieux/E3SM/commits/lnd/fates-refactor commit b8be65d, with FATES main branch, commit 905bf1141196 from two weeks ago and a 4x5 grid. I am attempting a two way parameter sensitivity analysis of growth and maintenance respiration parameters (grperc and fates_maintresp_leaf_atkin2017_baserate).

Twenty two out of thirty six runs failed, with lower combinations of growth and maintenance respiration parameters causing fails (details here). Growth respiration ranged from 0.05 to 0.35 and the maintenance respiration base rate, which varies by pft, was multiplied by scalars from 0.7 to 1.4. Possibly these parameter combinations are unrealistic?

The runs fail between 1 and 20 years due to TotalBalanceCheck (call index 1) in EDMainMod and all show that gpp_acc and aresp_acc are both NaN. e.g.

105: mass balance error detected
105: element type (see PRTGenericMod.F90): 1
105: error fraction relative to biomass stock: NaN
105: absolut error (flux in - change): NaN
105: call index: 1
105: Element index (PARTEH global): 1
105: net: NaN
105: dstock: 28.611680557110958
105: seed_in: 0.0000000000000000
105: net_root_uptake: 0.0000000000000000
105: gpp_acc: NaN
105: flux_generic_in: 0.0000000000000000
105: wood_product: 0.0000000000000000
105: error from patch resizing: 0.0000000000000000
105: burn_flux_to_atm: 0.0000000000000000
105: seed_out: 0.0000000000000000
105: flux_generic_out: 0.0000000000000000
105: frag_out: 3.3555525863802003
105: aresp_acc: NaN
105: error=net_flux-dstock: NaN
105: biomass 13831.452550710850
105: litter 1875.3429779103751
105: seeds 2397.4877834851391
105: total stock 18104.283312106363 r
105: previous total 18075.671631549252
105: lat lon -26.000000000000000 315.00000000000000

Extra calls to TotalBalanceCheck show that the fail occurs after the call to ed_integrate_state_variables which is where site_cmass%aresp_acc gets updated. During the call to LeafLayerMaintenanceRespiration_Atkin_etal_2017, veg_tempk is NaN which results in dark respiration being a NaN. Vegetation temperature is updated in CanopyFluxesMod. As far as I can tell from print statements, it looks like it goes to NaN at some point during the stability iterations - around here:

https://github.com/E3SM-Project/E3SM/blob/39b5b368c0560052deae22b5e16c6acb6910d968/components/elm/src/biogeophys/CanopyFluxesMod.F90#L1050C11-L1050C11

I'm continuing to investigate.

[JessicaNeedham]

I ran with DEBUG=TRUE which showed Subscript 3 of the array FTWEIGHT has value 31 which is greater than the upper bound of 30 at this line in EDSurfaceAlbedoMod. This was causing an NaN in weighted_dif_up which was propagated through to leaf temperature and was causing the NaN in gpp_acc and aresp_acc that was triggering the TotalBalanceCheck fail. Increasing nlevleaf fixed the problem. Although I’m not sure why it failed there, or what the connection was with the respiration parameters… lower respiration causing less trimming and more leaf layers filling up… ?

My parameter file has fates_vai_top_bin_width = 0.1 and fates_vai_width_increase_factor = 1.1. I wonder if nlevleaf should also be a parameter, or if it should be calculated from the vai parameters?

[JessicaNeedham]

We could make nlevleaf dependent on vai_top_bin_width and vai_width_increase_factor. Something like:

 integer :: l                                     ! leaf layer index
 integer  :: maxnl                           ! max number of leaf layers                                                                                                 
 real(r8), allocatable :: edges(:)  ! leaf layer edges                                                                                                       
 integer :: target_cumulative_lai ! used to determine how many leaf layers are needed                                                                       

 allocate(edges(maxnl))                                                                                                                                     
 target_cumulative_lai = 30                                                                                                                                 
 maxnl = 100 

do l = 1, maxnl                                                                                                                                         
     edges(l) = ED_val_vai_top_bin_width * ED_val_vai_width_increase_factor ** l                                                                           
     if(sum(edges(1:l)) .ge. target_cumulative_lai) then                                                                                                   
         this%nlevleaf = l                                                                                                                                  
         exit                                                                                                                                               
     end if                                                                                                                                                
end do

where target_cumulative_lai is 30 which is the current default cumulative lai with nlevleaf = 30 and bins of equal width of 1 (vai_top_bin_width and vai_width_increase_factor = 1) which has always seemed to work.

The problem is where to put this. I initially thought it should go in FatesParameterDerivedMod. But in EDParamsMod there are two variables that depend on nlevleaf (dinc_vai and dlower_vai) which is going to cause problems. Anyone have suggestions or preferences @rgknox @glemieux @adrifoster @ckoven @rosiealice ?

An alternative would be to have some warning message that the user needs to adjust nleveaf if the combination of their vai parameters is going to max out the default 30 leaf layers.

[rgknox]

Parameter derived should come right afterwards, we could move dinc_vai and dlower_vai to the same routine (as long as there isn't also dependency there...) I like your solution. The work that tries to liberate us from tracking nlevleaf all-together is not straightforward, so its good to have this more direct solution right now.

[JessicaNeedham]

Thanks @rgknox . I think moving dinc_vai and dlower_vai could work. But another place I’m not sure of is FatesCohortMod and FatesPatchMod where a bunch of vegetation structures need to know nlevleaf e.g.

https://github.com/NGEET/fates/blob/73f31abe06048d2c3bda34231ceacbc8059264ba/biogeochem/FatesCohortMod.F90#L166

and

https://github.com/NGEET/fates/blob/73f31abe06048d2c3bda34231ceacbc8059264ba/biogeochem/FatesPatchMod.F90#L99

In my initial attempts to get this working (here) they won’t accept param_derived%nlevleaf. If there is an easy fix great, but I don’t think this is worth a big refactor.

We could also just increase the hard coded value of nlevleaf to 36 which is what’s needed to prevent crashes with the vai params set to their SP-mode-calibrated values of 0.1 and 1.1.