Global shallow-water run

[miniufo]

Hi, it is a great model here. I've already make a default run run_speedy(). Now I need to customize it. Specifically, I need to:

• modify dynamical core to single-layer shallow-water model,
• change the resolution to a higher one,
• add forcing field (probably in lat/lon gridded nc file), and
• control the dissipation through regular Laplacian diffusion.

Is there any doc or notebook that I can follow? I guess I don't need to change much because everything was set to the earth's atmosphere (which is great to me).

[milankl]

In #main (v0.5 will be released soon) you can

• use the shallow water model by run_speedy(ShallowWater)
• use higher resolution by choosing run_speedy(ShallowWater,trunc=31), whereby 31 (also called T31) is the default resolution containing all spherical harmonics to degree 31, that's about 400km global resolution. Typical higher resolutions are 42,63,85,127,170,255,340,511,682,1023. The video in the readme is done with T1023.
• control the diffusion by setting diffusion=HyperDiffusion() whereby the latter takes arguments as shown here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/fa734053886d0432cb3aec674b95d73b9237690d/src/dynamics/define_diffusion.jl#L1-L20 You could do power=1 for a Laplacian but then you'd also need to increase the time scale otherwise it's overly diffusive. There's no problem in a spectral model to use higher powers, and for the shallow water model I recommend you to use the default parameters but change diffusion=HyperDiffusion(resolution_scaling=0,adaptive=false) (only in the primitive equation that's not sufficient diffusion).
• forcing: with #320 we're planning to change the interface for that. So please hold back for now unless you would like to contribute to resolving that issue. It's supposed to be as easy as forcing=MyForcing and similar to the structure outlined in that issue you would then be able to externally define your forcing MyForcing and just pass it on. For now: What kind of forcing where you considering?

[miniufo]

That is pretty convenient! I should play with these but I am not sure about:

1. What is resolution_scaling? What is T in constant with T? Temperature? I am expecting a layer thickness h for shallow-water model.
2. I use single layer shallow-water model so power_stratosphere should be irrelavent? Curious why it is less power than troposphere? Smoother in stratosphere?
3. Any forcing to momentum would be OK. I just want to input energy to the flow from rest and with some dissipation to reach a statistical steady state, just like the beautiful demo. Later I may need to calculate the energy input to the flow.

[milankl]

1. That's truncation T, in the shallow water model you don't have to change the time scale of the diffusion with resolution, hence resolution_scaling=0. 1.2 there's layer_thickness [km] (default 8.5 I think), as a kw argument for run_speedy, note that this should be chosen considerably higher than the orography, which is about 5-6km at most (himalayas aren't fully resolved obviously)
2. yes, for the primitive equation model one usually needs more diffusion in the stratosphere as a top boundary condition
3. What I anticipate you'll be able to do u_tend[ij] +=, v_tend[ij] += meaning you can just add whatever tendency you want to u,v at grid point ij and the model will then add the other terms to solve the shallow water equations.

[milankl]

because you said lat lon grid. You can run on such a grid, that's in our case called the FullClenshawGrid which you can pass on as an argument Grid=FullClenshawGrid. For an exact spectral transform you'd then need dealiasing=3, which increases the resolution in grid-point space relative to spectral space. For more performance you could check out the octahedral grids we have.

[miniufo]

I guess the grid is used only for the calculation of nonlinear advection term, and the rests are still in spectral space. So I do want a more efficient octahedral grid.

The latlon grid I used is only for specifying forcings or initial conditions, as well as postprocess (most people may get used to this). Is this the default design?

I tried:

Is this due to 1) v0.5 is still not released (how can I update after your release?) and 2) I made a mistake (I am pretty new to Julia)

[navidcy]

@miniufo, I'm not sure if the error is because of your version. But most probably you have the v0.4 installed. You can confirm this by

using Pkg; Pkg.status()

Until v0.5 is tagged you can install the version on main by

Pkg.add(name="SpeedyWeather", rev="main")

and use that?

[miniufo]

Thanks @navidcy, I have v0.4 installed. So waiting for v0.5...

[navidcy]

Cool. You don't need to wait though. Just use Pkg.add(name="SpeedyWeather", rev="main") :)

[miniufo]

Yes, that works fine! So now how to change the period for my run and how to ouput state variables to netcdf files?

[milankl]

n_days is the argument that controls the length of a run, output=true stores netcdf output, at the moment you can change the following output arguments

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/23f3bce634d9f9e8da39656887bb325fa7a826e8/src/default_parameters.jl#L273-L333

[milankl]

Cool. You don't need to wait though. Just use Pkg.add(name="SpeedyWeather", rev="main") :)

julia> ] add https://github.com/SpeedyWeather/SpeedyWeather.jl#main is another way of doing it.

[miniufo]

Everything goes fine. I need to specify the output variables. I tried output_vars=["u", "v", "vor"], but it seems a wrong expression. How to do it right (I'm a pythoner)? Also, for shallow water simulation, I need "h", thickness of the fluid, and is there a variable for this? or do I need other variables like "temp" to calculate "h"?

[navidcy]

output_vars::Vector{Symbol}

so that means it needs to be a vector of symbols. Try

output_vars = [:u, :v, :vor]

[milankl]

We don't use the layer thickness h as a prognostic variable, but it's anomaly with respect to layer_thickness, often called eta in the shallow water equations. However, for compatibility with the other models we have implemented eta is called pres because it takes on the role of pressure in the shallow water system. So

output_vars = [:u,:v,:vor,:pres]

will output eta in meters. you can then add layer_thickness to get h.

[miniufo]

So layer_thickness is 8.5km, defined in default_parameters.jl? Thanks to you guys for these kind help.

[miniufo]

When I only increase the resolution, these show up:

So do I need to change other parameters when increasing the resolution?

[milankl]

Can you paste your code format in ```julia code ```? Pasting a picture doesn't make it easier to read or to copy parts of it for reproducibility

[miniufo]

It is in jupyter notebook, so the output is just plain text. I try this:

MethodError: no method matching _divergence!(::SpeedyWeather.SpeedyTransforms.var"#kernel#42"{Bool, Bool}, ::LowerTriangularMatrix{ComplexF32}, ::LowerTriangularMatrix{ComplexF64}, ::LowerTriangularMatrix{ComplexF64}, ::SpectralTransform{Float32})

Closest candidates are:
  _divergence!(::Any, ::LowerTriangularMatrix{Complex{NF}}, ::LowerTriangularMatrix{Complex{NF}}, ::LowerTriangularMatrix{Complex{NF}}, ::SpectralTransform{NF}) where NF<:AbstractFloat
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:67

Stacktrace:
 [1] curl!(curl::LowerTriangularMatrix{ComplexF32}, u::LowerTriangularMatrix{ComplexF64}, v::LowerTriangularMatrix{ComplexF64}, S::SpectralTransform{Float32}; flipsign::Bool, add::Bool)
   @ SpeedyWeather.SpeedyTransforms C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:26
 [2] curl!(curl::LowerTriangularMatrix{ComplexF32}, u::LowerTriangularMatrix{ComplexF64}, v::LowerTriangularMatrix{ComplexF64}, S::SpectralTransform{Float32})
   @ SpeedyWeather.SpeedyTransforms C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:15
 [3] initial_conditions!(coefs::ZonalJet, progn::PrognosticVariables{Float32, SpeedyWeather.ShallowWaterModel{Float32, SpeedyWeather.CPUDevice}}, model::SpeedyWeather.ShallowWaterModel{Float32, SpeedyWeather.CPUDevice})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\dynamics\initial_conditions.jl:156
 [4] initial_conditions
   @ C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\dynamics\initial_conditions.jl:9 [inlined]
 [5] initialize_speedy(::Type{Float32}, ::Type{ShallowWater}; kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:trunc, :diffusion, :n_days, :output, :output_dt, :output_path, :output_vars), Tuple{Int64, HyperDiffusion, Int64, Bool, Int64, String, Vector{Symbol}}}})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:68
 [6] initialize_speedy
   @ C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:37 [inlined]
 [7] run_speedy(::Type{Float32}, ::Type{ShallowWater}; kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:trunc, :diffusion, :n_days, :output, :output_dt, :output_path, :output_vars), Tuple{Int64, HyperDiffusion, Int64, Bool, Int64, String, Vector{Symbol}}}})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:15
 [8] top-level scope
   @ In[5]:1

[navidcy]

That’s more helpful than screenshot. But I think @milankl was referring to the input, not the output.

Overall, most of times the problem you are facing might be not as simple to answer just by a glance. So people need to be able to rerun your code. It’s quite unlikely that someone might try to look at a screenshot and copy character-by-character to their machine. But if you put code in code mode like above one can copy and paste easily, reproduce your error and/or figure out what the issue was!

[miniufo]

Oops, here are the codes:

using SpeedyWeather
import SpeedyWeather.ShallowWater as ShallowWater
import SpeedyWeather.HyperDiffusion as HyperDiffusion
import SpeedyWeather.Float32 as Float32

run_speedy(Float32, ShallowWater, trunc=127,
           diffusion=HyperDiffusion(resolution_scaling=0,adaptive=false),
           n_days=15, output=true, output_dt=6,
           output_path="D:/Data/speedy", output_vars=[:u,:v,:vor,:pres])

with the outputs:

MethodError: no method matching _divergence!(::SpeedyWeather.SpeedyTransforms.var"#kernel#42"{Bool, Bool}, ::LowerTriangularMatrix{ComplexF32}, ::LowerTriangularMatrix{ComplexF64}, ::LowerTriangularMatrix{ComplexF64}, ::SpectralTransform{Float32})

Closest candidates are:
  _divergence!(::Any, ::LowerTriangularMatrix{Complex{NF}}, ::LowerTriangularMatrix{Complex{NF}}, ::LowerTriangularMatrix{Complex{NF}}, ::SpectralTransform{NF}) where NF<:AbstractFloat
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:67

Stacktrace:
 [1] curl!(curl::LowerTriangularMatrix{ComplexF32}, u::LowerTriangularMatrix{ComplexF64}, v::LowerTriangularMatrix{ComplexF64}, S::SpectralTransform{Float32}; flipsign::Bool, add::Bool)
   @ SpeedyWeather.SpeedyTransforms C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:26
 [2] curl!(curl::LowerTriangularMatrix{ComplexF32}, u::LowerTriangularMatrix{ComplexF64}, v::LowerTriangularMatrix{ComplexF64}, S::SpectralTransform{Float32})
   @ SpeedyWeather.SpeedyTransforms C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\SpeedyTransforms\spectral_gradients.jl:15
 [3] initial_conditions!(coefs::ZonalJet, progn::PrognosticVariables{Float32, SpeedyWeather.ShallowWaterModel{Float32, SpeedyWeather.CPUDevice}}, model::SpeedyWeather.ShallowWaterModel{Float32, SpeedyWeather.CPUDevice})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\dynamics\initial_conditions.jl:156
 [4] initial_conditions
   @ C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\dynamics\initial_conditions.jl:9 [inlined]
 [5] initialize_speedy(::Type{Float32}, ::Type{ShallowWater}; kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:trunc, :diffusion, :n_days, :output, :output_dt, :output_path, :output_vars), Tuple{Int64, HyperDiffusion, Int64, Bool, Int64, String, Vector{Symbol}}}})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:68
 [6] initialize_speedy
   @ C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:37 [inlined]
 [7] run_speedy(::Type{Float32}, ::Type{ShallowWater}; kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:trunc, :diffusion, :n_days, :output, :output_dt, :output_path, :output_vars), Tuple{Int64, HyperDiffusion, Int64, Bool, Int64, String, Vector{Symbol}}}})
   @ SpeedyWeather C:\Users\Yu-Kun Qian\.julia\packages\SpeedyWeather\CuEY5\src\run_speedy.jl:15
 [8] top-level scope
   @ In[5]:1

[milankl]

There was a type instability in the spectral_truncation/_interpolation functions, the PR should fix this.

[milankl]

There is currently no multi-threading for the shallow water system, one could definitely think about the best way to do this, but it's not as trivial.

[navidcy]

btw I doubt that line

import SpeedyWeather.Float32 as Float32

does anything. Float32 is a Julia-defined type, not something that SpeedyWeather defines.

[milankl]

Yes, it surprises me that this works! 🤣 Like we obv use Float32, but we don't actually define it, so I'm surprised that you can import it...

[milankl]

And note that

import SpeedyWeather.ShallowWater as ShallowWater
import SpeedyWeather.HyperDiffusion as HyperDiffusion

isn't necessary either because we already export ShallowWater and HyperDiffusion, so as long as you do using SpeedyWeather they'll imported into the current scope. Only if you'd do import SpeedyWeather everything will be behind SpeedyWeather.

[miniufo]

Not sure if this is somewhat strange in jupyter notebook. I am pretty new to julia so I just tried everything util it works. I should try using SpeedyWeather only later, as the progress bar has not finish yet.

[navidcy]

That’s nothing to do with being run via a notebook or not. Most probably something else was the matter if you had problems.

[miniufo]

I removed all the imports, and it works fine now. I don't know why it complains not finding ShallowWater at the beginning. Right now I start to look a bit closer to the results:

1. The default run seems to be a spindown run of an initial zonal jet stream. Just want to make sure there is no external forcing like radiation or restoring. Also, I would like to verify that the total mass is conserved by $M(t) = \iint p cos(\phi) d\lambda d\phi$, which gives:

indicating a large drop of total mass at first few hours. If there are some forcings or restoring things, how to turn those off?

2. The outputs are defined on Gaussian grid. Can I use uniform lat/lon grids as I need some finite difference calculations? I know I could do it offline, but it would be convenient if the model allows me to do it online. But I still prefer the nonlinear terms being calculated on more efficient grids (just output to uniform lat/lon grids). Is this possible?

[milankl]

1. Thanks for checking this, there's a lot of metrics I've never actually verified, good that someone actually is finally doing this! You are right that the global integral of the interface displacement should be conserved, and I'm happy to see that it's not drifting at least. But a 30% drop in the first hours is ... interesting ... could you check:

• Whether this changes with orography? orography=NoOrography() will disable orography and then the default initial conditions are the Galewsky zonal jet with perturbation that should get unstable

• Whether this changes with the implicit time integration? To filter out gravity waves we use a semi implicit scheme that may not conserve mass, I don't know. The reason I'm not necessarily surprised to see the largest changes in the first hours is that we hit the model at time = 0 with full orography which kicks off a crazy amount of gravity waves. The semi implicit schemes filters them out. You should see them vanishing quickly over the first hours/days when looking at a video of eta/u. You can change the semi implicit scheme from backwards implicit_α = 1 (typle \alpha and hit tab) which is the default to centered implicit implicit_α = 0.5 and see whether this makes a difference.

• If you use another grid for output than for the dynamical core the fields get interpolated and our default interpolation is not conservative. This is will introduce an error, but one that doesn't accumulate over time. So playing around with the grids may give us an idea how big the error is from interpolation, which brings me to the next point:

2. Using Grid= and output_Grid= you can control the grid that is used in the dynamical core and the one that's used for output. Setting output_nlat_half you can also output onto a finer/coarser grid, but default is the same grid size as in the dynamical core. For equi-angle grids you should set output_Grid=FullClenshawGrid, which you could also use for the computations, but note that the spectral transform isn't exact for that grid, you'd need to set dealiasing=3 if you wanted it exact (default is 2 for which only the Gaussian grids are exact).

[milankl]

To add: I have played around with interface relaxation in the shallow water system, but they are currently switched off by default. Meaning there is no forcing/restoring at the moment.

[miniufo]

I use the following codes (slightly modify some parameters) for quick tests:

using SpeedyWeather

run_speedy(Float32, ShallowWater, trunc=62,
           diffusion=HyperDiffusion(resolution_scaling=0,adaptive=false),
           n_days=60, output=true, output_dt=6,
           output_path="D:/Data/speedy", output_vars=[:u,:v,:vor,:pres],
           #orography=NoOrography(),
           implicit_α = 0.5,
           run_id="alpha0.5",
           )

It turns out the total mass for each is:

I feel that some fluctuations could be numerical errors. Maybe it is better to define a variation percentage like ($\Delta p$/8.5km $\times$ 100%)

PS: Is the initial condition also output to the netcdf file?

[milankl]

Yes, the initial condition is the first time steo in output.nc.

Yes, could you actually plot $\iint h dA$, i.e. the globally integrated layer thickness? Because that's actually the quantity that's supposed to be conserved. It's only because $\eta + H = h$ and $H$ the undisturbed height (layer_thickness=8.5km minus mountains) doesn't change with time that also $\eta$ is conserved.

it's good to see that the semi implicit on average doesn't seem to affect the mass (non)-conservation. it's more jagged because there's more gravity waves bouncing around with alpha=0.5. Could you try to run on a grid where you can also directly compute the integral on, that means we can circumvent the non-conservative interpolation in the output? e.g. Grid=FullClenshawGrid, dealiasing=3?

[miniufo]

I use the python codes to calculate the total thickness:

ds = xr.open_dataset('output.nc')
M = ( (ds.pres + 8500) * np.cos(np.deg2rad(ds.lat)) ).sum(['lat','lon'])

The updated plot of total thickness:

I think if the grid is not uniform (like previous Gaussian grid where dy is not uniform), the above sum may not be accurate (should be weighted by dy). For uniform ClenShaw grid, this gives:

Changed slightly, but still has some drop off from initial conditions. Maybe this is OK?

[milankl]

Really appreciated that you look into this!

I use the python codes to calculate the total thickness:

ds = xr.open_dataset('output.nc')
M = ( (ds.pres + 8500) * np.cos(np.deg2rad(ds.lat)) ).sum(['lat','lon'])

Unfortunately that's only correct in the NoOrography case! From Vallis' book

What's conserved in the shallow water system is

$$\iint h \cos(\theta) d\theta d\lambda$$

with the thickness $h$ being the sum of the undisturbed thickness $H$ (that's layer_thickness = 8.5km) plus the interface height $\eta$ (pres in our case) minus the mountains! So you need to subtract the mountains still. You can either set them to zero, or get them by

julia> using SpeedyWeather
julia> p,d,m = initialize_speedy(Grid=FullClenshawGrid,trunc=42);
julia> b = m.boundaries.orography.orography   # on the grid you are using
8064-element, 63-ring FullClenshawGrid{Float32}:
  -46.952774
  -40.0008
  -33.107384
  -26.367054
  -19.869013
  -13.694752
   -7.9159937
   -2.5929527
    2.2268085
    6.509178
    ⋮
 2454.1396
 2474.8992
 2496.397
 2518.4946
 2541.061
 2563.9773
 2587.1372
 2610.448
 2633.83
 2657.216

julia> Matrix(b)   # convert to matrix (full grids only)
128×63 Matrix{Float32}:
  -46.9528    -4.04344    24.5938    -21.7987   …  2519.62    2222.29   2358.36  2680.55
  -40.0008     7.07757    33.3884      2.74256     2650.87    2289.42   2410.5   2703.78
 ...

[milankl]

But your bottom left plot based on the FullClenshawGrid and NoOrography looks like that the mass is conserved to 1e-8 accuracy? Which sounds like machine precision to me with Float32, you could check how much that changes with Float64 too!

[miniufo]

Yes, I forgot to take the variation of orography into account. But here we only care about the temporal variation, not its magnitude. As orography does not change with time, the variation should be OK.

Variation seems to be 1e2, with total ~1e8, the accuracy would be 1e-6 (I guess)?

I tried Float64 (only change the first parameter of run_speedy), but it looks almost the same as before:

Maybe setting up a series of diagnostics for these domain integrals could be helpful in verifying the dynamic core of the model:

• total mass
• total momentum
• total angular momentum
• total energy
• total enstrophy all of which should be conserved when forcings and dissipations are absent.

[milankl]

Awesome, thanks for your work!

Could I ask you to do the same analysis but with Grid=FullClenshawGrid, dealiasing=3 ? Because only then the transforms are exact! I'm wondering whether the sudden drop on the first time step is because of transform errors with the initial conditions. But I'm already very happy ro see that there isn't much of a drift, which I would be concerned about.

[milankl]

Maybe setting up a series of diagnostics for these domain integrals could be helpful in verifying the dynamic core of the model:

* total mass

* total momentum

* total angular momentum

* total energy

* total enstrophy
  all of which should be conserved when forcings and dissipations are absent.

Mass and energy we can asses and I'd appreciate your efforts if you want to! Momentum/angular momentum is trickier, and I have little experience with. But we won't conserve enstrophy or energy as you cannot run the model without dissipation neither is this generally desired as we want some numerical dissipation to keep the model stable. Nevertheless an analysis to see how both change over time would be interesting to see!

[miniufo]

Could I ask you to do the same analysis but with Grid=FullClenshawGrid, dealiasing=3 ? Because only then the transforms are exact! I'm wondering whether the sudden drop on the first time step is because of transform errors with the initial conditions.

I used the following codes:

using SpeedyWeather

run_speedy(Float64, ShallowWater, trunc=62,
           diffusion=HyperDiffusion(resolution_scaling=0,adaptive=false),
           n_days=60, output=true, output_dt=6,
           output_path="D:/Data/speedy", output_vars=[:u,:v,:vor,:pres],
           #Grid=FullClenshawGrid, dealiasing=3,
           #output_Grid=FullGaussianGrid,
           orography=NoOrography(),
           implicit_α = 0.5,
           run_id="OrographyAlpha0.5Float64",
           )

So I should switched Grid=FullClenshawGrid, dealiasing=3 on or off at the same time. The plots with FullClenshawGrid in the title should be also with dealiasing=3 turned on.

Now I need to setup my initial conditions for thickness h and velocity u. Can I load regular lat-lon 2D data fields from netcdf files and passing them to the model?

(EDIT: milankl accidentally edited your post instead of adding a new one, see below. It should be back to what you posted initially)

[milankl]

The plots with FullClenshawGrid in the title should be also with dealiasing=3 turned on.

Awesome, but did you output on the FullGaussianGrid (meaning you uncommented both lines) or on the FullClenshawGrid (which would be the default if not specified) ? Because one includes an interpolation on the output step (from Clenshaw to Gaussian) the other one wouldn't...

Now I need to setup my initial conditions for thickness h and velocity u. Can I load regular lat-lon 2D data fields from netcdf files and passing them to the model?

With #327 we're currently working on restructuring the user interface completely. However, you can already do:

1) Load a file as array

julia> using NCDatasets

julia> ds = NCDataset("run-0001/output.nc");

julia> pres = Float32.(ds["pres"][:,:,end])
96×47 Matrix{Float32}:
 1042.5   1038.0    1043.75   1052.5    1037.75   1009.75   …   913.75    740.0    707.875  764.625  737.125
 1043.0   1040.5    1043.75   1044.0    1027.5    1008.25       892.875   714.125  685.75   750.0    730.125
 1043.0   1040.5    1038.75   1032.0    1018.62   1008.62       880.375   695.0    667.875  736.875  723.0
 1043.0   1039.0    1031.25   1020.0    1012.62   1009.5        876.375   682.125  654.0    725.0    715.875
 1042.5   1036.5    1023.0    1010.5    1009.88   1010.12       878.875   674.25   643.75   714.375  708.75
 1042.0   1033.0    1015.5    1004.88   1010.75   1011.0    …   883.625   669.5    636.375  705.0    701.625
 1041.25  1029.5    1010.0    1004.0    1014.75   1013.38       885.375   666.0    631.5    696.5    694.625
 1040.5   1026.5    1007.0    1007.38   1021.62   1017.75       879.625   662.0    628.375  688.875  687.75
 1039.5   1023.62   1006.62   1014.25   1030.25   1023.62       863.75    656.75   626.75   681.875  681.0
 1038.75  1021.62   1008.5    1023.25   1039.0    1029.5        838.5     650.625  626.25   675.5    674.375
    ⋮                                                ⋮      ⋱                                 ⋮
 1030.0    937.875   800.25    723.5     744.125   840.25      1065.5    1043.75   997.25   917.5    790.25
 1031.5    949.125   826.375   761.5     793.125   897.625     1063.5    1024.0    970.375  903.375  785.875
 1033.25   962.0     859.125   811.875   855.125   956.875     1059.0     998.0    939.25   887.375  781.0
 1034.75   976.125   896.0     869.375   920.0    1004.25   …  1051.75    966.5    905.125  870.125  775.625
 1036.5    990.25    933.75    927.125   977.125  1031.5       1040.0     930.125  869.25   852.0    770.0
 1038.0   1003.75    969.125   978.25   1019.0    1039.25      1022.5     890.5    832.875  833.5    763.875
 1039.5   1015.75    999.375  1017.75   1042.75   1034.0        998.75    849.5    797.625  815.125  757.5
 1040.75  1025.75   1022.25   1042.5    1050.25   1023.75       970.375   809.25   764.375  797.375  750.875
 1041.75  1033.25   1037.0    1053.25   1047.0    1014.75   …   940.75    772.125  734.25   780.375  744.125

I had to add a Float32.(...) here because the original data came as type Matrix{Union{Missing,Float32}}...

2) Wrap it into the grid <:AbstractGrid the data comes on. For lat-lon it would be FullClenshawGrid, etc

julia> pres_grid = FullClenshawGrid(pres)
4512-element, 47-ring FullClenshawGrid{Float32}:
 1042.5f0
 1043.0f0
 1043.0f0
 1043.0f0
 1042.5f0
 1042.0f0
 1041.25f0
 1040.5f0
 1039.5f0
 1038.75f0
    ⋮
  790.25f0
  785.875f0
  781.0f0
  775.625f0
  770.0f0
  763.875f0
  757.5f0
  750.875f0
  744.125f0

3) initialise the model and set the prognostic variables

julia> p,d,m = initialize_speedy();
julia> SpeedyWeather.set_pressure!(p,pres_grid)

[milankl]

At the moment it is not as straight forward to choose u as intial condition as the model works with vorticity, meaning you'd need to calculate vorticity from your u field and transform it to spectral space too. To make this more user-friendly I'll add some functionality in the next days, please bear with me.

[milankl]

Maybe #376 is of interest as it implements forcing and drag for the shallow water model to have a stable climatology

[milankl]

Are there any further issues that remain here? Otherwise I would close this issue. Feel free to reopen or add comments though. Many thanks for all the suggestion/inspirations!