dynamic root code divides by 0

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2015-10-31 05:15:52 -0600
Bugzilla Id: 2237
Bugzilla CC: andre@ucar.edu, dlawren@ucar.edu, gbisht@lbl.gov, rfisher@ucar.edu,

Runs that activate both use_dynroot and crop crash in the first time step if run in debug mode. This is due to a divide by 0 error in this line in CNRootDynMod:

                root_depth(p) = max(zi(c,2), min(hui(p)/huigrain(p)* root_dmx(ivt(p)), root_dmx(ivt(p))))

huigrain(p) can sometimes be 0; in the one case I looked at, hui was also 0. Also, in the first timestep of a cold start, huigrain is NaN.

So this needs:

(1) initialization of huigrain in cold start so that it isn't NaN: it should probably be initialized to 0

(2) appropriate handling of the case when huigrain = 0

[billsacks]

Erik Kluzek < erik@ucar.edu > - 2015-11-18 14:29:31 -0700

The fix that Beth Drewniak gave us for this was to initialize huigrain to zero everywhere, and then...

Index: components/clm/src/biogeochem/CNRootDynMod.F90
===================================================================
--- components/clm/src/biogeochem/CNRootDynMod.F90	(.../branch_tags/fun20exp_tags/fun20exp_n08_fun2_0_n09clm4_5_6_r153)	(revision 75190)
+++ components/clm/src/biogeochem/CNRootDynMod.F90	(.../branches/fun20exp)	(revision 75190)
@@ -134,7 +134,9 @@
             c = pcolumn(p)
             if (ivt(p) /= noveg) then
                 if((ivt(p)) >= npcropmin)then !skip generic crop types
-                    root_depth(p) = max(zi(c,2), min(hui(p)/huigrain(p)* root_dmx(ivt(p)), root_dmx(ivt(p))))
+                    if(huigrain(p) > 0._r8)then
+                        root_depth(p) = max(zi(c,2), min(hui(p)/huigrain(p)* root_dmx(ivt(p)), root_dmx(ivt(p))))
+                    end if
                 else
                 ! this can be changed to any depth (i.e. the maximum soil depth)
                     root_depth(p) = zi(c,nlevsoi)

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2015-11-18 14:52:50 -0700

Shouldn't there be an 'else' clause here, so that root_depth gets set explicitly if huigrain(p) == 0? Or is the intention really that it remains at whatever value it was at in the previous time step?

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2015-11-23 12:33:10 -0700

Gautam's suggestion to the 'else', approved by Beth, is:

How about adding root_depth(bounds%begp:bounds%endp) = 0._r8 on CNRootDynMod.F90#L111?

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2015-11-23 12:35:28 -0700

Oops, I forgot to finish Gautam's suggestion:

Then, you can get rid of the 'else' clause at line CNRootDynMod.F90#L130.

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2015-11-23 12:39:27 -0700

However, in response to comment 3: Based on the performance impact we have seen for CNZero, I am coming to feel that we should avoid code like:

root_depth(bounds%begp:bounds%endp) = 0._r8

instead either (a) initializing this in the loop over the filter (i.e., setting root_depth(p) = 0._r8 before the conditionals), or (b) keeping / adding explicit 'else' clauses to ensure it is set to something for any path through the conditionals.

This isn't a big deal for this particular change, but is something we should think about in general: As we have seen from CNZero, when you add up the total time we spend zeroing quantities, it amounts to a significant fraction of CLM run time.

[billsacks]

Erik Kluzek < erik@ucar.edu > - 2015-11-23 13:13:54 -0700

The first part of this was completed in clm4_5_6_r155.

[billsacks]

It looks like the fix from #38 (comment) was applied, but the suggested else clause still hasn't been put in. (The same seems true in E3SM.) So, if huigrain goes to 0, I believe root_depth remains at its previous value. The one place where I can find huigrain being set to 0 is here:

https://github.com/ESCOMP/ctsm/blob/57adf1ed1f3f5d4b4881277734c6d4b67a392c64/src/biogeochem/CNPhenologyMod.F90#L2233-L2236

I'm not sure whether we should leave things as is (because root depth may not matter in cases where huigrain has been set to 0 anyway), or do the (fairly easy) suggestion above (adding an else clause).

[billsacks]

Closing as a wontfix because dynamic roots needs some work anyway.