Why change anything?

So far when ICAR used a forcing dataset where the z-levels did not vary with time it did extrapolate some quantities downward from the lowest forcing level by keeping them constant. This could potentially lead to unwanted behavior as I’m going to show with an example in the following. I’ll refer to the ICAR version before the code modification as “ICAR old” and the version after the modification as “ICAR new”.

An example illustrating the problem

Let’s say the potential temperature in the forcing is defined at elevations of 500m, 1500m, 2500m and so on while the vertical levels of ICAR are defined at 25m, 87.5m, 175m, 325m, 525m, 725m and in 200m spacings from there on. That means there are 4 vertical levels in ICAR below the lowest forcing level.

In these four levels ICAR would then assign a quantity, such as potential temperature, the value of this quantity given at the lowest vertical level in the forcing (blue curve in the following plot). However, the values could also be extrapolated with linear interpolation, which for the reasons below should be the preferred way to do so.

Potential temperature in ICAR before the code modification (icar old), after the code modification (icar new) and in the forcin

For a dry atmosphere the Brunt Väisälä frequency is given by

N^2 = g \frac{d \ln\Theta(z)}{z}

where g is the acceleration due to gravity, \Theta the potential temperature and z the height above ground.

Evaluating this expression for the ICAR potential temperature profiles in the plot above and plotting them reveals the problem that occurs:

Profiles of the Brunt Väisälä frequency in ICAR before the code modification (icar old), after the code modification (icar new)

If there’s no gradient in \Theta then N is equal to zero! However, ICAR enforces a minimal value of N so it will never be actually zero but some value specified in the options file (the default is 3.2\cdot 10^{-4} \text{s}^{-1}). While this is certainly not wrong it’s probably not the behavior you’d expect.

Since the downscaled wind-field also depends on the value of N there are some small differences in the w field depending on how the extrapolation is handled. Here’s an example for an idealized simulation over a Witch of Agnesi ridge (U=10m/s, \Theta profile chosen so that N=0.01s^{-1} everywhere, a_0=1800m, a_1=13200m):

Cross section of the vertical wind-speed in ICAR for and idealized simulation with a Witch of Agnesi ridge topography. Left: ICAR with linear interpolation of quantities to lower levels, right: ICAR without linear interpolation

The structure of the w-field for this scenario remains similar. However, clearly differences are visible to the lee of the ridge (centered at x=0) where the ICAR version with downward linear extrapolation (left) is much smoother.

ICAR code modification

Changes in the ICAR code are limited to the generation of the lookup table for vertical interpolation. The routine is called vLUT and is found in src/utilities/vinterp.f90.

The code changes between ICAR old and ICAR new may be viewed on github. It will take a while to get the changes into the official ICAR repository at NCAR. For that to happen I’ll add additional options so that the applied change is optional and not the standard behavior. Currently I don’t have much time to do so, but it will happen eventually. I’ll add an edit here once it’s done. Until then you can pull the changes I’ve made from my repository of ICAR.