get_dist for VERTISSURFACE

[hkershaw-brown]

Moved from #756

location_mod get_dist "! Vert distance can only be done for like vertical locations units"
why can you not calculate vertical distance between two points one with VERTISSURFACE(m) and one with VERTISHEIGHT(m)?

DART/assimilation_code/location/threed_sphere/location_mod.f90

Lines 507 to 509 in e218864

	if(.not. comp_h_only) then
	! Vert distance can only be done for like vertical locations units
	if(loc1%which_vert /= loc2%which_vert) then

Originally posted by @hkershaw-brown in #756 (comment)

[hkershaw-brown]

copied comments over:

nancycollins commented 1 hour ago
it's because the vertical conversion is almost always model-dependent. the location mod expects the model mod to do the conversion before calling get_dist.

@nancycollins
Member
nancycollins commented 1 hour ago
p.s. in this case the location mod can't know the height of the surface to do the calculation.

@hkershaw-brown
Member
Author
hkershaw-brown commented 35 minutes ago
VERTISSURFACE is defined as surface value (value is surface elevation in m) so no conversion needed between VERTISSURFACE and VERTISHEIGHT

https://docs.dart.ucar.edu/en/latest/guide/creating-obs-seq-real.html#observation-locations

[hkershaw-brown]

this text in the docs:

DART/guide/creating-obs-seq-real.rst

Lines 333 to 337 in e218864

	When observations are declared to be ``VERTISSURFACE`` or ``VERTISUNDEF``
	it is not possible to compute a vertical distance between the observation and
	anything else. Consequently, the distance between that observation and everything
	else (state, other observations) is strictly a horizontal distance, and the observation
	will impact the entire column (all levels) within the horizontal localization radius.

Is this true for VERTISSURFACE? That you would want want anything surface to impact the whole column?
This is not what the code is doing.

DART/assimilation_code/location/threed_sphere/location_mod.f90

Lines 501 to 520 in e218864

	! If which_vert has no vertical definition for either location do only horizontal
	if(loc1%which_vert == VERTISUNDEF .or. loc2%which_vert == VERTISUNDEF) comp_h_only = .true.
	! If both verts are surface, do only horizontal
	if(loc1%which_vert == VERTISSURFACE .and. loc2%which_vert == VERTISSURFACE) comp_h_only = .true.
	! Add in vertical component if required
	if(.not. comp_h_only) then
	! Vert distance can only be done for like vertical locations units
	if(loc1%which_vert /= loc2%which_vert) then
	write(msgstring,*)'loc1%which_vert (',loc1%which_vert, &
	') /= loc2%which_vert (',loc2%which_vert,')'
	call error_handler(E_MSG, 'get_dist', msgstring, source)
	call write_location(logfileunit,loc1)
	call write_location(logfileunit,loc2)
	call write_location(0,loc1)
	call write_location(0,loc2)
	call error_handler(E_ERR, 'get_dist', &
	'Dont know how to compute vertical distance for unlike vertical coordinates', &
	source)
	endif

[hkershaw-brown]

mpas example https://github.com/hkershaw-brown/DART/tree/mpas_surf_height

[nancycollins]

i don't believe that DART itself requires that the vertical value for VERTISSURFACE be in meters. i think it's ignored. WRF and then MPAS were then free to require that observations set the value to be the elevation of the observing station so they could do rejections of obs at too different elevations. you might want to grep the code because i may be confused, but i don't think we require that the vertical value be anything when the vertical type is VERTISSURFACE. that doesn't mean we can't change it if that's useful.

[hkershaw-brown]

so what makes you say "i don't believe that DART itself requires that the vertical value for VERTISSURFACE be in meters'? Was the design originally no value, and people 'hacked' a use for it to reject surface observations?

Looking at the code and docs there are a bunch of places where VERTISSURFACE is m is assumed (or according to comments, defined).

DART/assimilation_code/location/threed_sphere/location_mod.f90

Line 63 in e218864

integer, parameter :: VERTISSURFACE = -1 ! surface value (value is surface elevation in m)

DART/assimilation_code/location/threed_sphere/location_mod.f90

Lines 896 to 902 in e218864

	else if(location%which_vert == VERTISSURFACE ) then
	! most 3d sphere users want height in meters, not pressure.
	! original code asked for pressure:
	!write(*, *) 'Vertical coordinate surface pressure (in hPa)'
	!location%vloc = 100.0 * location%vloc ! only applies to pressure
	write(*, *) 'Vertical coordinate height'
	read(*, *) location%vloc

https://docs.dart.ucar.edu/en/latest/guide/vertical-conversion.html

DART/assimilation_code/programs/obs_selection/obs_selection.f90

Lines 100 to 106 in e218864

	! max differences allowed when deciding a location is the same
	real(r8) :: latlon_tolerance = 0.000001 ! horizontal, degrees
	real(r8) :: surface_tolerance = 0.0001 ! vertical, meters
	real(r8) :: pressure_tolerance = 0.001 ! vertical, pascals
	real(r8) :: height_tolerance = 0.0001 ! vertical, meters
	real(r8) :: scaleheight_tolerance = 0.001 ! vertical, pressure ratio
	real(r8) :: level_tolerance = 0.00001 ! vertical, fractional model levels

It seems prone to bugs to have a converter put out VERTISURFACE where the vertical could be in pressure, height, no value, etc and the model can interpret that number differently.

[nancycollins]

how i remember it is that originally if a location was marked as VERTISSURFACE, the actual vertical value was undefined (meaning unused, not specified as needing to be set, etc). then the wrf folks decided to use it for their own purposes since it was unused. (this was 10-15 years ago.) at that point i probably should have made it part of the spec that the value needed to be vertical in meters, or something. but it was model dependent so i just left it alone. after that point, since i knew it was being used that way, i'm guessing some code migrated into the mainline dart code that assumed meters, without it ever being part of a spec. it wasn't good, but that's the history as i remember it. i agree it's bug prone for obs.

[hkershaw-brown]

ok thanks.
"observation location is model dependent." 🙃