Division by zero in lim_itd_me_ridgeprep (LIM3)

[ufla]

When running r4346 of dev_MERGE_2013 I get division by zero in line 682 of LIM_SRC_3/limitd_me.F90 when astar is used:

astari = 1.0/astar

This variable should be initialised from namelist namiceitdme in subroutine lim_itd_me_init, however, the above line is reached via the call tree:

sbc --> sbc_ice_lim
      --> lim_dyn
        --> lim_rhg
          --> lim_itd_me_icestrength
            --> lim_itd_me_ridgeprep

So, in that case, lim_itd_me_init is not called before before the use of astar (and other variables?).

[ufla]

I can add that in sbc_ice_lim, lines 337--342, lim_dyn is called first (and triggers the error), whereas lim_itd_me (which in turn calls lim_itd_me_init) is called afterwards, which is too late.

[clem]

I do not see your point here. Maybe I'm wrong but to me, astar is initialized at the very beginning of the module limitd_me.F90 (via the call to lim_itd_me_init). Therefore astari, which is used later on in the same module (in lim_itd_me_ridgeprep) should also be defined right?

[ufla]

Correct, astar is initialised in lim_itd_me_init, which is called from lim_itd_me. However, the code in sbcice_lim.F90 (lines 337..341) looks like this:

337    CALL lim_dyn( kt )
       [...]
341    CALL lim_itd_me

So before lim_itd_me is called (and initialises astar via lim_itd_me_init), there is a call to lim_dyn. This leads to a call to lim_itd_me_ridgeprep, which uses astar. This means that the call to lim_dyn bypasses lim_itd_me and, hence, lim_itd_me_init.

Put in another way: There are two places that call lim_itd_me_ridgeprep, one is lim_itd_me (which makes sure all is initialised) and the other is lim_itd_me_icestrength (which bypasses initialisation).

I run into the problem when I tested the code with floating point exceptions enabled. It says division by zero and using a debugger I found the astar being zero when used in lim_itd_me_ridgeprep.

[ufla]

Just as a test, i moved the line that initialises the limitd_me module:

IF( numit == nstart  ) CALL lim_itd_me_init ! Initialization (first time-step only)

from the lim_itd_me subroutine to lim_itd_me_icestrength and that makes to code work. Of course, I have no idea about side-effects, so I'm not assuming that this is a good fix. See attached patch.

[ufla]

Patch that moves the initialisation of the limitd_me. (Not suggested as permanent fix)

[flavoni]

fix in revision 4831 of trunk:

moved call to initialization of lim_itd_me_init in iceini.F90

(this solution works also for ticket #1308)

Simona