Grass cohorts causing Floating-point exception at log() operations, several places

[mccabete]

UPDATE: I have encountered several places where a a log() statement of a pft-based number casues a numerical error. Below in the allometry is the first time I encountered it. I have been substituting log(value) for log(max(value, 0.001)) when I encounter them.

I am assuming what is happening is that as new grass cohorts are being made, they are being made smaller and smaller as the system reaches its asymptote, and eventually something becomes negative.

• How many decimals after 0 make sense? I want something that is bigger than ED2's numeric errors, but small enough that it won't throw stuff off. Is there a standard way people have set these bounds?

• Is there a check for cohort size that I'm missing that would prevent these numerical issues entirely? I assume that if some of of these variables are becoming negative, it means that somewhere earlier in the code the cohort shouldn't have been made in the first place. Certainly don't want my max() statements to let the model keep adding cohorts indefinitely.

I have encountered this issue:

• in dbh in the allometry code.
• in hight in the canopy structure code.
• basal area in update_derived_utils code.

Slightly different version

• lai in radiation code.

First version of this issue

Runs will successfully grow grasses for a handful of months, but will eventually fail will a floating point exception error (See below) at this line:

dbh2h = exp (b1Ht(ipft) + b2Ht(ipft) * log(mdbh) )

ED2 will run the above line with problems several times before failing on it.

 - Simulating:   09/05/2013 00:00:00 UTC
default grass allometry
Entering allom default ED-2.1
default grass allometry
Entering allom default ED-2.1
Entering allom default ED-2.1
default grass allometry
Entering allom default ED-2.1
Entering allom default ED-2.1
default grass allometry
Entering allom default ED-2.1

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x2B3198FE46D7
#1  0x2B3198FE4D1E
#2  0x2B319B7463FF
#3  0x2B319931F0CC
#4  0xF26C67 in __allometry_MOD_dbh2h at allometry.f90:110 (discriminator 9)
#5  0xF23210 in __allometry_MOD_bl2h at allometry.f90:528 (discriminator 4)
#6  0x122BFCC in __growth_balive_MOD_dbalive_dt at growth_balive.f90:379 (discriminator 12)
#7  0xF1A27D in __vegetation_dynamics_MOD_veg_dynamics_driver at vegetation_dynamics.f90:103 (discriminator 4)
#8  0x547B59 in ed_model_ at ed_model.F90:497
#9  0x42E385 in ed_driver_ at ed_driver.F90:381
#10  0x429EAD in MAIN__ at edmain.F90:290
/var/spool/sge/scc-tb1/job_scripts/1286440: line 55:   814 Floating point exception"/projectnb/dietzelab/mccabete/ED/Ed2/ED/build/ed_2.2-dbg-serdp_2020_version-5145133" ""

[mpaiao]

It is unclear to me whether the FPE is caused by log or the exponential. Maybe it's worth adding some debugging if statements to understand what is causing the crash.

For example, right after this line mdbh = min(dbh,dbh_crit(ipft)) , I would temporarily add this code and see what happens:

if (mdbh < tiny_num) then
   write(unit=*,fmt='(a)') ' Zero DBH found!'
   write(unit=*,fmt='(a,1x,i6)') 'PFT = ', ipft
   write(unit=*,fmt='(a,1x,es12.5)') 'Actual DBH = ', dbh
   write(unit=*,fmt='(a,1x,es12.5)') 'Critical DBH = ', dbh_crit(ipft)
   write(unit=*,fmt='(a,1x,es12.5)') 'b1Ht = ', b1Ht(ipft)
   write(unit=*,fmt='(a,1x,es12.5)') 'b2Ht = ', b2Ht(ipft)
   call fatal_error('Invalid DBH','dbh2h','allometry.f90')
end if

(You will need to add use consts_coms, only : tiny_num to use the pre-defined variable).

If the run continues to crash and does not print the error message, then it's likely that the exponential is the culprit, not the log. In this case, you modify the dbh2h = exp (b1Ht(ipft) + b2Ht(ipft) * log(mdbh) ) line with the following:

lnexp = b1Ht(ipft) + b2Ht(ipft) * log(mdbh)
if (lnexp >= lnexp_min .and. lnexp <= lnexp_max) then
   dbh2h = exp(lnexp)
else
   write(unit=*,fmt='(a)') ' FPE when computing height!'
   write(unit=*,fmt='(a,1x,i6)') 'PFT = ', ipft
   write(unit=*,fmt='(a,1x,es12.5)') 'Actual DBH = ', dbh
   write(unit=*,fmt='(a,1x,es12.5)') 'Critical DBH = ', dbh_crit(ipft)
   write(unit=*,fmt='(a,1x,es12.5)') 'DBH used = ', mdbh
   write(unit=*,fmt='(a,1x,es12.5)') 'b1Ht = ', b1Ht(ipft)
   write(unit=*,fmt='(a,1x,es12.5)') 'b2Ht = ', b2Ht(ipft)
   call fatal_error('Dangerous exponential','dbh2h','allometry.f90')
end if

In general, I would advise against adding log(max(value, 0.001)) unless you are sure value could be indeed zero. In this case, the crash may be indicating some other problem, because DBH should never be zero in ED2.

[mccabete]

👍 print statements.

Tentatively, I think it's the log. When I put in the hack I described, things "worked", as in that line didn't cause an error, and the run ran longer.

Looking at the monthly file pre-allometry-crash:

DBH by cohort: 1.70492851734161, 1.70492851734161 0.953645884990692, 0.350115686655045
Basal Area by cohort: 2.2829806804657 2.2829806804657 0.714272856712341 0.0962748900055885
Hight by cohort: 1.5, 1.5, 1.00226330757141, 0.5

[mccabete]

Sorry it took me so long to circle back! Yep, it's the log().

 Zero DBH found!
PFT =       1
Actual DBH =   0.00000E+00
Critical DBH =   1.70493E+00
b1Ht =   3.52000E-02
b2Ht =   6.94000E-01

[mccabete]

I put some print statements into the bl2dbh() function, that I think is generating the dbh values. Seems like I get the FPE when leaf biomass drops to zero:

PFT =       1
bleaf  =   0.00000E+00
l1DBH =   1.39617E+01
l2DBH =   5.79043E-01
 Zero DBH found!
PFT =       1
Actual DBH =   0.00000E+00
Critical DBH =   1.70493E+00
b1Ht =   3.52000E-02
b2Ht =   6.94000E-01

[mccabete]

Seems like there just needs to be a check for bleaf > tiny_num before this line where grass updates its height daily.

Not sure if that would mess up a carbon balance term somewhere. Theoretically though, it makes sense that grasses might not have the biomass to grow taller every day.

[mccabete]

Second theory: Could the grass phenology cases be being referred to incorrectly?
I ask because these lines seems to change cohort leaf biomass.

[mccabete]

If I run with the phenology scheme set to -1 I do not get the errors above. I was running phenology 1 before.

[mpaiao]

Are you using the NL%IGRASS=1? This should force grasses to be evergreen, but maybe this capability was inadvertently lost at some point. If you are using NL%IGRASS=1, try running the deciduous with NL%IGRASS=0 and see if it crashes.

[mccabete]

Yeah IGRASS = 1

[mccabete]

I'll try a phenology = 1 and a IGRASS = 0 run

[mccabete]

Yes! If I run with is_grass = 0 and deciduous phenology, I get no problems.