Tagging @mdtoyNOAA @grantfirl @rhaesung @Qingfu-Liu @dustinswales

@climbfuji Thanks for identifying this. Sure seems like a bug to me.
@mdtoy Thoguhts?

This is definitely a bug. However, this code is in subroutine 'gwdps_v0' that is only called by the CCPP scheme 'cires_ugwp' when 'do_ugwp=true'. I see there are still some suites (e.g., FV3_GFS_v16, FV3_RRFS_v1beta) that include this scheme. The later suites use the 'unified_ugwp' and 'ugwpv1_gsldrag' schemes, which do not call 'gwdps_v0', but instead call 'gwdps_run' or 'drag_suite', which are free of the corresponding bug.

(One caveat is that the bug exists on lines 437-438 of module 'cires_ugwpv1_oro':
437 bnv2(i,k+1) = max( bvf2, bnv2min )
438 ri_n(i,k+1) = bnv2(i,k)/shr2 ! richardson number consistent with bnv2
But these lines are only called if 'do_ugwp_v1_orog_only' = true, which is false by default and has never been set to true in practice in any workflow that I know of.)

I am planning on doing a code clean-up eventually, and these subroutines that contain the bug will be removed.

@mdtoyNOAA Thanks. Unfortunately, the NRL code still uses this. I am seeing this a few lines further down:

      k = 1
      DO I = 1, npt
        bnv2(i,k) = bnv2(i,k+1)
      ENDDO

This fixes the NaN for bnv2 for k=1, but by that time it has already been assigned to RI_N(i,k+1) (and that one doesn't get fixed.

Would it be ok to move the calculation of RI_N to its own k-loop after bnv2 (i,1) is set to bnv2(i,2)?

@mdtoyNOAA Thanks. Unfortunately, the NRL code still uses this. I am seeing this a few lines further down:

      k = 1
      DO I = 1, npt
        bnv2(i,k) = bnv2(i,k+1)
      ENDDO

This fixes the NaN for bnv2 for k=1, but by that time it has already been assigned to RI_N(i,k+1) (and that one doesn't get fixed.

Would it be ok to move the calculation of RI_N to its own k-loop after bnv2 (i,1) is set to bnv2(i,2)?

Or, even better:

RI_N(i,k+1) = Bnv2(i,max(k,2))/SHR2

? This way:

RI_N(i,2) = Bnv2(i,2)/SHR2
RI_N(i,3) = Bnv2(i,2)/SHR2
RI_N(i,4) = Bnv2(i,3)/SHR2
...

@climbfuji It looks like

@mdtoyNOAA Thanks. Unfortunately, the NRL code still uses this. I am seeing this a few lines further down:

      k = 1
      DO I = 1, npt
        bnv2(i,k) = bnv2(i,k+1)
      ENDDO

This fixes the NaN for bnv2 for k=1, but by that time it has already been assigned to RI_N(i,k+1) (and that one doesn't get fixed.
Would it be ok to move the calculation of RI_N to its own k-loop after bnv2 (i,1) is set to bnv2(i,2)?

Or, even better:

RI_N(i,k+1) = Bnv2(i,max(k,2))/SHR2

? This way:

RI_N(i,2) = Bnv2(i,2)/SHR2
RI_N(i,3) = Bnv2(i,2)/SHR2
RI_N(i,4) = Bnv2(i,3)/SHR2
...

This wouldn't work because SHR2 is calculated at each level in the k-loop, so it would have to become a stored vector, i.e., SHR2(k)

yes, therefore the second suggestion to use max(k,2) in the current place in the code

@climbfuji Oops! Sorry, I got mixed up on the order of your suggestions. Yes, that's the way to do it:
RI_N(i,k+1) = Bnv2(i,max(k,2))/SHR2

Thanks! I'll make the change in the NRL fork, and if it all works as expected I'll create the PR for upstream (NCAR ccpp-physics main) to close this issue.