Improved algorithm to compute terrain horizon and Sky View Factor

[ChristianSteger]

Algorithm

Description:
The algorithm to compute terrain horizon (and derived Sky View Factor) for ICON is replaced by a C++ routine that used the high-performance ray tracing library Embree (and Intel TBB). The general idea is based on the this reference. I've already tested the new algorithm for the official MeteoSwiss ICON domains (1km & 2km resolution) and for various other Alpine domains with resolutions from ~100 m to 2km.

Improvements:

• Accuracy increased:
  
• Correct computation of terrain horizon close to ICON domain boundary with available topography data
  (but still no added margin to compute terrain horizon correctly close to the boundary):
  
• Azimuth shift corrected (before, the terrain horizon value for the first azimuth angle was not for 0.0 deg (North))
• Significant runtime increase (ca. 250 times faster for ~100 m ICON resolution)

Shortcomings of new algorithm:

• In the algorithm, spherical coodinate are converted to Cartesian coordinates (local tangent plane coordinate system) [m]. Ray tracing in Embree is resricted to single precision, implying that coordinates can become unprecise for very large domains (> 5000 km to global), which can cause issues with ray tracing. Increasing the value 'ray_org_elev' to e.g. 0.5 or 1.0 m can mitigate this problem but for extremely large domain (e.g., global), a domain decomposition would be necessary.

Remaining questions:

• make 'ray_org_elev' a customisable parameter in EXTPAR?
• Keep old Fortran code to compute horizon or remove?
• How do handle two new external dependencies (see below)? Always install them or only if the user needs the terrain horizon computation feature of EXTPAR?

Dependencies, compilation and linking of code

The new code depends on two external libraries, Intel Embree and TBB. For testing, I install the binaries of these dependencies manually but they can also be installed via various package managers (e.g., APT).

Install Embree

cd /scratch/mch/csteger/ExtPar/tbb_embree/
mkdir embree
cd embree
wget https://github.com/embree/embree/releases/download/v4.4.0/embree-4.4.0.x86_64.linux.tar.gz
tar xzf embree-4.4.0.x86_64.linux.tar.gz
rm embree-4.4.0.x86_64.linux.tar.gz
cd ..

embree_inc=/scratch/mch/csteger/ExtPar/tbb_embree/embree/include
embree_lib=/scratch/mch/csteger/ExtPar/tbb_embree/embree/lib

Install TBB

cd /scratch/mch/csteger/ExtPar/tbb_embree/
wget https://github.com/uxlfoundation/oneTBB/releases/download/v2022.1.0/oneapi-tbb-2022.1.0-lin.tgz
tar xzf oneapi-tbb-2022.1.0-lin.tgz
rm oneapi-tbb-2022.1.0-lin.tgz

tbb_inc=/scratch/mch/csteger/ExtPar/tbb_embree/oneapi-tbb-2022.1.0/include
tbb_lib=/scratch/mch/csteger/ExtPar/tbb_embree/oneapi-tbb-2022.1.0/lib/intel64/gcc4.8

For testing purposes, I also didn't properly adapt the Makefile yet and I compiled and linked the new code with the following workaround:

cd /scratch/mch/csteger/ExtPar/extpar

./configure.balfrin.gcc
source modules.env

# Compile C++ code (with Embree & TBB)
g++ -I${embree_inc} -L${embree_lib} -lembree4 -I${tbb_inc} -L${tbb_lib} -ltbb -c src/mo_lradtopo_horayzon.cpp -o src/mo_lradtopo_horayzon.o

make -j 16 # aborts during linking of 'extpar_topo_to_buffer.exe'

# Manual link 'extpar_topo_to_buffer.exe' (link library paths in executable)
gfortran -lstdc++ -I${embree_inc} -L${embree_lib} -Wl,-rpath,${embree_lib} -lembree4 -I${tbb_inc} -L${tbb_lib}  -Wl,-rpath,${tbb_lib} -ltbb -o bin/extpar_topo_to_buffer.exe -Imod -Jmod -Ibundled/cdi/src  -I/mch-environment/v6/linux-sles15-zen3/gcc-11.3.0/netcdf-fortran-4.5.4-qpkuzjpjr65anrkv5fenbitri4d3o72q/include  -cpp -Wall -pedantic -fbacktrace -O3 -g -ffree-line-length-256 -fopenmp -L/mch-environment/v6/linux-sles15-zen3/gcc-11.3.0/netcdf-c-4.8.1-ermcs67iouhyo7xlxk5xvetkzdyorqce/lib  -pthread -Wl,-rpath -Wl,/mch-environment/v6/linux-sles15-zen3/gcc-11.3.0/netcdf-c-4.8.1-ermcs67iouhyo7xlxk5xvetkzdyorqce/lib src/extpar_topo_to_buffer.o src/mo_lradtopo.o src/mo_topo_output_nc.o src/mo_var_meta_data.o src/mo_python_data.o src/mo_soil_data.o src/mo_terra_urb.o src/mo_agg_sgsl.o src/mo_sgsl_routines.o src/mo_agg_topo_icon.o src/mo_search_icongrid.o src/mo_target_grid_routines.o src/mo_read_extpar_namelists.o src/mo_icon_grid_routines.o src/mo_additional_geometry.o src/mo_icon_grid_data.o src/mo_icon_domain.o src/mo_preproc_for_sgsl.o src/mo_agg_topo_cosmo.o src/mo_physical_constants.o src/mo_bilinterpol.o src/mo_oro_filter.o src/mo_search_ll_grid.o src/mo_lradtopo_horayzon.o src/mo_topo_sso.o src/mo_cosmo_grid.o src/mo_target_grid_data.o src/mo_topo_routines.o src/mo_base_geometry.o src/mo_math_constants.o src/mo_utilities_extpar.o /scratch/mch/csteger/ExtPar/extpar/src/info_extpar.o src/mo_topo_data.o src/mo_topo_tg_fields.o src/mo_array_cache.o src/mo_util_mmap_cache.o src/util_mmap_cache.o src/mo_io_utilities.o src/mo_logging.o src/mo_grid_structures.o src/mo_io_units.o src/mo_kind.o bundled/cdi/src/.libs/libcdi_f2003.a bundled/cdi/src/.libs/libcdi.a -L/mch-environment/v6/linux-sles15-zen3/gcc-11.3.0/netcdf-fortran-4.5.4-qpkuzjpjr65anrkv5fenbitri4d3o72q/lib -Wl,-rpath -Wl,/mch-environment/v6/linux-sles15-zen3/gcc-11.3.0/netcdf-fortran-4.5.4-qpkuzjpjr65anrkv5fenbitri4d3o72q/lib -lnetcdff  -lnetcdf  -lm  -luuid

[ChristianSteger]

One additional small problem just came to my mind (not sure if it ever occurs but it could lead to strange behaviour). If a ICON grid cell circumcenter aligns perfectly with either the North or South Pole, it is obviously not possible to define the aspect (slope azimuth) angle for this point .This information is, together with the slope angle, required for the radiation-topography-correction scheme. I think there is a similar quantity in the SSO scheme that is expressed as an azimuth angle. However, it seems that for a global grid, the ICON vertices typically align with the Poles (and no the circumcenters) and for the regional domain, it is very unlikely that his happens - so it might not be necessary to address this issue...