CTSM/FATES enters an 'infinte loop' on Dec. 31 in trying to demote cohorts to balance areas

[mvertens]

Brief summary of bug

I have been working with @mvdebolskiy to try to understand a crash that happens at the end of year 1 on Dec. 31 with the following behavior.

code repository: 
   https://github.com/mvdebolskiy/CTSM.git
code branch: 
   updt-noresm-to-5.3.084
test case: 
   SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.cl\m-FatesColdNoComp.matvey/

• The code dies on Dec. 31 in EDCanopyStructureMod.F90 at first with an array out of bounds for arealayer(nclmax+5). It turns out that arealayer does not have to be an array - so I have made it a scalar to get past this point.
• Now it dies upon failure of the condition (patch_area_counter > max_patch_iterations .and. area_not_balanced)
• This happens at numerous gridcells - and I have picked just one gridcell failure to look at more carefully -
  lat, lon = 1.33058, 27.50000 and the crash happens with pft=13. I'm using 1024 tasks - and task 134: contains the problem gridcell.
• I have restart files written at both Dec. 30 and Dec. 31.
• On dec. 30:

134: DEBUG: lat, lon, year, mon, day, tod =       1.33058       27.50000  2000    12    30  1800
134:
134: DEBUG: area not balanced loop: pft, counter =   13     0
134: DEBUG: z1 =        3
134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    1        720.9159737710        714.2857142857          6.6302594853
134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   1     6.6302594853   631.1649323585
134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    2        724.1804081137        714.2857142857          9.8946938280
134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   1     9.8946938280   634.4244105726
134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    3        726.7865900232        714.2857142857         12.5008757375
134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1     0.6813581681     0.6813581681
134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   2    11.8195175694   716.2105380271
134: DEBUG: z2  =        4
13

On dec. 31:

 134: DEBUG: lat, lon, year, mon, day, tod =       1.33058       27.50000  2000    12    31  1800
 134:
 134: DEBUG: area not balanced loop: pft, counter =   13     0
 134: DEBUG: z1 =        3
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    1       8088.8092821351        714.2857142857       7374.5235678494
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   1  7374.5235678494  7661.6637274354
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    2      12887.8364490678        714.2857142857      12173.5507347821
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1    37.8177039835    37.8177039835
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   2  5475.4951772349  5475.4951772349
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   3  6660.2378535637  7374.5235678494
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    3      15777.8579507025        714.2857142857      15063.5722364168
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1     0.6853056261     0.6853056261
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   2    37.8177039835    37.8177039835
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   3  3603.6219102943  3603.6219102943
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   4  5475.4951772349  5475.4951772349
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   5  5945.9521392780  6660.2378535637
 134: DEBUG: z2  =        4
 134:
  134:
 134: DEBUG: area not balanced loop: pft, counter =   13     1
 134: DEBUG: z1 =        4
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    1        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    2        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    3        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   1     0.0000000000   714.2857142857
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    4      15063.5722364168        714.2857142857      14349.2865221311
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1     0.6853056261     0.6853056261
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   2    37.8177039835    37.8177039835
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   3  3603.6219102943  3603.6219102943
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   4  5475.4951772349  5475.4951772349
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   5  5231.6664249923  5945.9521392780
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   6     0.0000000000     0.0000000000
 134: DEBUG: z2  =        5
 134:
 134:
 134: DEBUG: area not balanced loop: pft, counter =   13     2
 134: DEBUG: z1 =        5
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    1        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    2        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    3        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    4        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    5      14349.2865221311        714.2857142857      13635.0008078454
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1     0.6853056261     0.6853056261
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   2    37.8177039835    37.8177039835
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   3  3603.6219102943  3603.6219102943
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   4  5475.4951772349  5475.4951772349
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   5  4517.3807107066  5231.6664249923
 134: DEBUG: z2  =        6
......
and demotion continues until the following last iteration of the loop which causes an abort
.......
 134:
 134: DEBUG: area not balanced loop: pft, counter =   13    10
 134: DEBUG: z1 =       13
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    1        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    2        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    3        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    4        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    5        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    6        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    7        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    8        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =    9        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =   10        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: target_area is less than nearzero - returning
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =   11        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =   12        714.2857142857        714.2857142857          0.0000000000
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   1     0.0000000000   482.6192892933
 134: DEBUG: canopy_structure i_lyr, arealayer, currentPatch%area, target_area =   13       8635.0008078455        714.2857142857       7920.7150935598
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   1     0.6853056261     0.6853056261
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   2    37.8177039835    37.8177039835
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   3  3603.6219102943  3603.6219102943
 134: DEBUG: demotion partial: cohort#, pd_area, cohort_area   4  4278.5901736559  4992.8758879416
 134: DEBUG: demotion total  : cohort#, pd_area, cohort_area   5     0.0000000000     0.0000000000
 134: DEBUG: z2  =       14
 134:
 134:  PATCH AREA CHECK NOT CLOSING

The problem seems to be that the arealayer for pft 14 coming into this routine on Dec. 31 is huge and there are not enough demotion calls to totally resolve it. I'm still trying to spin up on FATES - so maybe I am missing how this is happening. @mvdebolskiy - do you have other data to add since we last talked?

[mvdebolskiy]

• Running ctsm5.3.084 on derecho with NGEET/fates:main:HEAD paramter file, and the pr branch Updt noresm to 5.3.084 #170 and Update to API42 plus fixes for nbp fates#43 with that parameter file passes the SMS_D_Ld366.ne30pg3_*.I2000Fates.*.clm-FatesColdNoComp test.

• Running with parameter file from NorESMhub/fates:noresm:HEAD fails the test at RUN with both setups from above.

• There is a difference in fusion thresholds on the parameter files, and reverting these on noresm paramfile to values in ngeet one does not help.

• Merging Protections on promotion/demotion for near zero sized cohorts NGEET/fates#1509 does not help.

• Full fates with noresm parameters gives WaterBalance errors when running more than one month (we do not have global tests for that long).
  Resolution does not matter.

• Tests FatesColdNoCompFixedBioGeo with noresm parameter file fail at month 3 with fpe in the SoilTemperatureMod.F90:1007 where division by a tiny z_avg happens.

• In addition, though we are runing without landuse, there are always warnings:

Warning. small nocomp patch wasnt able to find another patch to fuse with.           0           0   8.8333450427411902E-003

flooding the cesm.log for both NoComp and NoCompFixedBioGeo.Usually this happens with currentPatch%nocomp_pft_label, currentPatch%land_use_label being zeros, however, in some cases it's actual pfts.

The error seems to occur due to a small area diff in canopy area within a layer and patch area, which does not get resolved within nclmax (which is increased to 3 in the ngeet/main) and blows up at around level 5 to huge values.

Tagging @rgknox.

[mvertens]

@rosiealice @kjetilaas @JessicaNeedham @maritsandstad - it would be helpful to have consensus as to what we should put together for noresm3_0_beta08. Should we have a quick meeting before Friday?

[mvdebolskiy]

• fpe overflow in SoilTemperatureMod.F90 also happen in noresm:HEAD but at much later date. for NoCompFixBioGeo.

[mvertens]

@mvdebolskiy - I think we should not have the condition at line 1005

            z_avg=frac_sno(c)*snow_depth(c)
            if (z_avg > 0._r8) then            <=======
               rho_avg=min(800._r8,h2osno_total(c)/z_avg)
            else
               rho_avg=200._r8
            endif

But replace it with a cutoff (say cutof = 1.e-10_r8)

            z_avg=frac_sno(c)*snow_depth(c)
            cuttoff = 1.e-10_r8
            if (z_avg > cuttoff) then            <=======
               rho_avg=min(800._r8,h2osno_total(c)/z_avg)
            else
               rho_avg=200._r8
            endif

I see several other places in this routine where you could potentially have an overflow where you divide by frac_sno and check beforehand if frac_sno > 0._r8. It's interesting that this has not overflowed before.

[kjetilaas]

@rosiealice @kjetilaas @JessicaNeedham @maritsandstad - it would be helpful to have consensus as to what we should put together for noresm3_0_beta08. Should we have a quick meeting before Friday?

Thanks for all this, @mvertens and @mvdebolskiy!

I agree that it would be good to have a quick meeting before Friday. From a quick look at the calendars of the cicero group, it seems tomorrow sometime between 10 and 11:30 could work. How does that sound?

[mvertens]

@kjetilaas - that time period will work for me. I can come up to Cicero if everyone will be meeting there.

[kjetilaas]

Great! I will set up the meeting!

You are very welcome to come over, although I wouldn't expect everyone to be here in person so it will anyway be hybrid.

[mvdebolskiy]

@kjetilaas I will try 10ish.

[rgknox]

Suspicious that it is on Dec 31st. Also, @mvertens , judging by the debug messages, this looks like pretty old code, perhaps before we refactored with tag sci.1.87.1 ?
See:
NGEET/fates@sci.1.87.0_api.41.0.0...sci.1.87.1_api.41.0.0

[mvertens]

@rgknox - the debug messages are my code that I inserted. If there is better output - let me know and I can use that instead. Part of this was actually me trying to understand the logic more.

[mvdebolskiy]

Running tests with ctsm5.3.045_noresm_v16:

  SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo (Overall: PASS) details:
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo CREATE_NEWCASE
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo XML
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo SETUP
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo SHAREDLIB_BUILD time=44
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo MODEL_BUILD time=53
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo SUBMIT
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo RUN time=956
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo MEMLEAK
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoCompFixedBioGeo SHORT_TERM_ARCHIVER
  SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp (Overall: PASS) details:
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp CREATE_NEWCASE
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp XML
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp SETUP
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp SHAREDLIB_BUILD time=284
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp MODEL_BUILD time=64
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp SUBMIT
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp RUN time=1439
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp MEMLEAK
    PASS SMS_D_Ld366_P1024.ne30pg3_ne30pg3_mtn14.I2000Clm60Fates.betzy_gnu.clm-FatesColdNoComp SHORT_TERM_ARCHIVER

So, it's recent changes to the paramfile that are causing z_avg FPE so early in the runs. Though it's not quite related to canopy_structure.

[rgknox]

I was able to reproduce the failure in canopy promotion/demotion with the test: SMS_D_Ld366.ne30pg3_g17.I2000Clm60Fates.derecho_intel.clm-FatesColdNoComp and a parameter file provided by @mvdebolskiy

I ran into the same issue, and tried the same fix that you recommended @mvertens: the arealayer array doesn't really need to be an array. (I'm re-running tests again now)

But this is telling of something really odd happening. At least 5 whole canopy layers were added on this day. I would like to think that the canopy promotion/demotion scheme is robust enough to cull 8+ layers of vegetation down to an orderly 3 layers, but its awfully strange that so many layers are being presented to it.

[rgknox]

For the canopy promotion/demotion area conservation crash... its possible the algorithm is actually working as intended. It will work its way down they canopy layers, pushing area down to the next lowest layer when overfull, and pulling area up from a lower layer when in deficit. On the first pass, the top most layer should have a bias of zero, with the second layer perhaps either gaining or losing area depending on the first layers needs. But lets say there was not enough area in the second layer to accomodate the first's needs, then we will need another pass. On the second pass through, area would have been moved up from lower layers so that the second layer can finish moving up to the top layer. I think what is happening here, is that since there are 14 canopy layers, it just takes many passes to get all the balancing done. First off, since we won't be keeping any layer lower than nclmax, there is no point in making sure those layers are balanced. As long as layers 1-nclmax are balanced, we can terminate those extra layers and move on. In my test run, layers 1-10 were all balanced near perfectly, and it was layers 11-14 that had not been balanced (hopefully according to my theory above, ie the balancer just had not gotten that far down yet).

[mvdebolskiy]

@mvertens that is because snow_depth is tiny, no? Though snow depth should be truncated. Or it is because their product goes into tiny(r8) range, I am not quite sure. In any case. Both of these are calculated after canopy interception that uses elai/esai.

[mvertens]

@mvdebolskiy - reverting the SoilTemperature.F90 change to not have a cutoff for z_avg and using the update FATES branch test_fusion with the terminate_cohorts inside the interation loop stil causes the model to blow up in Model Date: 2000-01-22. I'm going to revert SoilTemperature.F90 to have the cutoff for z_avg.