Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

debug compilation mode - Floating point exception in Boulac PBL scheme

LluisFB

Member
Hi,

Now it seems that in debug compilation mode, with WRF v4.5.1 (and the configuration I am using) model initializes properly. However, whilst integration, I step in a Floating point exception from the Boulac PBL scheme.

The error message (after 333 seconds of integration, 10 time steps for d01 in a 3 nested domain run) is the following one:
Code:
(...)
Caught signal 8 (Floating point exception: floating-point invalid operation)
==== backtrace (tid:3580340) ====
 0 0x0000000000012cf0 __funlockfile()  :0
 1 0x000000000a87ee78 module_bl_boulac_mp_length_bougeault_()  /ccc/work/cont003/gen6877/fitabllu/WRF/v451_git_FPS/intel/20-0-0/netcdf_4-7-4netcdff_4-5-3/dmpar/WRFdbg/phys/module_bl_boulac.f90:633
 2 0x000000000a8785db module_bl_boulac_mp_boulac1d_()  /ccc/work/cont003/gen6877/fitabllu/WRF/v451_git_FPS/intel/20-0-0/netcdf_4-7-4netcdff_4-5-3/dmpar/WRFdbg/phys/module_bl_boulac.f90:530

The corresponding line where the floating-point exception occurs is the following one in phys/module_bl_boulac.F in length_bougeault subroutine around line #632:
Code:
(...)
        do iz=kts,kte
          dld(iz)=min(dld(iz),dlg(iz))
          dls(iz)=sqrt(dlu(iz)*dld(iz))
          dlk(iz)=min(dlu(iz),dld(iz))

I suspect that the floating-point exception, might occur with negative values to perform the square root, so I modified the code introducing:
Code:
(...)
        do iz=kts,kte
          dld(iz)=min(dld(iz),dlg(iz))
          PRINT *,'  length_bougeault -- Lluis ', iz, ' dlu', dlu(iz), ' dld', dld(iz), ' prod', dlu(iz)*dld(iz)
          dls(iz)=sqrt(dlu(iz)*dld(iz))
          dlk(iz)=min(dlu(iz),dld(iz))

Running the simulation with this new line, just before the segmentation fault, I got:
Code:
   length_bougeault -- Lluis           18  dlu   715.2068      dld
  -18.28473      prod  -13077.36
Which as you can see the product within the SQRT is negative

According to comments in the code, this subroutine stands for the 'computation of the length scales for dissipation and eddy coefficients' found in boulac1D subroutine.

I am not an specialist in Boulac PBL scheme, so I would suggest that someone with a deeper knowledge of the scheme, if can purpose a solution in this case. What to do when dl[u/d] (up/down length scales) have different sign, and what does mean and what to do coherently with the parametrization. Because I suspect that it is not as simple as introduce an IF ...
 
Top