New ice salinity can be higher than ocean salinity

[barthele]

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 ...

[clem]

fixed in branch dev_r4028_CNRS_LIM3 to be included in the next merge party (Nov. 2013)