Floating point exceptions related to GLS_MIXING and PERFECT_RESTART

[colucix]

Following the discussion in the ​forum, I'd like to report a bug triggered by the PERFECT_RESTART option and related to vertical diffusion processes for the Generic Length Scale formulation.

The variables tke and gls, together with diffusion coefficients are initialized to their background values (see subroutine initialize_mixing), but subsequently they are read from the restart NetCDF and the values on land areas (matching _FillValue in NetCDF) are set to 0 (see subroutine nf_fread4d).

The tke and gls values enter as denominators for calculations inside subroutine gls_corstep and the 0 values generate a floating point exception.

The resulting NaN are propagated in particular (but not only) to the vertical diffusion coefficient for tracers (Akt) and this causes a model blow-up in subsequent calculations.

I suggest to set the flag SetFillValue to .FALSE. (in: subroutine def_rst and subroutine wrt_rst) for variables tke and gls, in order to preserve the background values over land points. Maybe the same should be applied to the other variables added to the restart file when PERFECT_RESTART is active (Akp, Akk and Lscale) but they don't seem to be involved in potential floating point exceptions.

Please see ​post in the ROMS Bugs forum for a bit more detailed description of the issue. Thank you, Alex

[arango]

Yes, thank you for debugging this carefully and giving us such nice details. You definitely won a lot of points in my book. Few users give us such details about code problems.

This was a tricky fix and difficult bug to find. We also need to remove the water_points attribute for these variables. I also fixed this problem in the history and averages files for the case that we want to initialize the model with these files. I also made the same corrections to the tangent and adjoint version of the history file.

I don't think that we need the same treatment for the Akp, Akk, and Lscale variables since they are initialized to zero in initialize_mixing. Their values are not used in the denominator of an expression of power factors.

I also corrected a compiler problem in routine gasdev.F. The IFORT compiler fails in this routine. Obviously, this is a bug in the IFORT compiler. To be safe I forced the call to the implicit function PACK to be real:

          CALL array_copy (REAL(PACK(v1(ng:n), mask(ng:n)),r8),        &
     &                     v1(ng:), nn, m)

Many thanks to Brian Powell for tracking this problem.