bug: when localizing in height, the mpas_atm model_mod may compute an incorrect vertical location

[nancycollins]

Describe the bug

In reviewing the mpas_atm model_mod code, the subroutine find_vert_indices() calls find_vert_level() which calls find_height_bounds(). find_height_bounds uses either the zGridCenter array or the zGridEdge for vertical locations. it returns the lower and upper model level numbers that enclose that vertical location, along with the fraction of the level.

in convert_vert_distrib() if the requested outgoing vertical is VERTISHEIGHT, after computing the level numbers based on cell centers, it indexes into the zGridFace array for the heights instead of zGridCenter or zGridEdge.

this seems wrong and needs to be tested - possibly with a version of model_mod_check which includes a call to convert_vert.

note that if W interpolation is added, this code will need to distinguish between cell centered data in the vertical and data on the cell faces, where the existing code might be correct for the new case.

To reproduce/test the bug

set the namelist so it will localize in the vertical on height.
test a location of a variable located on the cell center (T, Q) in the vertical on a known level (VERTISLEVEL)
look at the location after vertical conversion. the height should match the height of the cell center, not the cell face.

Error Message

This will not provoke any error messages, but would increase the distances used in localization by half a vertical cell.

Which model(s) are you working with?

mpas_atm

Version of DART

This code is in the main branch

Have you modified the DART code?

there is a modified version that includes other changes which is not in the dart repo.

Build information

This was discovered by inspecting the code. The tests can be run anyplace but based on the size of the mpas test cases, it probably needs to be run on cheyenne.

[hkershaw-brown]

@nancycollins are these the lines with the bug?

DART/models/mpas_atm/model_mod.f90

Lines 4689 to 4690 in f996f26

	fdata(i, :) = zGridFace(k_low(i, :),c(i))*(1.0_r8 - fract(i, :)) + &
	zGridFace(k_up (i, :),c(i))*fract(i, :)

If not can you give a link to the lines:

[nancycollins]

yep, those are them.

[nancycollins]

i'm sorry i didn't make this clear in the bug report. i am intending to do the test once cheyenne comes back up (all the test case code is there).

[hkershaw-brown]

Reproducer for this available at
/glade/derecho/scratch/hkershaw/DART/Bugs/bug_524/runs_524
one obs in height, localizing in height

You can run this on your laptop, on Derecho you'll need an interactive job (killed for memory on login nodes).

Code to show the bug:
https://github.com/hkershaw-brown/DART/tree/bug_524_reproducer

This code skips the early return for ztypein == ztypeout == VERTISHIEGHT

zout should == zin (converting height to height)

 PE 0: filter: Ready to assimilate up to       1 observations
 HK skipping early return ztypein == ztypeout == VERTISHIEGHT
 HK verttype           3
 HK location lon   309.111074247756     
 HK location lat   38.6883872615110     
 HK location vert   12045.3370000000     
 HK convert_vert_distrib k_low, k_up          35          35          35
          36          36          36
 HK zGridFace(k_low(i, :)   11317.0703125000     
 HK zGridFace(k_up(i, :)   12045.3369140625     
 HK zGridCenter(k_low(i, :)   11681.2036132812     
 HK zGridCenter(k_up(i, :)   12423.4858398438     
 HK zin(:)    12045.3370000000     
 HK zout(:)   11674.3282054662     
 PE 0: filter_assim: Processed       1 total observations

with zGridCenter

@@ -4700,8 +4700,8 @@ print*, "HK zGridCenter(k_up(i, :)", zGridCenter(k_up(1, :), c(1))
    fdata = 0.0_r8
    do i = 1, n
       where (istatus == 0)
-         fdata(i, :) = zGridFace(k_low(i, :),c(i))*(1.0_r8 - fract(i, :)) + &
-                       zGridFace(k_up (i, :),c(i))*fract(i, :)
+         fdata(i, :) = zGridCenter(k_low(i, :),c(i))*(1.0_r8 - fract(i, :)) + &
+                       zGridCenter(k_up (i, :),c(i))*fract(i, :)
       end where
    enddo

HK skipping early return ztypein == ztypeout == VERTISHIEGHT
 HK verttype           3
 HK location lon   309.111074247756     
 HK location lat   38.6883872615110     
 HK location vert   12045.3370000000     
 HK convert_vert_distrib k_low, k_up          35          35          35
          36          36          36
 HK zGridFace(k_low(i, :)   11317.0703125000     
 HK zGridFace(k_up(i, :)   12045.3369140625     
 HK zGridCenter(k_low(i, :)   11681.2036132812     
 HK zGridCenter(k_up(i, :)   12423.4858398438     
 HK zin(:)    12045.3370000000     
 HK zout(:)   12045.3370000000     
 PE 0: filter_assim: Processed       1 total observations

This is obs_seq.one, can change to other locations as needed to test:

observation location is cell 35, zgrid lev 35 = 12045.337m
lon =  5.395006 rad
lat =  0.67523974 rad
lon =  309.11105 deg
lat =  38.688385 deg
height zgrid[lev-1] = 11317.07
height zgrid[lev] = 12045.337
height zgrid[lev+1] = 12801.635