At line 394 in limthd_dh.F90, an iterative procedure is used to compute (among others things) the salinity of the new ice created at the bottom of the existing layer. A zfracs is computed and used as :

s_i_new (salinity of new ice) = zfracs * sss_m (salinity of ocean)

In some cases (that seem to be extremely rare but, to be honest, not fully understood), the procedure gives zfracs > 1. The new ice salinity is then higher than ocean salinity, things go wrong and the model crashes.

Following a discussion with Martin Vancoppenolle, we've solved the problem by adding an upper bound to zfracs :

418 zds = zfracs * sss_m(ii,ij) - s_i_new(ji) 419 s_i_new(ji) = zfracs * sss_m(ii,ij) 420 ENDIF ! fc_bo_i

becomes :

418 zds = zfracs * sss_m(ii,ij) - s_i_new(ji) ----- zfracs = MIN( 0.5 , zfracs ) 419 s_i_new(ji) = zfracs * sss_m(ii,ij) 420 ENDIF ! fc_bo_i

PS : there was already a limit on the new ice salinity a few lines below :

427 ! New ice salinity must not exceed 15 psu 428 s_i_new(ji) = MIN( s_i_new(ji), s_i_max )

but :

• s_i_max is by default at 20 psu, not 15 psu (the comment should be corrected to avoid confusion)
• this limit appears to be insufficient to avoid problems
• this limit may not be necessary anymore, if the correction above is included ...

We have been seeing crashes in NEMO+LIM3 for a while and suspect uninitialised allocatable arrays in LIM3. Here is the result of a test that confirms this for the variables vt_i and vt_s:

Both variables are allocated in ice_alloc()

ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     &
   &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
   &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     &
   &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) )

If we add these lines

vt_i = HUGE(vt_i)
vt_s = HUGE(vt_s)

and run the code with floating point exceptions enabled, it gets stopped in lim_sbc_init at this point:

snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )

and using a debugger we find that both vt_i and vt_s contain the HUGE value set after allocation.

We do suspect at least one other variable (e_i), but see it as likely that there are more. In fact, with the old LIM3 version as present in NEMO 3.3.1 we have added cpp key controlled initialisation of all NEMO/LIM variables to either zero or HUGE. The model runs stable with zero-based initialisation but crashes reproducible with initialisation to HUGE.

Variable a_i is not defined in sbcice_if.F90 when key_coupled is set together with key_lim3. I've seen that in the dev_MERGE_2013 branch this was solved by using

#if defined key_coupled && defined key_lim2

but I wonder what happens for LIM3. The a_i array is defined in the "ice" module of LIM3, is it appropriate to use this?