mpaiao
commented
2 years ago 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 ifIn 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
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)forlog(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:
Slightly different version
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.