Available resolutions

[milankl]

While it might be nice to have a wide range of resolutions available, in practice this might be difficult for speedy. There might be stability constraints such that we may want to agree on a small subset (n=3...8?) of available resolutions, e.g. T30, T60, ... Speedy's default is

• T30, L8, 96x48 grid points (3.75˚ resolution), 40min time step

It seems that ECMWF used operationally in the past decades

• 1985: T106L16
• 1986: T106L19
• 1991: T213L31
• 1998: T319L31
• 1999: T319L50

Which sounds like a good subset of resolutions to aim for. @tomkimpson @samhatfield what do you think, and what are other constraints we should be worried about?

[samhatfield]

Yes this sounds good to me. I think if you want something above T319L50 then Speedy isn't the right model for you.

The parametrizations will have been tuned to work best at T30 and might not be so great at T319. However, you will still see a general circulation broadly corresponding to the real world and I think for the potential use-cases I can envisage for Speedy.jl (data assimilation, ML emulation, GCM design and training, idealised climate projections etc.) this is fine.

[milankl]

The parametrizations will have been tuned to work best at T30 and might not be so great at T319.

Yes, I agree. Once we reach that stage we could think about adapting some parameterizations, but I believe that can be left for later as currently the focus is not on getting the climate bias-free.

[milankl]

Okay, then, how do we get all these parameters?

T	L	nlon	nlat	dt
30	8	98	46	40min
106	19
213	31
319	50

(and we can then probably mix them a little bit, say T319L31). The model currently defines the vertical levels as

σ_half = [0.0, 0.05, 0.14, 0.26, 0.42, 0.6, 0.77, 0.9, 1.0]

which looks like this

So for more vertical levels should we just fit something sinusoidal to it and pick more values in between? On that note, we also probably need to know which model levels are considered to be "stratosphere"

[samhatfield]

For the horizontal there is some flexibility in the choice of mapping between spectral and gridpoint space. As far as I remember the existing Fortran Speedy code can in principle work with any choice of lat/lon and truncation, provided that nlon = 2*nlat. This is something we can play around with, but as a rough guide I would suggest

T	nlon	nlat
30	96	48
106	320	160
213	640	320
319	800	400

I was inspired by this (note N80 refers to a Gaussian grid with 80 latitudes between pole and equator).

As for the vertical, perhaps you could use this as inspiration.

[milankl]

You probably wanted to swap nlon,nlat for T106/213/319, right?

[samhatfield]

Yes 😫

[milankl]

Okay, even in Untch and Hortal, 2004 which describes the finite element vertical discretisation how the actual choice of vertical levels is done is not discussed apart from

Near the surface the layer thickness is as small as 20 m, and itincreases gradually throughout the troposphere and lower stratosphere to reach 1.5 kmnear 50 hPa. It is kept constant at this value up to about 5 hPa and then grows againtowards the top of the model. The distribution of these 60 levels is shown in Fig. 4.

So I literally believe, that it's someone's experience + wild guess. This is the resolution distribution of a handful of operational setups

(xaxis 0 = model top, xaxis=1 surface; sigma = p/p_surface, i.e. sigma = 0 is pressure 0, sigma = 1 is surface pressure)

For L60/L91/137 you can see that many levels were used for the stratosphere (i.e. low pressures) and surface (i.e. high pressure) such that the curve is more S-like looking. The earlier vertical coordinates had a more evenly distributed vertical, such that the curves looked more like straight lines.

I'm suggesting to find a function (generalised logistic probably) that resembles L31 as default but will put some adjustable parameters such that we can later if needed also have a vertical setup that looks like L60/91/137 (but with nlev levels). I reckon we shouldn't have a big focus on stratosphere/boundary layer in speedy anyway and then we can choose an arbitrary number of vertical levels, which might be nice for the parallelisation. Does this make sense?

[samhatfield]

Yes, makes sense. Provided the code is flexible to the number of levels and the sigma values of each level (and I believe the existing Fortran is flexible, though I never tried any other configuration than the 8 levels given), you can also change your mind later right?

[milankl]

Yes. A 7-parameter fit of a generalised logistic function yields this

which looks good, but also seems a bit like an overkill given the range of vertical coordinates that were used in L31/60/91/137. So maybe something simpler is better?

[milankl]

Adressed with c0e7d5c. This changes speedy's default 8 levels, defined by their 9 sigma half levels, to

#	new	old
0.5	0.0	0.0
1.5	0.073686	0.05
2.5	0.162284	0.14
3.5	0.272612	0.26
4.5	0.408969	0.42
5.5	0.573089	0.6
6.5	0.75429	0.77
7.5	0.913631	0.9
8.5	1.0	1.0

[milankl]

Regarding the horizontal, it seems in the past other decisions have been made, but since T255 the number of grid points in the longitude is always twice the max wavenumber T, and choosing T=2^n-1 makes sense for the FFT

So why don't we just let the user pick only T and then nlon = 2(T+1), nlat=T+1 such that we'll have

T	nlon	nlat
31	64	32
63	128	64
127	256	128
255	512	256
511	1024	512

Questions

• Why is speedy's default T30 and not T31?
• Is there anything different because we use a regular Gaussian grid?
• How do we pick the time step dt?

[samhatfield]

That would also work. I think the deviation from the pattern of 2(T+1), nlat=T+1 shown by the operational system is something to do with the choice of "linear", "quadratic" or "cubic" grid, but I have to admit I don't understand this aspect of the IFS's grid so well. I'm not sure this is relevant to us because Speedy uses a "full Gaussian grid" --- nlon = 2*nlat for all parallels (i.e. not reduced).

• I don't know. There is probably a historical reason.
• See above - I think what you propose now is exactly what's appropriate for a regular ("full") Gaussian grid
• I'm not sure. For now you could pick 40 minutes for T31, as it is in Speedy's T30, and then half it for every doubling in spectral resolution?

[milankl]

Okay it seems there's some leakage in the transformation when changing to T31 and 64x32 nlonxlat. Leakage meaning that the transformation-back-transformation isn't as exact as it was with T30 and 96x48, I'll investigate that. But, the transforms already work with any T and nlon x nlat combinations which makes me happy. At the moment any mismatching resolution spectral vs gridded is just padded with zeros, which I guess makes sense.

[samhatfield]

Interesting. You'd always expect some difference on the order of the machine epsilon right?

[tomkimpson]

Where is the operation that transforms between T31 and 64x32 nlonxlat? If it involves a Fourier transform it is probably just a spectral leakage from FT a discrete signal. Can use some windowing filter to compensate

[milankl]

To spice up this discussion, I found in Bourke 1972 (which seems to have provided a large mathematical background to speedy)

Apparently, they combined T30 with 96 in longitude but 80 in latitude. This apparently comes from the requirement that nlat >= (5J + 1)/2 and nlon >= 3J + 1 for "alias free evaluation of the quadratic terms" (Bourke 1974), with J being the spectral truncation (for T30, J=30). While speedy's default resolution sticks to the longitude condition, it violates the latitudinal from these papers, such that I believe the Legendre transform has improved since Bourke 1972/1974?

[samhatfield]

This is probably related to the choice of truncation - SPEEDY (and all spectral models nowadays) uses a triangular truncation instead of rhomboidal.

[samhatfield]

David Randall has a great book chapter on this by the way: page 409.

[milankl]

Maybe we should split this discussion into horizontal and vertical. For the vertical I just came across this

From Lindzen and Fox-Rabinovitz, 1989 such that it seems that aiming for a scaling whereby 1˚ of horizontal resolution should aapprox correspond to O(1km) in the vertical. This means that speedy's default resolution with O(4˚) and 8 vertical levels has a rather high vertical resolution. I don't know how this relates to spectral vs grid resolution in the truncation between them, but if that means we can actually go to fairly high resolution with fever vertical levels that might be nice.

[milankl]

The horizontal resolution is with #35 now variable and tested successfully for

T	nlon	nlat
30 (old default)	96	48
31 (new default)	96	48
42	128	64
85	256	128
170	512	256
341	1024	512
682	2048	1024

obeying the constraints of triangular truncation: nlon >= 3T+1, nlat >= (3T+1)/2.

nlon,nlat are rounded up from 3T+1 by finding the next larger integer that can be represented as 2^i*3^j, with i>0 but j=0,1 for FFT efficiency. E.g. while T=30 => 3T+1 = 91, nlon=96 is chosen because 96=2^5*3

[milankl]

Available horizontal resolutions described in docs with #65

[milankl]

Manually defined sigma half levels are being implemented with #139

[milankl]

Updated the figure with the vertical levels from Frierson 2006

z = collect(range(1,0,nlev+1))
σ_half = @. exp(-5*(0.05*z + 0.95*z^3))

[milankl]

Created with

L31 = CSV.read(joinpath(path,"L31_phalf_levels.txt"),DataFrame)
L40 = CSV.read(joinpath(path,"L40_phalf_levels.txt"),DataFrame)
L60 = CSV.read(joinpath(path,"L60_phalf_levels.txt"),DataFrame)
L91 = CSV.read(joinpath(path,"L91_phalf_levels.txt"),DataFrame)
L137 = CSV.read(joinpath(path,"L137_phalf_levels.txt"),DataFrame);

p_half31 = Vector(Array(L31)[:,1])
p_half40 = Vector(Array(L40)[:,1])
p_half60 = Vector(Array(L60)[:,1])
p_half91 = Vector(Array(L91)[:,1])
p_half137 = Vector(Array(L137)[:,1]);

# fit based on L?
p0 = 1013.25                       # surface pressure
y = p_half31/p0                    # dependent data
x = Vector(range(0,1,length(y)))   # independent data

# generalised logistic
p0 = Float64[0,1,1,1,1,0,1]

function model(x,p)
    A,K,C,Q,B,M,ν = p
    return @. A + (K-A)/(C+Q*exp(-B*(x-M)))^(1/ν)
end

fit = curve_fit(model,x,y,p0)
pfit = coef(fit)
x = range(0,1,9)
yfit = model(x,Float16.(pfit))
Float16.(pfit)   # round for shorter inclusion in the default parameters

# Frierson 2006 spacing, higher resolution in surface boundary layer and in stratosphere
nlev = 64
z = collect(range(1,0,nlev+1))
σ_half = @. exp(-5*(0.05*z + 0.95*z^3));

 for (p,t) in zip([p_half31,p_half40,p_half60,p_half91,p_half137],
                    ["L31","L40","L60","L91","L137"])
    plot(range(0,1,length(p)),p/p[end],label=t,alpha=.7)
end

yfit[1] = 0
yfit[end] = 1

plot(x,yfit,"ko--",label="GL fit",alpha=.3)
plot(reverse(z),σ_half,"k",label="Frierson, 2006")

xlabel("level/nlev")
ylabel("σ")
grid(alpha=.3)

legend()
tight_layout()
# savefig("vertical_levels.png")

[milankl]

input_data/vertical_coordinate folder removed with b8e49ae