Use of uninitialized values in `fv_mapz_mod.cs_profile` when `iv == -3`

[spencerkclark]

Describe the bug

Debugging the source of non-reproducibility for a configuration of SHiELD turned up the use of uninitialized values in a local array in the cs_profile subroutine of the fv_mapz module. The issue occurs only when iv == -3, which corresponds to one of two possible modes for use when remapping the vertical velocity, the other being iv == -2. I have verified this with additional means, but the issue can be identified by inspecting the lines of code clipped below, considering the case when iv == -3:

GFDL_atmos_cubed_sphere/model/fv_mapz.F90

Lines 1988 to 2071 in e94dc4b

	real gam(i1:i2,km)
	real q(i1:i2,km+1) ! interface values
	real d4(i1:i2)
	real bet, a_bot, grat
	real pmp_1, lac_1, pmp_2, lac_2, x0, x1
	integer i, k, im
	if ( iv .eq. -2 ) then
	do i=i1,i2
	gam(i,2) = 0.5
	q(i,1) = 1.5*a4(1,i,1)
	enddo
	do k=2,km-1
	do i=i1, i2
	grat = delp(i,k-1) / delp(i,k)
	bet = 2. + grat + grat - gam(i,k)
	q(i,k) = (3.*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
	gam(i,k+1) = grat / bet
	enddo
	enddo
	do i=i1,i2
	grat = delp(i,km-1) / delp(i,km)
	q(i,km) = (3.*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / &
	(2. + grat + grat - gam(i,km))
	q(i,km+1) = qs(i)
	enddo
	do k=km-1,1,-1
	do i=i1,i2
	q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
	enddo
	enddo
	else ! all others
	do i=i1,i2
	grat = delp(i,2) / delp(i,1) ! grid ratio
	bet = grat*(grat+0.5)
	q(i,1) = ( (grat+grat)*(grat+1.)*a4(1,i,1) + a4(1,i,2) ) / bet
	gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
	enddo
	if (iv.eq.-3) then !LBC for vertical velocities
	do k=2,km-1
	do i=i1,i2
	d4(i) = delp(i,k-1) / delp(i,k)
	bet = 2. + d4(i) + d4(i) - gam(i,k-1)
	q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
	gam(i,k) = d4(i) / bet
	enddo
	enddo
	do i=i1,i2
	! a_bot = 1. + d4(i)*(d4(i)+1.5)
	!q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
	! / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
	d4(i) = delp(i,km-1) / delp(i,km)
	bet = 2. + d4(i) + d4(i) - gam(i,km-1)
	grat = delp(i,km-1) / delp(i,km)
	q(i,km) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - grat*qs(i) - q(i,k-1) )/bet
	q(i,km+1) = qs(i)
	enddo
	else ! all others
	do k=2,km
	do i=i1,i2
	d4(i) = delp(i,k-1) / delp(i,k)
	bet = 2. + d4(i) + d4(i) - gam(i,k-1)
	q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
	gam(i,k) = d4(i) / bet
	enddo
	enddo
	do i=i1,i2
	a_bot = 1. + d4(i)*(d4(i)+1.5)
	q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
	/ ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
	enddo
	endif
	do k=km,1,-1
	do i=i1,i2
	q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
	enddo
	enddo
	endif

The issue specifically stems from the local gam array, which is declared in line 1988 without any initial assignment. When iv == -3 it gets values for level 1 in line 2025 and levels 2 through km - 1 (inclusive) in line 2034, but is not assigned values for level km prior to being used in line 2068. This leads to undefined / non-reproducible behavior at runtime.

To Reproduce

Run SHiELD in non-hydrostatic mode with fv_core_nml.kord_wz less than zero and abs(fv_core_nml.kord_wz) > 7, e.g. fv_core_nml.kord_wz = -9. This tells the model to use iv == -3 and cs_profile when remapping the vertical velocity:

GFDL_atmos_cubed_sphere/model/fv_mapz.F90

Lines 411 to 424 in e94dc4b

	!3) Map W and density; recompute delz; limit w if needed
	if ( .not. hydrostatic ) then
	! Remap vertical wind:
	if (kord_wz < 0) then
	call map1_ppm (km, pe1, w, ws(is,j), &
	km, pe2, w, &
	is, ie, j, isd, ied, jsd, jed, &
	-3, abs(kord_wz))
	else
	call map1_ppm (km, pe1, w, ws(is,j), &
	km, pe2, w, &
	is, ie, j, isd, ied, jsd, jed, &
	-2, abs(kord_wz))
	endif

Expected behavior

I would expect values to be defined for gam(:,km) when iv == -3 prior to it being used in line 2068.

System Environment

Describe the system environment, include:

• OS: Ubuntu 20.04
• Compiler(s): GNU version 10.5
• MPI type, and version: MPICH 3.3.2
• netCDF Version: netCDF-C 4.7.3; netCDF-Fortran 4.5.2
• Configure options: 64bit

[lharris4]

@spencerkclark you are correct. gam(i,km) was incorrectly not specified in the iv .eq. -2 branch. I believe it should be specified the same way as in the interior (gam(i,k) = d4(i) / bet). This is an oversight on my part.

[lharris4]

@laurenchilutti I haven't been assigned an issue before. Should I test the code and submit a patch? Thanks.

[laurenchilutti]

@lharris4 I did not look closely before assigning this - I had thought this was a PR and I was adding you as a reviewer. Please disregard the assignment.

[lharris4]

@laurenchilutti Thank you.

[spencerkclark]

Thanks @lharris4 — I went ahead and made a PR with your suggested fix, and confirmed that it fixed the reproducibility issue: #302. I'm leaning on your expertise regarding the scientific correctness of this change.