Divide by zero error in phys/module_sf_mynn.F sub zolri()

[bbrashers]

In subroutine zolri() in phys/module_sf_mynn.F, the denominator (fx2-fx1) is not protected from becoming zero. I have verified it's a divide-by-zero floating point exception in very stable conditions with shallow (< 4m) eta levels near the surface.

My suggested fix is below. The value of small is open for debate.

  REAL function zolri(ri,za,z0,zt,zol1)

  ! This iterative algorithm was taken from the revised surface layer
  ! scheme in WRF-ARW, written by Pedro Jimenez and Jimy Dudhia and
  ! summarized in Jimenez et al. (2012, MWR). This function was adapted
  ! to input the thermal roughness length, zt, (as well as z0) because
  ! zt is necessary input for the Dyer-Hicks functions used in MYNN.

  IMPLICIT NONE
  REAL, INTENT(IN) :: ri,za,z0,zt,zol1
  REAL :: x1,x2,fx1,fx2,fxdiff
  REAL, parameter :: small = 1.e-6
  INTEGER :: n

  if (ri.lt.0.)then
     x1=zol1 - 0.02  !-5.
     x2=0.
  else
     x1=0.
     x2=zol1 + 0.02 !5.
  endif

  n=0
  fx1=zolri2(x1,ri,za,z0,zt)
  fx2=zolri2(x2,ri,za,z0,zt)
  if (abs(fx2-fx1) < small) then
     fxdiff = small
  else
     fxdiff = fx2 - fx1
  endif
  Do While (abs(x1 - x2) > 0.01 .and. n < 5)
    if(abs(fx2).lt.abs(fx1))then
      x1=x1-fx1/fxdiff*(x2-x1)
      fx1=zolri2(x1,ri,za,z0,zt)
      zolri=x1
    else
      x2=x2-fx2/fxdiff*(x2-x1)
      fx2=zolri2(x2,ri,za,z0,zt)
      zolri=x2
    endif
    n=n+1
  enddo

  if (n==5 .and. abs(x1 - x2) >= 0.01) then
     !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri
     !Tests results: fails convergence ~ 0.07 % of the time
     !set approximate values:
     if (ri.lt.0.)then
        zolri=ri*5.
     else
        zolri=ri*8.
     endif
  !else
  !   print*,"iter OK, n=",n," Ri=",ri," z/L=",zolri
  endif

  return
  end function

[bbrashers]

See https://forum.mmm.ucar.edu/phpBB3/viewtopic.php?f=40&t=9946 for more information - sorry I forgot to reference that in the original post.

[davegill]

@joeolson42 @weiwangncar @dudhia @bbrashers Bart, I have included the WRF physics team and the scheme developer in this conversation. There may be a few requests for example print out, description of your case, etc. I will step aside now.

We very much appreciate when a user tracks down a bug for us, very much appreciate it. Dave

[joeolson42]

This seems like a reasonable patch. I'll add it to my evolving module and try to get it checked in soon.

Note that there has been a lot of work on this surface layer scheme in the past year, almost entirely within CCPP - not WRF. The CCPP version has recently been updated and committed to DTC's authoritative repository and is being tested in a development version of the RRFS (FV3-based successor to the HRRR). The CCPP and WRF versions have diverged considerably. I'll try to minimize that divergence as I prepare code for WRF.

Thanks for your suggested bug fix.

-joe

On Tue, Jan 26, 2021 at 1:46 PM Dave Gill notifications@github.com wrote:

@joeolson42 https://github.com/joeolson42 @weiwangncar https://github.com/weiwangncar @dudhia https://github.com/dudhia @bbrashers https://github.com/bbrashers Bart, I have included the WRF physics team and the scheme developer in this conversation. There may be a few requests for example print out, description of your case, etc. I will step aside now.

We very much appreciate when a user tracks down a bug for us, very much appreciate it. Dave

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/wrf-model/WRF/issues/1386#issuecomment-767817199, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLRR3TUEGYSKOO4XPCJT4TS34S3FANCNFSM4WUAUDOQ .

-- Joseph Olson Meteorologist Environmental Prediction Advancement Division Global Systems Laboratory NOAA-Earth System Research Labs Boulder, Colorado

[joeolson42]

Bart,

Actually, I take it back. That suggested mod will not work (as is) because fx1 and fx2 (and therefore fxdiff) need to be updated within the do-while loop. With the suggested mod, I don't think it will converge all of the time so it is probably falling back to the approximate solution. I'll have to add the same fxdiff check inside the do-while loop as well.

I wonder if the original method Pedro coded can avoid this problem because it starts with a larger difference between x1 and x2 in the first iteration. I purposely start with a smaller difference (when z/L is near 0) in an attempt to converge faster, but that may be causing a problem here.

Note that one of the changes I made to this scheme over the past year is the addition of an alternative iterative solver. It's more of a brute-force method as opposed to this two-point secant method (I think). The new brute-force method certainly avoids this problem.

-joe

On Wed, Jan 27, 2021 at 1:19 PM Joseph Olson - NOAA Federal < joseph.b.olson@noaa.gov> wrote:

This seems like a reasonable patch. I'll add it to my evolving module and try to get it checked in soon.

Note that there has been a lot of work on this surface layer scheme in the past year, almost entirely within CCPP - not WRF. The CCPP version has recently been updated and committed to DTC's authoritative repository and is being tested in a development version of the RRFS (FV3-based successor to the HRRR). The CCPP and WRF versions have diverged considerably. I'll try to minimize that divergence as I prepare code for WRF.

Thanks for your suggested bug fix.

-joe

On Tue, Jan 26, 2021 at 1:46 PM Dave Gill notifications@github.com wrote:

@joeolson42 https://github.com/joeolson42 @weiwangncar https://github.com/weiwangncar @dudhia https://github.com/dudhia @bbrashers https://github.com/bbrashers Bart, I have included the WRF physics team and the scheme developer in this conversation. There may be a few requests for example print out, description of your case, etc. I will step aside now.

We very much appreciate when a user tracks down a bug for us, very much appreciate it. Dave

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/wrf-model/WRF/issues/1386#issuecomment-767817199, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLRR3TUEGYSKOO4XPCJT4TS34S3FANCNFSM4WUAUDOQ .

-- Joseph Olson Meteorologist Environmental Prediction Advancement Division Global Systems Laboratory NOAA-Earth System Research Labs Boulder, Colorado

-- Joseph Olson Meteorologist Environmental Prediction Advancement Division Global Systems Laboratory NOAA-Earth System Research Labs Boulder, Colorado

[bbrashers]

Hi Joe,

If there's some sample code available that I can try soon, I would be very grateful - and so would the State of Alaska. . In addition that the problem earlier in this thread about zolri(), I am getting floating point exception errors in psim_unstable().

As a reminder, this is wintertime in Fairbanks, Alaska; work to support their next round of PM2.5 State Implementation Plan revisions.

Previous work (Gaudet & Stauffer) used several 3.4m deep lowest layers, which I must use to maintain regulatory consistency. 12/4/1.33km domains, one-way nesting, initialized from ERA5.

I have debug_level = 500, ran configure -D (for debugging) before compiling. One of the 32 threads ends with:

d02 2019-12-24_12:00:06 SST_UPDATE is on d02 2019-12-24_12:00:06 in MYNNSFC [c31:32462] Process received signal [c31:32462] Signal: Floating point exception (8) [c31:32462] Signal code: Invalid floating point operation (7) [c31:32462] Failing at address: 0x3dc284c [c31:32462] [ 0] /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../lib64/libpthread.so.0(+0xf5f0)[0x7f8a670665f0] [c31:32462] [ 1] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_sf_mynn_psimunstable+0x1c)[0x3dc284c] [c31:32462] [ 2] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_sf_mynnzolri2+0xb4)[0x3dc1f14] [c31:32462] [ 3] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_sf_mynnzolri+0x1b8)[0x3dc1dd8] [c31:32462] [ 4] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_sf_mynn_sfclay1dmynn+0x83e1)[0x3db8301] [c31:32462] [ 5] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_sf_mynn_sfclaymynn+0x3c75)[0x3dafde5] [c31:32462] [ 6] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_surface_driver_surfacedriver+0x12436)[0x2e89ef6] [c31:32462] [ 7] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_first_rk_step_part1_first_rk_steppart1+0x243da)[0x1f7597a] [c31:32462] [ 8] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(solveem+0x8873)[0x15caa23] [c31:32462] [ 9] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(solveinterface+0x2587)[0x13c22f7] [c31:32462] [10] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_integrateintegrate+0x34a)[0x4e44aa] [c31:32462] [11] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_integrateintegrate+0xa5a)[0x4e4bba] [c31:32462] [12] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(module_wrf_top_wrfrun+0x27)[0x48c937] [c31:32462] [13] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(MAIN_+0x35)[0x48c4a5] [c31:32462] [14] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe(main+0x44)[0x48c444] [c31:32462] [15] /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../lib64/libc.so.6(__libc_start_main+0xf5)[0x7f8a66436505] [c31:32462] [16] /usr/local/src/wrf/WRF-4.2.2-debug/main/wrf.exe[0x48c339] [c31:32462] End of error message

This was the same compilation that gave a different error message for zolri(), per my posthttps://github.com/wrf-model/WRF/issues/1386 on Github.

I suspect the ridiculously shallow layers in extremely cold and stable conditions are causing some issues with ALOG() or ATAN().

Any help would be greatly appreciated, as ADEC wants deliverables.

Bart

From: Joseph Olson @.> Sent: Wednesday, January 27, 2021 1:15 PM To: wrf-model/WRF @.> Cc: Bart Brashers @.>; Mention @.> Subject: Re: [wrf-model/WRF] Divide by zero error in phys/module_sf_mynn.F sub zolri() (#1386)

Bart,

Actually, I take it back. That suggested mod will not work (as is) because fx1 and fx2 (and therefore fxdiff) need to be updated within the do-while loop. With the suggested mod, I don't think it will converge all of the time so it is probably falling back to the approximate solution. I'll have to add the same fxdiff check inside the do-while loop as well.

I wonder if the original method Pedro coded can avoid this problem because it starts with a larger difference between x1 and x2 in the first iteration. I purposely start with a smaller difference (when z/L is near 0) in an attempt to converge faster, but that may be causing a problem here.

Note that one of the changes I made to this scheme over the past year is the addition of an alternative iterative solver. It's more of a brute-force method as opposed to this two-point secant method (I think). The new brute-force method certainly avoids this problem.

-joe

On Wed, Jan 27, 2021 at 1:19 PM Joseph Olson - NOAA Federal < @.**@.>> wrote:

This seems like a reasonable patch. I'll add it to my evolving module and try to get it checked in soon.

Note that there has been a lot of work on this surface layer scheme in the past year, almost entirely within CCPP - not WRF. The CCPP version has recently been updated and committed to DTC's authoritative repository and is being tested in a development version of the RRFS (FV3-based successor to the HRRR). The CCPP and WRF versions have diverged considerably. I'll try to minimize that divergence as I prepare code for WRF.

Thanks for your suggested bug fix.

-joe

On Tue, Jan 26, 2021 at 1:46 PM Dave Gill @.**@.>> wrote:

@joeolson42 https://github.com/joeolson42 @weiwangncar https://github.com/weiwangncar @dudhia https://github.com/dudhia @bbrashers https://github.com/bbrashers Bart, I have included the WRF physics team and the scheme developer in this conversation. There may be a few requests for example print out, description of your case, etc. I will step aside now.

We very much appreciate when a user tracks down a bug for us, very much appreciate it. Dave

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/wrf-model/WRF/issues/1386#issuecomment-767817199, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLRR3TUEGYSKOO4XPCJT4TS34S3FANCNFSM4WUAUDOQ .

-- Joseph Olson Meteorologist Environmental Prediction Advancement Division Global Systems Laboratory NOAA-Earth System Research Labs Boulder, Colorado

-- Joseph Olson Meteorologist Environmental Prediction Advancement Division Global Systems Laboratory NOAA-Earth System Research Labs Boulder, Colorado

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fwrf-model%2FWRF%2Fissues%2F1386%23issuecomment-768581623&data=04%7C01%7Cbbrashers%40ramboll.com%7C6607b75d35204f9ee06408d8c3089c56%7Cc8823c91be814f89b0246c3dd789c106%7C0%7C0%7C637473789035186123%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Qpk8ut%2B8RMIdVJQnr%2B9pFgrHxEQm2QnOVjOHyk8xuXQ%3D&reserved=0, or unsubscribehttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAGAKAC3Y4TQLYXNRISRRDO3S4B65LANCNFSM4WUAUDOQ&data=04%7C01%7Cbbrashers%40ramboll.com%7C6607b75d35204f9ee06408d8c3089c56%7Cc8823c91be814f89b0246c3dd789c106%7C0%7C0%7C637473789035186123%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=3Dddftye5UBArzLfXmPyAjTTdKM7mxzS0F3F7%2BC4Rvw%3D&reserved=0.