bug: wrf geopotential_height using k array rather than single k

🐛 Your bug may already be reported!

Describe the bug

1. build wrf-dart with default compile flags (02)
   use model_mod_check to interpolate QTY_GEOPOTENTIAL_HEIGHT
2. Take a look at the interpolation calculation for fld(1,:), the array k being used to get phb(i,j,k) rather than
   a single value k(1)
3. What actually happened?
   fld(1,: )gives different answers if you have k vs k(1)

Error Message

No error, just different results if you pass k or k(1)
Here is printing zloc, k, uniquek, fld(1,:) with k vs k(1)

interpolating at  Lon/Lat(deg): 358.99400000  43.74100000  Vert:    7.8068000 km
interpolating for "QTY_GEOPOTENTIAL_HEIGHT"
-------------------------------------------------------------

 helen zloc=   22.384799882548812       k=          22 size(k)=           1
 helen uniquek=          22
 helen fld(1,:)=   7598.3411235953572     
 helen fld(1,:)=   7598.3411235953554     
 helen fld(2,:)=   8111.1110924475479     
member     1, SUCCESS with value    ::    7795.6549473842369

using array k

https://github.com/NCAR/DART/blob/a32ab0341f2f3ec1107e6796c498bf636f33f370/models/wrf/model_mod.f90#L2799-L2803

later using scalar k(1)+1

https://github.com/NCAR/DART/blob/a32ab0341f2f3ec1107e6796c498bf636f33f370/models/wrf/model_mod.f90#L2816-L2821

Which model(s) are you working with?

wrf

Screenshots

If applicable, add screenshots to help explain your problem.

Here is the diff (note the lines are off because I have other fixes in for bitwise testing:

@@ -2755,7 +2755,7 @@ else
          ! Adjust zloc for staggered ZNW grid (or W-grid, as compared to ZNU or M-grid)
          zloc = zloc + 0.5_r8
          k = max(1,int(zloc))  ! Only 1 value of k across the ensemble?
-
+         print*, 'helen zloc=', zloc, ' k=', k, 'size(k)=', size(k)
          deallocate(uniquek)
          ! Re-find the unique k values
          ksort = sort(k)
@@ -2775,7 +2775,7 @@ else
                uk = uk + 1
             endif
          enddo
-
+print*, 'helen uniquek=', uniquek
          ! Check to make sure retrieved integer gridpoints are in valid range
          if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_gz ) .and. &
               boundsCheck( j, wrf%dom(id)%polar,      id, dim=2, type=wrf%dom(id)%type_gz ) .and. &
@@ -2801,7 +2801,18 @@ else
                              dx *wrf%dom(id)%phb(lr(1), lr(2), k) ) + &
                        dy *( dxm*wrf%dom(id)%phb(ul(1), ul(2), k)   + &
                              dx *wrf%dom(id)%phb(ur(1), ur(2), k) ) )  / gravity
-            
+print*, 'helen fld(1,:)=', fld(1,:)         
+
+
+
+           fld(1,:) = ( dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) + &
+                       dym*( dxm*wrf%dom(id)%phb(ll(1), ll(2), k(1))   + &
+                             dx *wrf%dom(id)%phb(lr(1), lr(2), k(1)) ) + &
+                       dy *( dxm*wrf%dom(id)%phb(ul(1), ul(2), k(1))   + &
+                             dx *wrf%dom(id)%phb(ur(1), ur(2), k(1)) ) )  / gravity
+print*, 'helen fld(1,:)=', fld(1,:)            

-00 will give same answer for fld(1,:)
-02 will give different answers.

Note there is counting of unique ks but k is used as an index to phb which is the same across the ensemble.
I don't think this is causing problems but is confusing reading the code.

Version of DART

Which version of DART are you using?
You can find the version using git describe --tags
v11.19.0-33-gc625c2a25

Have you modified the DART code?

Yes

Build information

Please describe:

1. mac M
2. GNU Fortran (MacPorts gcc13 13.3.0_2+stdlib_flag) 13.3.0

Start in models/wrf/model_mod.f90 around the interpolation code at the referenced lines, then reproduce the QTY_GEOPOTENTIAL_HEIGHT calculation with model_mod_check under the default and -O2 compile flags. Compare the fld(1,:) results when phb is indexed with k and k(1); done means the indexing is consistent and the observed differing results are explained or corrected.