Possible NaN is introduced

[weiyuan-jiang]

GEOSldas/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90

Line 2018 in f8e314e

cat_diagS_ensstd(ii) = cat_diagS_sqrt( cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii)) )

Because of roundoff error, the result from the subtraction can become very small negative number when the two should be the same and the result should be zero. Then the function sqrt would produce NaN ( if in debug mode, it crashes). I am not sure if that affects @saraqzhang . I insert some codes :

         cat_diagS_tmp1 = cat_diagS_ensstd(ii)/Nm1 
         cat_diagS_tmp2 = NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii))
         cat_diagS_tmp = cat_diagS_tmp1 - cat_diagS_tmp2
         do i = 1, 3
            if (cat_diagS_tmp%tp(i) < 0) then
              print*,"small-big: ", cat_diagS_tmp1%tp(i), cat_diagS_tmp2%tp(i)          
            endif
         enddo

An the print out:

small-big: 659.9205 659.9206
small-big: 713.0347 713.0348
small-big: 679.2969 679.2969

@gmao-rreichle

[gmao-rreichle]

Thanks, @weiyuan-jiang, for figuring this out!
In practice, this probably has no impact on science results. When we're running an ensemble, we're running with the usual set of perturbations, and when we're running with such perturbations, the soil moisture, Tsurf, and top-layer soil temperature (tp1) all have non-negligible standard deviations within the first few time steps.
The deeper-layer soil temperatures are also included in cat_diagS and may have truly zero ens stdv for some time, until the perturbations propagate deeper down. But these deeper-layer temps are not assigned to export variables. Consequently, NaNs that may be generated by this roundoff error for deeper layer temps don't make it into the output.

It's of course a problem that the code crashes when run in debug mode because this prevents us from debugging.
The obvious fix would be to apply a max operator to the difference:

max( 0, cat_diagS_ensstd(ii)/Nm1 - NdivNm1*(cat_diagS_ensavg(ii)*cat_diagS_ensavg(ii)) )

before taking the square root. The annoying bit is that I'm not sure we have such a max operator defined for the cat_diagS type. But I suppose it wouldn't be difficult to write one.
I'm cautiously optimistic that such a fix might even be 0-diff because we may not have the exports from these calcs in the HISTORY of any of the test cases.
Either way, we'll have to fix this as soon as we can find the time.

[gmao-rreichle]

addressed with #679