Missing `L_qg(k) = .true.` line in non-hail-aware Thompson

[sael9740]

WRF/phys/module_mp_thompson.F

Lines 3266 to 3303 in f52c197

	if (is_hail_aware) then
	do k = kts, kte
	if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
	L_qg(k) = .true.
	rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
	ng(k) = MAX(R2, (ng1d(k) + ngten(k)*DT)*rho(k))
	rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT)
	rb(k) = MIN(rg(k)/rho(k)/rho_g(1), rb(k))
	idx_bg(k) = MAX(1,MIN(NINT(rg(k)/rho(k)/rb(k) *0.01)+1,NRHG))
	else
	rg(k) = R1
	ng(k) = R2
	rb(k) = R1/rho(k)/rho_g(NRHG)
	idx_bg(k) = idx_bg1
	L_qg(k) = .false.
	endif
	enddo
	else
	do k = kte, kts, -1
	idx_bg(k) = idx_bg1
	enddo
	do k = kte, kts, -1
	if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
	rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
	ygra1 = alog10(max(1.E-9, rg(k)))
	zans1 = 3.0 + 2./7.*(ygra1+8.)
	zans1 = MAX(2., MIN(zans1, 6.))
	N0_exp = 10.**(zans1)
	lam_exp = (N0_exp*am_g(idx_bg(k))*cgg(1,1)/rg(k))**oge1
	lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
	ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
	rb(k) = rg(k)/rho(k)/rho_g(idx_bg(k))
	else
	rg(k) = R1
	ng(k) = R2
	rb(k) = R1/rho(k)/rho_g(NRHG)
	L_qg(k) = .false.
	endif

I was browsing this code and noticed we don't have a L_qg(k) = .true. line for the non-hail aware condition. I'm not sure how impactful this is but it looks like a bug to me just from a software standpoint.

[weiwangncar]

@sael9740 Would L_qg(k) = .false. on line 1948 count?

[sael9740]

The following line would set L_qg(k) =.true. at the start of mp_thompson based on the input values of qg:

WRF/phys/module_mp_thompson.F

Lines 1915 to 1917 in f52c197

	if (qg1d(k) .gt. R1) then
	no_micro = .false.
	L_qg(k) = .true.

But later on it looks like it should be updated based on the tendency of qg:

WRF/phys/module_mp_thompson.F

Lines 3266 to 3303 in f52c197

	if (is_hail_aware) then
	do k = kts, kte
	if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
	L_qg(k) = .true.
	rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
	ng(k) = MAX(R2, (ng1d(k) + ngten(k)*DT)*rho(k))
	rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT)
	rb(k) = MIN(rg(k)/rho(k)/rho_g(1), rb(k))
	idx_bg(k) = MAX(1,MIN(NINT(rg(k)/rho(k)/rb(k) *0.01)+1,NRHG))
	else
	rg(k) = R1
	ng(k) = R2
	rb(k) = R1/rho(k)/rho_g(NRHG)
	idx_bg(k) = idx_bg1
	L_qg(k) = .false.
	endif
	enddo
	else
	do k = kte, kts, -1
	idx_bg(k) = idx_bg1
	enddo
	do k = kte, kts, -1
	if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
	rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
	ygra1 = alog10(max(1.E-9, rg(k)))
	zans1 = 3.0 + 2./7.*(ygra1+8.)
	zans1 = MAX(2., MIN(zans1, 6.))
	N0_exp = 10.**(zans1)
	lam_exp = (N0_exp*am_g(idx_bg(k))*cgg(1,1)/rg(k))**oge1
	lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
	ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
	rb(k) = rg(k)/rho(k)/rho_g(idx_bg(k))
	else
	rg(k) = R1
	ng(k) = R2
	rb(k) = R1/rho(k)/rho_g(NRHG)
	L_qg(k) = .false.
	endif

The line corresponding to line 3269 is missing from the else condition when is_hail_aware = .false.. You can see this same pattern for assigning L_q[rics] in the code right above this:

WRF/phys/module_mp_thompson.F

Lines 3214 to 3264 in f52c197

	do k = kts, kte
	if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
	rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
	nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
	if (.NOT. is_aerosol_aware) nc(k) = Nt_c
	L_qc(k) = .true.
	else
	rc(k) = R1
	nc(k) = 2.
	L_qc(k) = .false.
	endif
	if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
	ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
	ni(k) = MAX(R2, (ni1d(k) + niten(k)*DT)*rho(k))
	L_qi(k) = .true.
	else
	ri(k) = R1
	ni(k) = R2
	L_qi(k) = .false.
	endif
	if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
	rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
	nr(k) = MAX(R2, (nr1d(k) + nrten(k)*DT)*rho(k))
	L_qr(k) = .true.
	lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
	mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
	if (mvd_r(k) .gt. 2.5E-3) then
	mvd_r(k) = 2.5E-3
	lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
	nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
	elseif (mvd_r(k) .lt. D0r*0.75) then
	mvd_r(k) = D0r*0.75
	lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
	nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
	endif
	else
	rr(k) = R1
	nr(k) = R2
	L_qr(k) = .false.
	endif
	if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
	rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
	L_qs(k) = .true.
	else
	rs(k) = R1
	L_qs(k) = .false.
	endif
	enddo

This would then be used for conditional execution of code later on in the routine. Example:

WRF/phys/module_mp_thompson.F

Lines 3744 to 3782 in f52c197

	if (ANY(L_qg .eqv. .true.)) then
	nstep = 0
	do k = kte, kts, -1
	vtg = 0.
	if (rg(k).gt. R1) then
	if (is_hail_aware) then
	xrho_g = MAX(rho_g(1),MIN(rg(k)/rho(k)/rb(k),rho_g(NRHG)))
	afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff
	afall = afall * visco(k)**(1.0-2.0*b_coeff)
	else
	afall = av_g_old
	bfall = bv_g_old
	endif
	vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bfall
	vtgk(k) = vtg
	! Goal: less prominent size sorting
	! the ELSE section below is technically (mathematically) correct:
	if (mu_g .eq. 0) then
	vtg = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) * ilamg(k)**bfall
	else
	vtg = rhof(k)*afall*cgg(8,idx_bg(k))*ogg2 * ilamg(k)**bfall
	endif
	vtngk(k) = vtg
	else
	vtgk(k) = vtgk(k+1)
	vtngk(k) = vtngk(k+1)
	endif
	if (vtgk(k) .gt. 1.E-3) then
	ksed1(4) = MAX(ksed1(4), k)
	delta_tp = dzq(k)/vtgk(k)
	nstep = MAX(nstep, INT(DT/delta_tp + 1.))
	endif
	enddo
	if (ksed1(4) .eq. kte) ksed1(4) = kte-1
	if (nstep .gt. 0) onstep(4) = 1./REAL(nstep)
	endif
	endif

[weiwangncar]

@sael9740 So you're saying that L_qg = .true. should be added, say, somewhere between lines 3288 - 3297? I see it now.

[sael9740]

@sael9740 So you're saying that L_qg = .true. should be added, say, somewhere between lines 3288 - 3297? I see it now.

Yes I believe so. I just wanted to confirm with the microphysics folks.

[weiwangncar]

@sael9740 I just checked the code in v4.3, and a section of code resembling those for snow and ice, etc. is there for graupel. So we are missing L_qg = .true. in this section after double-moment for graupel is added. Do you want to make a PR?

[sael9740]

Yes I can. I noticed one other thing:

WRF/phys/module_mp_thompson.F

Lines 1349 to 1447 in f52c197

	do k = kts, kte
	qv(i,k,j) = qv1d(k)
	qc(i,k,j) = qc1d(k)
	qi(i,k,j) = qi1d(k)
	qr(i,k,j) = qr1d(k)
	qs(i,k,j) = qs1d(k)
	qg(i,k,j) = qg1d(k)
	ni(i,k,j) = ni1d(k)
	nr(i,k,j) = nr1d(k)
	th(i,k,j) = t1d(k)/pii(i,k,j)
	#if ( WRF_CHEM == 1 )
	IF ( wetscav_on ) THEN
	rainprod(i,k,j) = rainprod1d(k)
	evapprod(i,k,j) = evapprod1d(k)
	ENDIF
	#endif
	if (qc1d(k) .gt. qc_max) then
	imax_qc = i
	jmax_qc = j
	kmax_qc = k
	qc_max = qc1d(k)
	elseif (qc1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (qr1d(k) .gt. qr_max) then
	imax_qr = i
	jmax_qr = j
	kmax_qr = k
	qr_max = qr1d(k)
	elseif (qr1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (nr1d(k) .gt. nr_max) then
	imax_nr = i
	jmax_nr = j
	kmax_nr = k
	nr_max = nr1d(k)
	elseif (nr1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (qs1d(k) .gt. qs_max) then
	imax_qs = i
	jmax_qs = j
	kmax_qs = k
	qs_max = qs1d(k)
	elseif (qs1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (qi1d(k) .gt. qi_max) then
	imax_qi = i
	jmax_qi = j
	kmax_qi = k
	qi_max = qi1d(k)
	elseif (qi1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (qg1d(k) .gt. qg_max) then
	imax_qg = i
	jmax_qg = j
	kmax_qg = k
	qg_max = qg1d(k)
	elseif (qg1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (ni1d(k) .gt. ni_max) then
	imax_ni = i
	jmax_ni = j
	kmax_ni = k
	ni_max = ni1d(k)
	elseif (ni1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	endif
	if (qv1d(k) .lt. 0.0) then
	write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
	' at i,j,k=', i,j,k
	CALL wrf_debug(150, mp_debug)
	if (k.lt.kte-2 .and. k.gt.kts+1) then
	write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
	CALL wrf_debug(150, mp_debug)
	qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
	else
	qv(i,k,j) = 1.E-7
	endif
	endif
	enddo

On line 1350 we set qv(i,k,j) = qv1d(k). Later on in the loop we check for negative values and reassign them to qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) (line 1442). The weird thing here is that qv(i,k+1,j) has not been assigned yet so we would be using this value from the previous time step. So we are using the updated value of qv from the lower level and the old one from the higher level. I can't imagine that this was intentional. Could we change this to:

qv(i,k,j) = MAX(1.E-7, 0.5*(qv1d(k-1) + qv1d(k+1)))

[weiwangncar]

@sael9740 It does look strange. Maybe it is in rare cases when qv is < 0? I will get other people to take a look.

[sael9740]

Maybe it is in rare cases when qv is < 0?

Yes this is probably a rare case and the change might not have much of an impact but it is inconsistent and I think it would be good to fix anyways. It does complicate things on our end in terms of keeping the GPU and CPU versions consistent with each other.

[weiwangncar]

@mefrediani @gthompsnWRF Can you review the issues outlined by @sael9740?

[mefrediani]

These two changes (lines 1442 and ~3297) make sense to me, @weiwangncar. Do you agree @gthompsnWRF?