WRF - debug: Floating point exception on p-level interpolation

LluisFB · Feb 13, 2024

Hi,

On using a WRF v4.5.1 debug compiled version, I encounter the following Floating point exception error message:

Code:

d01 2022-12-02_12:00:14  Input data is acceptable to use:
[irene2433:1139612:0:1139612] Caught signal 8 (Floating point exception: floating-p
oint divide by zero)
==== backtrace (tid:1139612) ====
 0 0x0000000000012cf0 __funlockfile()  :0
 1 0x0000000009676c7b module_diag_pld_mp_pld_()  /ccc/work/cont003/gen6877/fitabllu
/WRF/v451_git_FPS/intel/20-0-0/netcdf_4-7-4netcdff_4-5-3/dmpar/WRFdbg/phys/module_d
iag_pld.f90:235
 2 0x000000000940b5bc module_diagnostics_driver_mp_diagnostics_driver_()  /ccc/work
/cont003/gen6877/fitabllu/WRF/v451_git_FPS/intel/20-0-0/netcdf_4-7-4netcdff_4-5-3/d
mpar/WRFdbg/phys/module_diagnostics_driver.f90:867
 3 0x0000000005a0aff7 module_after_all_rk_steps_mp_after_all_rk_steps_()  /ccc/work
/cont003/gen6877/fitabllu/WRF/v451_git_FPS/intel/20-0-0/netcdf_4-7-4netcdff_4-5-3/d
mpar/WRFdbg/dyn_em/module_after_all_rk_steps.f90:171

Looking into the lines #230-#236 of phys/module_diag_pld.f90:

Code:

(...)
                     eu = qu * pu * 0.01 / ( eps + qu )
                     ed = qd * pd * 0.01 / ( eps + qd )

                     eu = max(eu, 0.001)
                     ed = max(ed, 0.001)

                     du = t_kelvin + ( s1 / ((s2 / log(eu/s3)) - 1.0) )
                     dd = t_kelvin + ( s1 / ((s2 / log(ed/s3)) - 1.0) )

I suspect that the error might come from a 0 divisor or a negative logarithm, thus I printed the values to verify it by introducing into phys/module_diag_pld.F (be aware, I modify the .F file, which is the one used when re-compiling!):

Code:

                     !  5. Dewpoint (K) - Use Bolton's approximation

                     eu = qu * pu * 0.01 / ( eps + qu ) ! water vapor press (mb)
                     ed = qd * pd * 0.01 / ( eps + qd ) ! water vapor press (mb)
                     eu = max(eu, 0.001)
                     ed = max(ed, 0.001)
                     PRINT *,'Lluis eu:', eu, ' s3:', s3, ' eu/s3:', eu/s3, ' s2:', s2
                     PRINT *,'   Lluis 2 s2 / log(eu/s3)', s2 / log(eu/s3)
                     PRINT *,'   Lluis 2 (s2 / log(eu/s3))-1.0', (s2 / log(eu/s3)) -1.0
                     du = t_kelvin + ( s1 / ((s2 / log(eu/s3)) - 1.0) )
                     dd = t_kelvin + ( s1 / ((s2 / log(ed/s3)) - 1.0) )

In doing that, one obtains:

Code:

(...)
    Lluis 2 s2 / log(eu/s3)  -73.07224
    Lluis 2 (s2 / log(eu/s3))-1.0  -74.07224
 Lluis eu:   6.112000      s3:   6.112000      eu/s3:   1.000000      s2:   17.67000

Where eu/s3: 1.000000 would gave a log(eu/s3) = 0. therfore s2 / log(eu/s3)

eu = qu * pu * 0.01 / ( eps + qu ): water vapour press at the level above (k+1), with eps = 0.622
s3 = svp1*10., where svp1 = 0.6112 : constant for saturation vapor pressure calculation (dimensionless) (from share/module_model_constants.F)

A quick, non trully 100% physic trick, could be to perturb the water vapour pressure for these cases within the error margin. In Bolton's article is cited an error of 0.1 %, to in the code this could be done:

Code:

                     IF (eu == s3) eu = eu*1.001
                     IF (ed == s3) ed = ed*1.001
                     du = t_kelvin + ( s1 / ((s2 / log(eu/s3)) - 1.0) )
                     dd = t_kelvin + ( s1 / ((s2 / log(ed/s3)) - 1.0) )

WRF - debug: Floating point exception on p-level interpolation

LluisFB

Member