Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix several bugs in fv_regional_bc.F90 relating to uninitialized or incorrectly initialized memory. #219

Merged
merged 13 commits into from
Oct 17, 2022
Merged

Fix several bugs in fv_regional_bc.F90 relating to uninitialized or incorrectly initialized memory. #219

Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
5963599
fixes and workarounds for uninitialized memory in fv_regional_bc
SamuelTrahanNOAA Sep 18, 2022
3766d58
remove workarounds and fix remaining known bugs in ps_reg
SamuelTrahanNOAA Sep 20, 2022
12c9bd1
a few more surface pressure bug fixes; now the test case runs in debu…
SamuelTrahanNOAA Sep 20, 2022
73b7e59
bug fixes to my bug fixes
SamuelTrahanNOAA Sep 20, 2022
6dd393a
Merge remote-tracking branch 'origin/dev/emc' into bugfix/rrfs-debug-…
SamuelTrahanNOAA Sep 20, 2022
5800359
Merge remote-tracking branch 'origin/dev/emc' into bugfix/rrfs-debug-…
SamuelTrahanNOAA Sep 20, 2022
7dce0d7
Merge branch 'bugfix/rrfs-debug-mode' of ssh://github.com/SamuelTraha…
SamuelTrahanNOAA Sep 20, 2022
8c67895
workarounds and bug fixes from gnu compiler testing
SamuelTrahanNOAA Sep 21, 2022
6d04c31
remove -9999999 commented-out code
SamuelTrahanNOAA Sep 21, 2022
229dbe5
quiet the NaNs passed to Atmp%ps
SamuelTrahanNOAA Sep 21, 2022
b923e9f
simplify comments and explain snan
SamuelTrahanNOAA Sep 21, 2022
68a833f
use i-1 & j-1 for two-point averages, when available
SamuelTrahanNOAA Sep 27, 2022
9824ef3
Replace many changes with PR #220
SamuelTrahanNOAA Oct 3, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 106 additions & 23 deletions model/fv_regional_bc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,67 @@ module fv_regional_mod
logical :: data_source_fv3gfs
contains


!-----------------------------------------------------------------------
!
logical function is_not_finite(val)
!
!-----------------------------------------------------------------------
!*** This routine is equivalent to ".not. ieee_is_finite(val)"
!*** which returns .true. for infinite and Not a Number (NaN), or
!*** .false. otherwise. It's here as a workaround for this gfortran bug:
!***
!*** https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82207
!***
!*** The compiler must use IEEE-standard floating point for this to work
!-----------------------------------------------------------------------
!
use, intrinsic :: iso_c_binding, only: c_int32_t, c_int64_t
implicit none
!
!-----------------------------------------------------------------------
! Portability note: shiftr() is part of Fortran 2008, but it is widely
! supported in older compilers.
!-----------------------------------------------------------------------
!
intrinsic shiftr, transfer, iand ! <--- declare intrinsic to help older compilers
!
!-----------------------------------------------------------------------
! Use value-based argument passing instead of reference-based to avoid
! signaling a NaN on conversion to addressable storage.
!-----------------------------------------------------------------------
!
real, value :: val ! <-- bit pattern to test for infinity or NaN
!
!-----------------------------------------------------------------------
! Bit manipulation constants for testing 32-bit floating-point
! non-finite values.
!-----------------------------------------------------------------------
!
#ifdef OVERLOAD_R4
integer(c_int32_t), parameter :: check = 255 ! <-- all bits on, size of exponent (8 bits)
integer, parameter :: shift = 23 ! <-- number of mantissa bits except sign
!
!-----------------------------------------------------------------------
! Bit manipulation constants for testing 64-bit floating-point
! non-finite values.
!-----------------------------------------------------------------------
!
#else
integer(c_int64_t), parameter :: check = 2047 ! <-- all bits on, size of exponent (11 bits)
integer, parameter :: shift = 52 ! <-- number of mantissa bits except sign
#endif
!
!-----------------------------------------------------------------------
! For IEEE standard floating point numbers, non-finite values follow
! a mandatory bit pattern. They have the mantissa sign bit on, and all
! exponent bits on, except the exponent sign which can be on or off.
!-----------------------------------------------------------------------
!
is_not_finite = iand(shiftr(transfer(val,check),shift),check)==check
!
end function is_not_finite
!
!-----------------------------------------------------------------------
!
subroutine setup_regional_BC(Atm &
Expand Down Expand Up @@ -765,8 +826,8 @@ subroutine setup_regional_BC(Atm &
!*** reference pressure profile. Compute it now.
!-----------------------------------------------------------------------
!
allocate(pref(npz+1))
allocate(dum1d(npz+1))
allocate(pref(npz+1)) ; pref=real_snan
allocate(dum1d(npz+1)) ; dum1d=real_snan
!
ps1=101325.
pref(npz+1)=ps1
Expand Down Expand Up @@ -951,7 +1012,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
regional_bc_bounds%ie_north_uvs=ied
!
regional_bc_bounds%js_north_uvs=jsd
regional_bc_bounds%je_north_uvs=nrows_blend+1
regional_bc_bounds%je_north_uvs=nrows_blend
!
regional_bc_bounds%is_north_uvw=isd
regional_bc_bounds%ie_north_uvw=ied+1
Expand All @@ -968,7 +1029,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
regional_bc_bounds%is_south_uvs=isd
regional_bc_bounds%ie_south_uvs=ied
!
regional_bc_bounds%js_south_uvs=jed-nhalo_model-nrows_blend+1
regional_bc_bounds%js_south_uvs=jed-nhalo_model-nrows_blend+2
regional_bc_bounds%je_south_uvs=jed+1
!
regional_bc_bounds%is_south_uvw=isd
Expand Down Expand Up @@ -1028,7 +1089,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
regional_bc_bounds%je_west_uvs=jed-nhalo_model+1
endif
!
regional_bc_bounds%is_west_uvw=ied-nhalo_model-nrows_blend+1
regional_bc_bounds%is_west_uvw=ied-nhalo_model-nrows_blend+2
regional_bc_bounds%ie_west_uvw=ied+1
!
regional_bc_bounds%js_west_uvw=jsd
Expand Down Expand Up @@ -1307,8 +1368,8 @@ subroutine start_regional_cold_start(Atm, ak, bk, levp &
enddo
enddo
!
allocate (ak_in(1:levp+1)) !<-- Save the input vertical structure for
allocate (bk_in(1:levp+1)) ! remapping BC updates during the forecast.
allocate (ak_in(1:levp+1)) ; ak_in=real_snan !<-- Save the input vertical structure for
allocate (bk_in(1:levp+1)) ; bk_in=real_snan ! remapping BC updates during the forecast.
do k=1,levp+1
ak_in(k)=ak(k)
bk_in(k)=bk(k)
Expand Down Expand Up @@ -1402,9 +1463,9 @@ subroutine start_regional_restart(Atm &
,isd, ied, jsd, jed &
,Atm%npx, Atm%npy )
!
allocate (wk2(levp+1,2))
allocate (ak_in(levp+1)) !<-- Save the input vertical structure for
allocate (bk_in(levp+1)) ! remapping BC updates during the forecast.
allocate (wk2(levp+1,2)) ; wk2=real_snan
allocate (ak_in(levp+1)) ; ak_in=real_snan !<-- Save the input vertical structure for
allocate (bk_in(levp+1)) ; bk_in=real_snan ! remapping BC updates during the forecast.
if (Atm%flagstruct%hrrrv3_ic) then
if (open_file(Grid_input, 'INPUT/hrrr_ctrl.nc', "read", pelist=pes)) then
call read_data(Grid_input,'vcoord',wk2)
Expand Down Expand Up @@ -1908,7 +1969,8 @@ subroutine regional_bc_data(Atm,bc_hour &
!*** the integration levels.
!-----------------------------------------------------------------------
!
allocate(ps_reg(is_input:ie_input,js_input:je_input)) ; ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed
allocate(ps_reg(is_input:ie_input,js_input:je_input)) !<-- Sfc pressure in domain's boundary region derived from BC files
ps_reg=real_snan !<-- detect access of uninitialized pressures
!
!-----------------------------------------------------------------------
!*** We have the boundary variables from the BC file on the levels
Expand Down Expand Up @@ -3545,6 +3607,12 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
je=js_bc+nhalo_data+nrows_blend-1
endif
!
! Ensure uninitialized memory isn't used
pn0 = real_snan
pn1 = real_snan
gz_fv = real_snan
gz = real_snan
pn = real_snan
allocate(pe0(is:ie,km+1)) ; pe0=real_snan
allocate(qn1(is:ie,npz)) ; qn1=real_snan
allocate(dp2(is:ie,npz)) ; dp2=real_snan
Expand Down Expand Up @@ -3575,13 +3643,14 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
pn(k) = 2.*pn(km+1) - pn(l)
enddo

pst = real_snan
do k=km+k2-1, 2, -1
if( phis_reg(i,j).le.gz(k) .and. phis_reg(i,j).ge.gz(k+1) ) then
pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-phis_reg(i,j))/(gz(k)-gz(k+1))
go to 123
exit
endif
enddo
123 ps(i,j) = exp(pst)
ps(i,j) = exp(pst)

enddo ! i-loop

Expand All @@ -3594,10 +3663,10 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
!*** the Atm object.
!---------------------------------------------------------------------------------
!
is=lbound(Atm%ps,1)
ie=ubound(Atm%ps,1)
js=lbound(Atm%ps,2)
je=ubound(Atm%ps,2)
is=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),is))
ie=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),ie))
js=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),js))
je=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),je))
Comment on lines +3666 to +3669
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I doubt that the modification here is really necessary.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Without these changes, the model uses data outside the boundary regions that are overwritten with real_snan in the following loop. This error was how I found the ps_reg=-9999999 and later bugs. (Maybe this is caused by a region of the code you don't typically use?)

However, even if we assume that the changes are not necessary, it is then a harmless bounds check on copying to Atm%ps. As with initializing data to real_snan, bounds checks are reasonable safeguard.

!
do j=js,je
do i=is,ie
Expand Down Expand Up @@ -4864,6 +4933,9 @@ subroutine bc_time_interpolation(array &
do j=j1_blend,j2_blend
factor_dist=exp(-(blend_exp1+blend_exp2*(j-j_bc-1)*rdenom)) !<-- Exponential falloff of blending weights.
do i=i1_blend,i2_blend
if(is_not_finite(array(i,j,k))) then
cycle ! Outside boundary
endif
Comment on lines +4936 to +4938
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don’t run the blending boundary condition so I cannot test it. This walkaround looks fine to me although I do think that revising the indices of the blending domain is the ultimate fix.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're probably right about the indices, but I don't trust my knowledge of the code. If someone else knows the blending code, hopefully they'll review and suggest a better solution.

blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated
+(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
!
Expand All @@ -4883,6 +4955,9 @@ subroutine bc_time_interpolation(array &
do j=j1_blend,j2_blend
factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom)) !<-- Exponential falloff of blending weights.
do i=i1_blend,i2_blend
if(is_not_finite(array(i,j,k))) then
cycle ! Outside boundary
endif
blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated
+(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
Expand All @@ -4900,6 +4975,9 @@ subroutine bc_time_interpolation(array &
do k=1,ubnd_z
do j=j1_blend,j2_blend
do i=i1_blend,i2_blend
if(is_not_finite(array(i,j,k))) then
cycle ! Outside boundary
endif
!
blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated
+(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
Expand All @@ -4921,6 +4999,9 @@ subroutine bc_time_interpolation(array &
do k=1,ubnd_z
do j=j1_blend,j2_blend
do i=i1_blend,i2_blend
if(is_not_finite(array(i,j,k))) then
cycle ! Outside boundary
endif
!
blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated
+(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
Expand Down Expand Up @@ -5727,7 +5808,7 @@ subroutine dump_field_3d (domain, name, field, isd, ied, jsd, jed, nlev, stag)
nyg = npy + 2*halo + jext
nz = size(field,dim=3)

allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) )
allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) ) ; glob_field=real_snan

isection_s = is
isection_e = ie
Expand Down Expand Up @@ -5846,7 +5927,7 @@ subroutine dump_field_2d (domain, name, field, isd, ied, jsd, jed, stag)
nxg = npx + 2*halo + iext
nyg = npy + 2*halo + jext

allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) )
allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) ) ; glob_field=real_snan

isection_s = is
isection_e = ie
Expand Down Expand Up @@ -6154,6 +6235,7 @@ subroutine write_full_fields(Atm)
nz=size(fields_core(nv)%ptr,3)
!
allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
global_field=real_snan
!
!-----------------------------------------------------------------------
!*** What is the local extent of the variable on the task subdomain?
Expand Down Expand Up @@ -6258,6 +6340,7 @@ subroutine write_full_fields(Atm)
nz=size(fields_tracers(1)%ptr,3)
!
allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
global_field=real_snan
!
!-----------------------------------------------------------------------
!*** What is the local extent of the variable on the task subdomain?
Expand Down Expand Up @@ -6614,10 +6697,10 @@ subroutine exch_uv(domain, bd, npz, u, v)
! buf1 and buf4 must be of the same size (sim. for buf2 and buf3).
! Changes to the code below should be tested with debug flags
! enabled (out-of-bounds reads/writes).
allocate(buf1(1:24*npz))
allocate(buf2(1:36*npz))
allocate(buf3(1:36*npz))
allocate(buf4(1:24*npz))
allocate(buf1(1:24*npz)) ; buf1=real_snan
allocate(buf2(1:36*npz)) ; buf2=real_snan
allocate(buf3(1:36*npz)) ; buf3=real_snan
allocate(buf4(1:24*npz)) ; buf4=real_snan

! FIXME: MPI_COMM_WORLD

Expand Down