Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reanalysis data as initial conditions #630

Open
rocroc2017 opened this issue Dec 5, 2024 · 22 comments · Fixed by #634
Open

Reanalysis data as initial conditions #630

rocroc2017 opened this issue Dec 5, 2024 · 22 comments · Fixed by #634
Labels
bug 🐞 Something isn't working initial conditions 🚥 Initial conditions of model simulations question ❓ ❔ We got answers!!

Comments

@rocroc2017
Copy link
Contributor

rocroc2017 commented Dec 5, 2024

Dear Dr. Klöwer,

I am writing to express my admiration for your work on the SpeedyWeather.jl model, which I came across recently.

Your development of SpeedyWeather.jl has truly revolutionized the way we think about traditional climate models. The modular approach you have adopted allows researchers like myself to easily adapt and explore various aspects of climate model. This innovation is particularly commendable, as it significantly enhances the flexibility and usability of climate modeling tools.

I am particularly interested in using reanalysis data as the initial data in SpeedyWeather.jl. However, I am facing some challenges in implementing this. I would greatly appreciate your guidance on how to effectively incorporate reanalysis data into the model as an initial data.

Thank you very much for considering my request. I look forward to the possibility of engaging with your expertise and potentially contributing to the ongoing discourse in climate modeling.

@milankl
Copy link
Member

milankl commented Dec 5, 2024

Thanks for reaching out and your interest! Yes, setting the initial conditions from ERA5 / some form of data assimilation is currently missing in SpeedyWeather and would require further attention. As you may have noticed SpeedyWeather’s current focus is more on idealised simulations, interactivity than on actually being good at the weather forecasting challenge. This can be changed but requires some effort. Many thanks for the praise, good to hear that a different approach to user and developer interface is appreciated.

I guess a list from easy to complex of required features to make this possible

  • exists already: use the set!(simulation, temp=…) function to provide some ERA5 data but interpolate that onto lon-lat-sigma coordinates for all our prognostic variables. This is basically what many models in WeatherBench do too but the successful ones would use some encoder (like NeuralGCM) to allow some machine-learned filtering of the initial conditions because just interpolating things may kick off gravity waves and other perturbations that would reduce the skill of a forecast. Check the documentation if you want to get started with this.

  • Use NeuralGCMs encoder although this would restrict us to use that exact same resolution that encoder is trained on. I think this was TL127, 32 vertical layers or TL255 also 32 levels.

  • Use that same encoder architecture and relearn it for SpeedyWeather. Even with a non-ML model having a learned encoder can provide better initial conditions and a functionality to map between them and a future state from ERA5.

  • Actually write a proper data assimilation algorithm. This can be data on a grid like from ERA5 but also more sparse minimising a loss function between model state and data. Variational data assimilation like at ECMWF would be one option, but that would require a differentiable model (in the works) or Ensemble Kalman filter would be another option of particle filter (just chatted with some people on Monday about this)

@milankl milankl added initial conditions 🚥 Initial conditions of model simulations question ❓ ❔ We got answers!! labels Dec 5, 2024
@rocroc2017
Copy link
Contributor Author

rocroc2017 commented Dec 5, 2024

Thank you for your enthusiastic response. I plan to start with the easy way as your suggestion.
I have implemented a simple code to convert from grid points to complex coefficients with the help of artificial intelligence.
My code is quite simple and does not consider truncation. The input is a real matrix, and the output is a complex matrix of the same size.
After reading the introduction in your document, it seems that you have already achieved similar functionality.
https://speedyweather.github.io/SpeedyWeather.jl/dev/grids/
Could you give me an example.
Thank you very much
Best wishes

using FFTW
function grid_to_spectral(grid_data)
nx, ny = size(grid_data)
spectral_coeffs = fft(grid_data)
return spectral_coeffs
end
grid_data = rand(8, 8)
spectral_coeffs = grid_to_spectral(grid_data)
println("Spectral Coefficients Matrix Size:")
println(size(spectral_coeffs))

@milankl
Copy link
Member

milankl commented Dec 5, 2024

Yes we use transform for that. Please have a look at the SpeedyTransforms section in the documentation

@rocroc2017
Copy link
Contributor Author

rocroc2017 commented Dec 10, 2024

Dear Dr. Klöwer
I successfully replaced the initialized vor with a new vor generated using randomly generated arrays by converting wind speed data, and the program is now functional. However, I am uncertain whether this approach is correct or if there are more convenient methods available. Additionally, I have several questions that I would appreciate your guidance on:

  1. Coordinate Conversion: Typically, reanalysis data is provided in pressure (p) coordinates, while SpeedyWeather employs sigma coordinates. Is there a module in Julia that facilitates the conversion from p coordinates to sigma coordinates?
  2. Divergence Calculation: Is there a function similar to curl for calculating divergence (div) in SpeedyWeather?
    Thank you very much for your anticipated response. I look forward to your insights.
using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, Grid=FullGaussianGrid)
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)

#############################################
#change the arrays in simulation.prognostic_variables
##############################################

u = rand(2, 48, 96, 8)  #[nstep, lat ,lon, lev]
v = rand(2, 48, 96, 8)
u = permutedims(u, (1, 3, 2, 4))
v = permutedims(v, (1, 2, 3, 4))
vor= simulation.prognostic_variables.vor

for i in 1:2
   for k in 1:8
      u_grid =FullGaussianGrid(u[i,:,:,k], input_as=Matrix)
      v_grid = FullGaussianGrid(v[i,:,:,k], input_as=Matrix)
vor[i][:,k] = curl(u_grid, v_grid, radius=spectral_grid.radius)
   end
end

run!(simulation)

@milankl
Copy link
Member

milankl commented Dec 10, 2024

1. Coordinate Conversion: Typically, reanalysis data is provided in pressure (p) coordinates, while SpeedyWeather employs sigma coordinates. Is there a module in Julia that facilitates the conversion from p coordinates to sigma coordinates?

We don't have such a functionality yet, you are more than welcome to implement it, you essentially only need to reverse this #268

2. Divergence Calculation: Is there a function similar to curl for calculating divergence (div) in SpeedyWeather?
   Thank you very much for your anticipated response. I look forward to your insights.

Yes, it's call divergence. see here https://speedyweather.github.io/SpeedyWeather.jl/dev/speedytransforms/#SpeedyWeather.SpeedyTransforms.divergence-Union{Tuple{Grid},%20Tuple{Grid,%20Grid}}%20where%20Grid%3C:AbstractGridArray

u = rand(2, 48, 96, 8) #[nstep, lat ,lon, lev]

Arrays in Julia are column-major so we generally use lon-lat-X-Y in that order where X, Y can be any non-horizontal dimensions. E.g. vertical or time. If you want to create random data on a grid you can just use rand(FullGaussianGrid, nlat_half) where nlat_half is the number of latitudes on one hemisphere (equator included) which is our resolution parameter for any grid. The 96x48 full Gaussian grid has 24 latitudes on one hemisphere.

In your case of setting initial conditions, you only have to provide initial conditions for the first leapfrog time step.

nlat_half = 24
set!(simulation, u = rand(FullGaussianGrid, 24, 8), v = rand(FullGaussianGrid, 24, 8))

this will leave the 2nd leapfrog time step at zero, as the first time step of the model is a small Euler step anyway, it would overwrite the 2nd leapfrog time step anyway.

for i in 1:2
for k in 1:8
u_grid =FullGaussianGrid(u[i,:,:,k], input_as=Matrix)
v_grid = FullGaussianGrid(v[i,:,:,k], input_as=Matrix)
vor[i][:,k] = curl(u_grid, v_grid, radius=spectral_grid.radius)
end
end

This then just becomes

julia> u, v = rand(FullGaussianGrid, 24, 8), rand(FullGaussianGrid, 24, 8);

julia> divergence(u, v)
560×8 LowerTriangularArray{ComplexF64, 2, Matrix{ComplexF64}}:
  ...

or transform that again.

@rocroc2017
Copy link
Contributor Author

rocroc2017 commented Dec 11, 2024

Dear Dr. Klöwer
I still don't understand how to use set!. When I use Code Segment A, the program runs normally, but when I use Code Segment B, the program throws an error.
Code Segment A:

using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, Grid=FullGaussianGrid)
#48X96
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)
set!(simulation, u = rand(FullGaussianGrid, 24, 8), v = rand(FullGaussianGrid, 24, 8))

Segment B

using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, Grid=FullGaussianGrid)
#48X96
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)
set!(simulation, u = rand(FullGaussianGrid, 24, 8), v = rand(FullGaussianGrid, 24, 8), temp = rand(FullGaussianGrid, 24, 8))

@milankl
Copy link
Member

milankl commented Dec 11, 2024

(please put code into
```julia
(code)
```
blocks for readability and so that one can just copy it with the button in the top right corner)

@milankl
Copy link
Member

milankl commented Dec 11, 2024

Code Segment A

Mhhh that's a bug, you currently need to provide the same number format Float32 as inspectral_grid. This shouldn't be the case, we need to fix this. You can check that it propagated u, v into the prognostic variables by checking the vorticity and divergence (which are the actual prognostic variables inside the model, not u, v.

julia> set!(simulation, u = rand(FullGaussianGrid{Float32}, 24, 8), v = rand(FullGaussianGrid{Float32}, 24, 8))

julia> simulation.prognostic_variables.vor[1]
560×8 LowerTriangularArray{ComplexF32, 2, Matrix{ComplexF32}}:
         0.0+0.0im                 0.0+0.0im                   0.0+0.0im                 0.0+0.0im
  3.77433f-7+0.0im          3.77419f-7+0.0im             3.80322f-7+0.0im          3.82084f-7+0.0im
  1.29421f-8+0.0im         9.57814f-11+0.0im             1.35997f-9+0.0im         3.46717f-10+0.0im
  2.21733f-7+0.0im          2.18971f-7+0.0im             2.21314f-7+0.0im          2.04825f-7+0.0im
  -1.4568f-8+0.0im          1.14279f-8+0.0im            -1.21748f-8+0.0im         -3.39043f-9+0.0im
  1.52779f-7+0.0im          1.80171f-7+0.0im            1.77427f-7+0.0im          1.72589f-7+0.0im

and similar for simulation.prognostic_variables.div[1], 1 is the first leapfrog time step, you don't have to set the second. That stays 0 (overwritten anyway).

Code Segment B

same here

@milankl milankl added the bug 🐞 Something isn't working label Dec 11, 2024
@milankl
Copy link
Member

milankl commented Dec 11, 2024

To bug originates from here as we don't convert grids to the number format of the transform S

specs = isnothing(S) ? transform(grids) : transform(grids, S)

I suggest to add the following in the line before

grids = convert.(eltype(S), grids)

and define

Base.eltype(S::SpectralTransform{NF}) where NF = NF

generally in src/SpeedyTransforms/spectral_transform.jl. Would you like to create a pull request with that?

@milankl
Copy link
Member

milankl commented Dec 11, 2024

Ideally we would have a fully number format flexible code but unfortunately FFTW isn't happy with that and so we need to iron those bits out where it otherwise would complain

@milankl milankl linked a pull request Dec 12, 2024 that will close this issue
@rocroc2017
Copy link
Contributor Author

Ideally we would have a fully number format flexible code but unfortunately FFTW isn't happy with that and so we need to iron those bits out where it otherwise would complain

Thank you for your help. Now It can work using following code

using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, Grid=FullGaussianGrid)
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)
set!(simulation, u = Float32.(rand(FullGaussianGrid, 24, 8)), v = Float32.(rand(FullGaussianGrid, 24, 8)), temp = Float32.(rand(FullGaussianGrid, 24, 8)))

@rocroc2017
Copy link
Contributor Author

Code Segment A

Mhhh that's a bug, you currently need to provide the same number format Float32 as inspectral_grid. This shouldn't be the case, we need to fix this. You can check that it propagated u, v into the prognostic variables by checking the vorticity and divergence (which are the actual prognostic variables inside the model, not u, v.

julia> set!(simulation, u = rand(FullGaussianGrid{Float32}, 24, 8), v = rand(FullGaussianGrid{Float32}, 24, 8))

julia> simulation.prognostic_variables.vor[1]
560×8 LowerTriangularArray{ComplexF32, 2, Matrix{ComplexF32}}:
         0.0+0.0im                 0.0+0.0im                   0.0+0.0im                 0.0+0.0im
  3.77433f-7+0.0im          3.77419f-7+0.0im             3.80322f-7+0.0im          3.82084f-7+0.0im
  1.29421f-8+0.0im         9.57814f-11+0.0im             1.35997f-9+0.0im         3.46717f-10+0.0im
  2.21733f-7+0.0im          2.18971f-7+0.0im             2.21314f-7+0.0im          2.04825f-7+0.0im
  -1.4568f-8+0.0im          1.14279f-8+0.0im            -1.21748f-8+0.0im         -3.39043f-9+0.0im
  1.52779f-7+0.0im          1.80171f-7+0.0im            1.77427f-7+0.0im          1.72589f-7+0.0im

and similar for simulation.prognostic_variables.div[1], 1 is the first leapfrog time step, you don't have to set the second. That stays 0 (overwritten anyway).

Code Segment B

same here

When I add code

set!(simulation, u = Float32.(rand(FullGaussianGrid, 24, 8)), v = Float32.(rand(FullGaussianGrid, 24, 8)))

And check using code

simulation.prognostic_variables.vor[1]

I find vor[1] has changed.

I assume that during initialization, I only need to modify u and v in set! and do not need to calculate vor and div separately. I would like to ask if my idea is correct.

@rocroc2017
Copy link
Contributor Author

Dear Dr. Klöwer
I write a code to run using ncep data, using T10, the code can work, but the result looks not right. Could you give me some suggestion. Thank you very much.

# Initialize spectral grid and model
spectral_grid = SpectralGrid(trunc=10, Grid=FullGaussianGrid)
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)

# Load data from NetCDF files
ds_u = NCDataset("uwnd.2023.nc")
u1 = ds_u["uwnd"][:,:,:,1]

ds_v = NCDataset("vwnd.2023.nc")
v1 = ds_v["vwnd"][:,:,:,1]

ds_t = NCDataset("air.2023.nc")
t1 = ds_t["air"][:,:,:,1]

ds_rh = NCDataset("rhum.2023.nc")
rh1 = ds_rh["rhum"][:,:,:,1]

ds_pres = NCDataset("pres.sfc.2023.nc")
pres1 = ds_pres["pres"][:,:,1]


function bilinear_interpolation(source::AbstractMatrix, target_size::Tuple{Int,Int})
    src_rows, src_cols = size(source)
    target_rows, target_cols = target_size
    target = zeros(eltype(source), target_rows, target_cols)
    row_scale = (src_rows - 1) / (target_rows - 1)
    col_scale = (src_cols - 1) / (target_cols - 1)
    
    for i in 1:target_rows
        for j in 1:target_cols
            src_row = (i - 1) * row_scale
            src_col = (j - 1) * col_scale
            
            r1 = floor(Int, src_row) + 1
            r2 = min(r1 + 1, src_rows)
            c1 = floor(Int, src_col) + 1
            c2 = min(c1 + 1, src_cols)
            
            w_r = src_row - (r1 - 1)
            w_c = src_col - (c1 - 1)
            
            target[i, j] = (1 - w_r) * (1 - w_c) * source[r1, c1] +
                           w_r * (1 - w_c) * source[r2, c1] +
                           (1 - w_r) * w_c * source[r1, c2] +
                           w_r * w_c * source[r2, c2]
        end
    end
    
    return target
end


source = pres1
target_size = (32, 16)
pres2 = bilinear_interpolation(source, target_size)

u2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    u2[:, :, i] = bilinear_interpolation(u1[:, :, i], target_size)
end

v2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    v2[:, :, i] = bilinear_interpolation(v1[:, :, i], target_size)
end

t2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    t2[:, :, i] = bilinear_interpolation(t1[:, :, i], target_size)
end

rh2 = zeros(Float32, 32, 16, 8)
for i in 1:8
    rh2[:, :, i] = bilinear_interpolation(rh1[:, :, i], target_size)
end


using Interpolations

p_c=[1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
sigma_c=collect(0.9375:-0.125:0.0625)

function p_to_sigma(var_p, p_levels, pres_surf, sigma_levels)
    var_p = convert(Array, var_p)
    p_levels = convert(Array, p_levels)
    pres_surf = convert(Array, pres_surf)
    sigma_levels = convert(Array, sigma_levels)
    nlon, nlat, nplev = size(var_p)
    nsiglev = length(sigma_levels)
    var_sigma = zeros(nlon, nlat, nsiglev)
    p_levels_sorted = reverse(p_levels)
    
    for i in 1:nlon
        for j in 1:nlat
            p_sigma = sigma_levels .* pres_surf[i, j]
            var_p_sorted = reverse(var_p[i, j, :])
            itp = interpolate((p_levels_sorted,), var_p_sorted, Gridded(Linear()))
            
            for k in 1:nsiglev
                if p_sigma[k] < p_levels_sorted[1]
                    var_sigma[i, j, k] = var_p_sorted[1]
                elseif p_sigma[k] > p_levels_sorted[end]
                    var_sigma[i, j, k] = var_p_sorted[end]
                else
                    var_sigma[i, j, k] = itp(p_sigma[k])
                end
            end
        end
    end
    
    return var_sigma
end


u3 = p_to_sigma(u2, p_c, pres2, sigma_c)
v3 = p_to_sigma(v2, p_c, pres2, sigma_c)
t3 = p_to_sigma(t2, p_c, pres2, sigma_c)
rh3 = p_to_sigma(rh2, [1000, 925, 850, 700, 600, 500, 400, 300], pres2, sigma_c)

u_grid = Float32.(FullGaussianGrid(u3, input_as=Matrix))
v_grid = Float32.(FullGaussianGrid(v3, input_as=Matrix))
t_grid = Float32.(FullGaussianGrid(t3, input_as=Matrix))
rh_grid = Float32.(FullGaussianGrid(rh3, input_as=Matrix))
pres_grid = Float32.(FullGaussianGrid(pres2, input_as=Matrix))

set!(simulation, u=u_grid, v=v_grid, temp=t_grid, humid=rh_grid, pres=pres_grid)

run!(simulation)

The NCEP data can be download from
https://downloads.psl.noaa.gov/Datasets/ncep.reanalysis/Dailies/pressure/

@milankl
Copy link
Member

milankl commented Dec 12, 2024

p_c=[1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
sigma_c=collect(0.9375:-0.125:0.0625)

the vertical coordinate in SpeedyWeather (and also in ECMWF's model) is sorted top to bottom, so k=1 is the uppermost layer, k=nlayers is the lowermost layer.

t_grid = Float32.(FullGaussianGrid(t3, input_as=Matrix))
rh_grid = Float32.(FullGaussianGrid(rh3, input_as=Matrix))
pres_grid = Float32.(FullGaussianGrid(pres2, input_as=Matrix))

temperature has to be in K, relative humidity in kg/kg, pressure is actually the logarithm of surface pressure. Maybe we should convert that automatically but currently we don't. So put a log.(...) around one of the arrays.

@milankl
Copy link
Member

milankl commented Dec 12, 2024

I assume that during initialization, I only need to modify u and v in set! and do not need to calculate vor and div separately. I would like to ask if my idea is correct.

Yes, if you set u, v then it computes the curl and divergence for you. that happens here

curl!(vor, u_, v_, S; add, radius=geometry.radius)
divergence!(div, u_, v_, S; add, radius=geometry.radius)

@rocroc2017
Copy link
Contributor Author

Dear Dr. Klöwer
Thank you for your help。
I have flipped the vertical axis, replaced relative humidity with specific humidity, taken the logarithm of surface pressure, and the temperature from NCEP is in K, while the wind speed is in m/s, so no unit conversion is needed. There is one more thing: the unit of surface pressure in NCEP is Pa. Should it be converted to hPa? After making these modifications, the model can run, but the results are still incorrect. Could you please provide some more suggestions or suggest where to start checking?

using SpeedyWeather
using NCDatasets

# Initialize spectral grid and model
spectral_grid = SpectralGrid(trunc=10, Grid=FullGaussianGrid)
model = PrimitiveWetModel(spectral_grid)
simulation = initialize!(model)

# Load data from NetCDF files
ds_u = NCDataset("uwnd.2023.nc")
u1 = ds_u["uwnd"][:,:,:,1]

ds_v = NCDataset("vwnd.2023.nc")
v1 = ds_v["vwnd"][:,:,:,1]

ds_t = NCDataset("air.2023.nc")
t1 = ds_t["air"][:,:,:,1]

ds_sh = NCDataset("shum.2023.nc")
sh1 = ds_sh["shum"][:,:,:,1]

ds_pres = NCDataset("pres.sfc.2023.nc")
pres1 = ds_pres["pres"][:,:,1]

function bilinear_interpolation(source::AbstractMatrix, target_size::Tuple{Int,Int})
    src_rows, src_cols = size(source)
    target_rows, target_cols = target_size
    target = zeros(eltype(source), target_rows, target_cols)
    row_scale = (src_rows - 1) / (target_rows - 1)
    col_scale = (src_cols - 1) / (target_cols - 1)
    
    for i in 1:target_rows
        for j in 1:target_cols
            src_row = (i - 1) * row_scale
            src_col = (j - 1) * col_scale
            
            r1 = floor(Int, src_row) + 1
            r2 = min(r1 + 1, src_rows)
            c1 = floor(Int, src_col) + 1
            c2 = min(c1 + 1, src_cols)
            
            w_r = src_row - (r1 - 1)
            w_c = src_col - (c1 - 1)
            
            target[i, j] = (1 - w_r) * (1 - w_c) * source[r1, c1] +
                           w_r * (1 - w_c) * source[r2, c1] +
                           (1 - w_r) * w_c * source[r1, c2] +
                           w_r * w_c * source[r2, c2]
        end
    end
    
    return target
end

source = pres1
target_size = (32, 16)
pres2 = bilinear_interpolation(source, target_size)

u2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    u2[:, :, i] = bilinear_interpolation(u1[:, :, i], target_size)
end

v2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    v2[:, :, i] = bilinear_interpolation(v1[:, :, i], target_size)
end

t2 = zeros(Float32, 32, 16, 17)
for i in 1:17
    t2[:, :, i] = bilinear_interpolation(t1[:, :, i], target_size)
end

sh2 = zeros(Float32, 32, 16, 8)
for i in 1:8
    sh2[:, :, i] = bilinear_interpolation(sh1[:, :, i], target_size)
end

using Interpolations

p_c=[1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
sigma_c=collect(0.9375:-0.125:0.0625)

function p_to_sigma(var_p, p_levels, pres_surf, sigma_levels)
    var_p = convert(Array, var_p)
    p_levels = convert(Array, p_levels)
    pres_surf = convert(Array, pres_surf)
    sigma_levels = convert(Array, sigma_levels)
    nlon, nlat, nplev = size(var_p)
    nsiglev = length(sigma_levels)
    var_sigma = zeros(nlon, nlat, nsiglev)
    p_levels_sorted = reverse(p_levels)
    
    for i in 1:nlon
        for j in 1:nlat
            p_sigma = sigma_levels .* pres_surf[i, j]
            var_p_sorted = reverse(var_p[i, j, :])
            itp = interpolate((p_levels_sorted,), var_p_sorted, Gridded(Linear()))
            
            for k in 1:nsiglev
                if p_sigma[k] < p_levels_sorted[1]
                    var_sigma[i, j, k] = var_p_sorted[1]
                elseif p_sigma[k] > p_levels_sorted[end]
                    var_sigma[i, j, k] = var_p_sorted[end]
                else
                    var_sigma[i, j, k] = itp(p_sigma[k])
                end
            end
        end
    end
    
    return var_sigma
end

u3 = p_to_sigma(u2, p_c, pres2, sigma_c)
v3 = p_to_sigma(v2, p_c, pres2, sigma_c)
t3 = p_to_sigma(t2, p_c, pres2, sigma_c)
sh3 = p_to_sigma(sh2, [1000, 925, 850, 700, 600, 500, 400, 300], pres2, sigma_c)

u5 = sort(u3, dims=3, rev=true)
v5 = sort(v3, dims=3, rev=true)
t5 = sort(t3, dims=3, rev=true)
sh5 = sort(sh3, dims=3, rev=true)
pres5 = log.(pres2)

u_grid = Float32.(FullGaussianGrid(u5, input_as=Matrix))
v_grid = Float32.(FullGaussianGrid(v5, input_as=Matrix))
t_grid = Float32.(FullGaussianGrid(t5, input_as=Matrix))
sh_grid = Float32.(FullGaussianGrid(sh5, input_as=Matrix))
pres_grid = Float32.(FullGaussianGrid(pres5, input_as=Matrix))

set!(simulation, u=u_grid, v=v_grid, temp=t_grid, humid=sh_grid, pres=pres_grid)

run!(simulation)
p_c=[1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
sigma_c=collect(0.9375:-0.125:0.0625)

the vertical coordinate in SpeedyWeather (and also in ECMWF's model) is sorted top to bottom, so k=1 is the uppermost layer, k=nlayers is the lowermost layer.

t_grid = Float32.(FullGaussianGrid(t3, input_as=Matrix))
rh_grid = Float32.(FullGaussianGrid(rh3, input_as=Matrix))
pres_grid = Float32.(FullGaussianGrid(pres2, input_as=Matrix))

temperature has to be in K, relative humidity in kg/kg, pressure is actually the logarithm of surface pressure. Maybe we should convert that automatically but currently we don't. So put a log.(...) around one of the arrays.

@milankl
Copy link
Member

milankl commented Dec 13, 2024

  1. For your horizontal bilinear_interpolation, we have the RingGrids module for that with interpolation explained in the documentation. e.g.
using SpeedyWeather
grid24 = rand(FullClenshawGrid, 24)
grid12 = RingGrids.interpolate(FullClenshawGrid, 12, grid)

would interpolate from a 96x47 grid to a 48x23.

  1. sort sorts the values not the coordinates! what you want is to but bottom to top and top to bottom, which is reverting the 2nd dimension, so swapping the left column with the right here (the horizontal dimension is flattened into the first dimension as a vector)
julia> grid = rand(FullClenshawGrid, 2, 3)
24×3, 3-ring FullClenshawArray{Float64, 2, Matrix{Float64}}:
 0.933543   0.301979   0.571856
 0.21249    0.951522   0.893032
 0.67024    0.555407   0.791283
 0.50855    0.802504   0.250079
 0.358565   0.933129   0.222025
 0.833646   0.143867   0.285908
 0.275013   0.982804   0.785195
 0.0837254  0.924352   0.300314
 0.258964   0.857231   0.221561
 0.888978   0.0955061  0.452848
 0.709174   0.59125    0.597493
 0.986523   0.107532   0.60472
 0.471154   0.417332   0.705191
 0.0497819  0.682505   0.144499
 0.984909   0.906613   0.299998
 0.302946   0.077889   0.88048
 0.233419   0.700507   0.782197
 0.846485   0.255638   0.935966
 0.341638   0.355736   0.541445
 0.0896589  0.022498   0.383306
 0.339466   0.853529   0.969174
 0.618453   0.755872   0.610549
 0.155821   0.0678199  0.596454
 0.137743   0.0420721  0.351357

julia> grid[:, end:-1:1]
24×3 Matrix{Float64}:
 0.571856  0.301979   0.933543
 0.893032  0.951522   0.21249
 0.791283  0.555407   0.67024
 0.250079  0.802504   0.50855
 0.222025  0.933129   0.358565
 0.285908  0.143867   0.833646
 0.785195  0.982804   0.275013
 0.300314  0.924352   0.0837254
 0.221561  0.857231   0.258964
 0.452848  0.0955061  0.888978
 0.597493  0.59125    0.709174
 0.60472   0.107532   0.986523
 0.705191  0.417332   0.471154
 0.144499  0.682505   0.0497819
 0.299998  0.906613   0.984909
 0.88048   0.077889   0.302946
 0.782197  0.700507   0.233419
 0.935966  0.255638   0.846485
 0.541445  0.355736   0.341638
 0.383306  0.022498   0.0896589
 0.969174  0.853529   0.339466
 0.610549  0.755872   0.618453
 0.596454  0.0678199  0.155821
 0.351357  0.0420721  0.137743
  1. you can always check what you did by plotting the variables. To get them from spectral to grid just run a period=Day(0) simulation. Then use the visualisation tools as outlined in the documentation, e.g. https://speedyweather.github.io/SpeedyWeather.jl/dev/initial_conditions/

@milankl
Copy link
Member

milankl commented Dec 13, 2024

julia> grid[:, end:-1:1]
24×3 Matrix{Float64}:

Mhhh this escapes the grid, just wrap it back in one, like FullClenshawArray(grid[:, end:-1:1])

@milankl milankl changed the title Inquiry on Using Reanalysis Data as Initial data in SpeedyWeather.jl Reanalysis data as initial conditions Dec 13, 2024
@rocroc2017
Copy link
Contributor Author

  1. For your horizontal bilinear_interpolation, we have the RingGrids module for that with interpolation explained in the documentation. e.g.
using SpeedyWeather
grid24 = rand(FullClenshawGrid, 24)
grid12 = RingGrids.interpolate(FullClenshawGrid, 12, grid)

would interpolate from a 96x47 grid to a 48x23.

  1. sort sorts the values not the coordinates! what you want is to but bottom to top and top to bottom, which is reverting the 2nd dimension, so swapping the left column with the right here (the horizontal dimension is flattened into the first dimension as a vector)
julia> grid = rand(FullClenshawGrid, 2, 3)
24×3, 3-ring FullClenshawArray{Float64, 2, Matrix{Float64}}:
 0.933543   0.301979   0.571856
 0.21249    0.951522   0.893032
 0.67024    0.555407   0.791283
 0.50855    0.802504   0.250079
 0.358565   0.933129   0.222025
 0.833646   0.143867   0.285908
 0.275013   0.982804   0.785195
 0.0837254  0.924352   0.300314
 0.258964   0.857231   0.221561
 0.888978   0.0955061  0.452848
 0.709174   0.59125    0.597493
 0.986523   0.107532   0.60472
 0.471154   0.417332   0.705191
 0.0497819  0.682505   0.144499
 0.984909   0.906613   0.299998
 0.302946   0.077889   0.88048
 0.233419   0.700507   0.782197
 0.846485   0.255638   0.935966
 0.341638   0.355736   0.541445
 0.0896589  0.022498   0.383306
 0.339466   0.853529   0.969174
 0.618453   0.755872   0.610549
 0.155821   0.0678199  0.596454
 0.137743   0.0420721  0.351357

julia> grid[:, end:-1:1]
24×3 Matrix{Float64}:
 0.571856  0.301979   0.933543
 0.893032  0.951522   0.21249
 0.791283  0.555407   0.67024
 0.250079  0.802504   0.50855
 0.222025  0.933129   0.358565
 0.285908  0.143867   0.833646
 0.785195  0.982804   0.275013
 0.300314  0.924352   0.0837254
 0.221561  0.857231   0.258964
 0.452848  0.0955061  0.888978
 0.597493  0.59125    0.709174
 0.60472   0.107532   0.986523
 0.705191  0.417332   0.471154
 0.144499  0.682505   0.0497819
 0.299998  0.906613   0.984909
 0.88048   0.077889   0.302946
 0.782197  0.700507   0.233419
 0.935966  0.255638   0.846485
 0.541445  0.355736   0.341638
 0.383306  0.022498   0.0896589
 0.969174  0.853529   0.339466
 0.610549  0.755872   0.618453
 0.596454  0.0678199  0.155821
 0.351357  0.0420721  0.137743
  1. you can always check what you did by plotting the variables. To get them from spectral to grid just run a period=Day(0) simulation. Then use the visualisation tools as outlined in the documentation, e.g. https://speedyweather.github.io/SpeedyWeather.jl/dev/initial_conditions/

Dear Dr. Klöwer
I followed your suggestion and found that when I run the simulation with period=Day(), the results display normally when Day is between 0 and 6. However, when Day reaches 7, the results seem to become abnormal. I believe the issue with using initial values has been mostly resolved, and I agree with closing the issue. The next step I need to address is how to ensure that the initial values I set can be integrated over a long period

@milankl
Copy link
Member

milankl commented Dec 16, 2024

I believe the issue with using initial values has been mostly resolved, and I agree with closing the issue.

Sorry this was automatically closed because I linked the pull request, we can continue the discussion here, might be better

@milankl milankl reopened this Dec 16, 2024
@milankl
Copy link
Member

milankl commented Dec 16, 2024

From #643

Dear Dr. Klöwer
I plot the simulation.prognostic_variables(temp,humid,pres,vor and div) in layer8
temp in day 0 and day 6
t0
t6
humid in day 0 and day 6
humid0
humid6
pres in day 0 and day6
pres0
pres6
vor in day 0 and day6
vor0
vor6
div in day 0 and day6
div0
div6

I find there are small changes in temp, humid and pres, but there are large changes in vor and div.
However, I am unsure of the next steps to take in order to identify the issue and ensure the model can run stably for long-term integration. Could you please provide some suggestions? Thank you very much.

@rocroc2017
Copy link
Contributor Author

Dear Dr. Klöwer
I use T10 (32x16) and
set!(simulation, u=u_grid, v=v_grid, temp=t_grid, humid=sh_grid)
It can run 1000 days, and the result looks right.
But when I change the pres using
set!(simulation, u=u_grid, v=v_grid, temp=t_grid, humid=sh_grid, pres=pres_grid)
It can only run 6 days.
I plot the pres using default and change to reanalysis data
default in day0
pres0day0
reanalysis data in day0
pres1day0

There seems to be no significant difference in magnitude between the two graphs.
Could you give me some suggestions? Thank you very much

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐞 Something isn't working initial conditions 🚥 Initial conditions of model simulations question ❓ ❔ We got answers!!
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants