Sites output suggestion

[mvieno]

Hi,
the sites.dat file assume that the EMEP model vertical layers are 20 (KMAX_MID=20).

for example for Auchencorth Moss:
Auchencorth_Moss 55.79323 -3.244781 20

This may create an issue as it is only OK if the model vertical layers are 20,

One option could be to add a "-999" and then assign the KMAX_MID:
Auchencorth_Moss 55.79323 -3.244781 -999! Default surface

I added this as a test near line 338 in Sites_mod.f90
s_gx(n) = ix
s_gy(n) = iy
if(lev==-999) lev=KMAX_MID
s_gz(n) = lev

A better solution may be to add -999 for the surface and the actual altitude w
Auchencorth_Moss 55.79323 -3.244781 -999 ! Default surface
Other_site 55.79323 -3.244781 200 ! Find the nearest layer at that altitude in m

[gitpeterwind]

I thing Dave and Alvaro have been developing the sites. The height are for now hardcoded for 20 levels, yes.

[gitpeterwind]

In addition, the levels are probably defined for the old vertical levels, which are no more in use!

[mifads]

Yes, this is strange now. Ideally I would prefer to change to an altitude-based system, or relative altitude system. I won't manage to do this any time soon though, since I find the coding in that module to be extremely difficult to follow. (The move to netcdf made everything too obscure in my view.)

The idea of using -999 seems much easier though, although sites with high altitude would still need some hard-coded iz values in sites.dat.

[mifads]

In addition, the levels are probably defined for the old vertical levels, which are no more in use!

I am not sure what is meant here. Doesn't the code just use the iz index, and make use of KMID_MAX? Nothing about sigma, eta, or specific pressure levels?

[gitpeterwind]

I am not sure what is meant here.

The surface is "easy", (one could simply replace 20 with kmax, as is done for columns). However the other levels depend on the vertical distribution of levels. We shifted from sigma to EC/IFS levels a few years ago, but I do not think the levels in sites has been updated.

[mifads]

In the f90 code KMAX_MID is used, but think the "20" and other iz values only apply to the sites.dat file, and yes, those were calculated for earlier vertical levels.

Specifying altitude would be a very good alternative, but it is so easily mis-used and mis-understood too. A site at 500m might need vertical profile and/or deposition correction if sitting on a 500m high plateau (hence it has a relative altitude of 0m, and should use iz=KMAX_MID), but if sitting on an isolated mountain in terrain of height 0m, the relative altitude would indeed to 500m and we should pick from some model level which represents that. Tricky! (Best might be to give station altitude and relative altitude in separate columns, so it is explicit at least.)

[gitpeterwind]

For now I have implemented Massimos proposition (but not altitude in meters):

    if(lev<=0) lev = KMAX_MID
    if(lev>KMAX_MID)then
       write(*,*)'WARNING: sites.dat found vertical level out of range. Setting to ',KMAX_MID
    endif
    lev=min(lev,KMAX_MID)```

[mifads]

Hi Peter,
We used to have a system where iz==0 meant use KMAX_MID, but ignore the deposition gradient. So, one could have both 20 and 0 and get different answers. Does that still work? (I am wondering if your <=0 above should just be <0=).

[gitpeterwind]

I did not know that. Now I removed the "=", so at least we have the option for the future.

[mifads]

HI @mvieno
The new rv3.36 allows some new options for the sites output, allowing altitude or relative altitude in meters. I need to update the docs to explain these, but if interested just let me know.

[mvieno]

HI @mvieno
The new rv3.36 allows some new options for the sites output, allowing altitude or relative altitude in meters. I need to update the docs to explain these, but if interested just let me know.

Hi @mifads yes please. 👍

[mifads]

Hi Massimo,
I'll update the docs soon, but not today. Briefly though, one can now set the "Coords" parameter (in the header of the sites.dat input file, used to be just "LatLong") to be one of:

1. LatLonKdown  - same as older LatLong, with lat/lon coordinates, then k-number with ground level being e.g. 20
2. LatLonZm - give lat, long, then altitude in metres. The model will then compare that altitude with the model's topography, to estimate a relative altitude (Hrel). It then calculates which model layer corresponds to that altitude.
3. LatLonHrel - give lat, long, and then your own preferred relative altitude (Hrel)
4. IJKdown - does everything in model's i, j and k coordinates

I actually prefer (3) now, since relative altitudes as usually calculated compared to local topography (e.g. < 5km), and so best calculated using an off-line digital elevation model. Option (3) also allows the same Hrel to be used regardless of the underlying meteo resolution. One can also get Hrel from the TOAR database for some sites.

However, I find that nothing is perfect, so it is good to have the chance to tweak the Hrel so that e.g. the diurnal variations come out best. (Sites with little day-night variation in O3 are probably on hills or similar, and might need higher Hrel to capture that. Beware coastal sites though which can also have no day-night changes despite being at sea-level.)