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

Interpolation on pressure levels #268

Open
milankl opened this issue Mar 23, 2023 · 8 comments
Open

Interpolation on pressure levels #268

milankl opened this issue Mar 23, 2023 · 8 comments
Labels
output 📤 NetCDF or JLD2 output, their writers, formats, etc vertical ⬆️ Affecting the vertical dimension

Comments

@milankl
Copy link
Member

milankl commented Mar 23, 2023

We currently output on $\sigma = p/p_s$ levels, however at some point we may want to interpolate onto pressure levels, the following function outlines how this can be done

function interpolate_pressure_levels(A::AbstractArray{T,3},pₛ::AbstractMatrix::AbstractVector,p::Real) where T
    nx,ny,nlev = size(A)
    B = zeros(T,nx,ny)
    
    for i in 1:nx
        for j in 1:ny
            σij = p/pₛ[i,j]            
            if σij < σ[1]
                B[i,j] = A[i,j,1]
            elseif σij >= σ[end]
                B[i,j] = A[i,j,end]
            else
                k = 1
                while (σ[k] < σij) k += 1 end
                Δσ = σ[k+1] - σ[k]       # level thickness
                w1 = (σ[k+1] - σij)/Δσ   # weight for k level
                w2 = (σij - σ[k])/Δσ     # weight for k+1 level
                B[i,j] = w1*A[i,j,k] + w2*A[i,j,k+1]
            end
        end
    end
    
    return B
end 
@milankl milankl added output 📤 NetCDF or JLD2 output, their writers, formats, etc vertical ⬆️ Affecting the vertical dimension labels Mar 23, 2023
@milankl milankl added this to the v0.6 milestone Mar 23, 2023
@milankl
Copy link
Member Author

milankl commented Mar 31, 2023

That's what Fortran SPEEDY does, probably makes sense regardless of the vertical levels?

image

@milankl milankl modified the milestones: v0.6, v0.7 Aug 30, 2023
@milankl milankl removed this from the v0.7 milestone Mar 19, 2024
@weathermanbarnes
Copy link

Hi SpeedyWeather developers. I have just stumbled across this repo and it looks really great! I am new to Julia, but will definitely be keen to play around with this more! Interpolation and output to pressure levels would be really great. I would also be interested in output to isentropic levels and an option for a potential vorticity variable.

@milankl
Copy link
Member Author

milankl commented Dec 4, 2024

I think there are important considerations for online vs offline interpolation on pressure levels that have pros and cons

  • Interpolate to pressure levels on the fly: For 3D variable the NetCDFOutput has a full 3D variable on model levels on a full grid, say FullGaussianGrid which data from the internal grid, say OctahedralGaussianGrid is first interpolated onto before writing to disk. We could then interpolate in the vertical onto another 3D variable with pressure levels before writing to disk. This would add some compute but would reduce the size of the netCDF file if the number of model levels is higher than the number of pressure levels to output onto.

  • Interpolate at post-processing: The current output on model levels could continue during the simulation but one could interpolate onto pressure levels at the finalize! step. We currently do this to convert accumulated precipitation to precipitation rate, simply be reading in the netCDF variable again, computing a difference and adding a new variable to that same netCDF file. We could do the same for interpolate onto pressure levels and this could also just be implemented as a callback where all interesting stuff happens at finalize!. This would give us the option to "add" pressure levels next to model levels or then also delete the model levels if we don't want to. The extra compute is then only added at the end of a simulation but no storage is saved during the simulation and only at the end if the model level data is deleted.

@milankl
Copy link
Member Author

milankl commented Dec 4, 2024

I currently reckon that post-processing is easier as one could just define a InterpolateToPressureLayers <: SpeedyWeather.AbstractCallback (I believe "layer" is somewhat less ambiguous than "levels" even though that's more commonly used) with a finalize! method similar to

function output!(
output::NetCDFOutput,
variable::AbstractRainRateOutputVariable,
acc_variable::AbstractOutputVariable,
)
# use .var to prevent Union{Missing, Float32} that NCDatasets uses
accumulated = output.netcdf_file[acc_variable.name].var[:, :, :]
# rate is defined as average precip since last output step, so first step is 0
# convert from accumulated [m] to [mm/hr] rain rate over output time step (e.g. 6hours)
s = Hour(1)/output.output_dt
nx, ny = size(accumulated)
rate = cat(zeros(eltype(accumulated), nx, ny), diff(accumulated, dims=3), dims=3)
rate .*= s
# DEFINE NEW NETCDF VARIABLE AND WRITE
define_variable!(output.netcdf_file, variable, eltype(rate))
output.netcdf_file[variable.name][:, :, :] = rate
return nothing
end

and

# at finalize step postprocess the convective precipitation to get the rate
finalize!(output::NetCDFOutput, variable::LargeScalePrecipitationOutput, args...) = output!(output, variable.rate, variable)

but implented as a Callback as that the interface would simply be

add!(model, InterpolateToPressureLayers(spectral_grid, layers=(30, 100, 200, 300, 500, 700, 850, 925))

I guess temperature is extrapolated beyond model layers using potential temperature, but what about humidity? Constant relative humidity of that potential temperature? What about vorticity, divergence, u, v?

@milankl
Copy link
Member Author

milankl commented Dec 4, 2024

Feel free to give this a go and post any progress/problems here? Happy to provide guidance. I've never interpolated to isentropic coordinates, but I guess that would work the same way. Which definition of potential vorticity would you like to output? The latter here from Wikipedia I guess?

image

@weathermanbarnes
Copy link

Thanks for your reply @milankl! I agree that the interpolation as a postprocessing step is probably easiest. I would however suggest that the developers consider a switch to output to either model levels or user-defined pressure levels or isentropic levels. Although it would use some compute, it would save a lot of storage for experiments which require high vertical resolution and which only require particular pressure levels as output. Both options would work and to have both options would be amazing!

In terms of PV, equation 13 is the one you want. Equation 14 is specific to isentropic PV. So I think the procedure would be to calculate it generally (13) and then interpolate to which ever coordinate system you specify.

@milankl
Copy link
Member Author

milankl commented Dec 5, 2024

I would however suggest that the developers consider a switch to output to either model levels or user-defined pressure levels or isentropic levels. Although it would use some compute, it would save a lot of storage for experiments which require high vertical resolution and which only require particular pressure levels as output. Both options would work and to have both options would be amazing!

By developers you mean me? Because it sounds like you're volunteering =) You know there's somewhat a philosophy in Julia around users being developers and developers being users! You find everything in here https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/main/src/output/netcdf_output.jl (most of it is defining how to output various variables and creating a modular interface for it ...

I guess we'd need

  • NetCDFOutput with additional field pressure_levels::Vector{Float64} = [] (default no pressure levels = model levels)
  • probably also const grid3D_pressure_levels::Grid3D to store those before writing to disk
  • on_pressure_levels(output::NetCDFOutput) = length(output.pressure_levels) > 0 as check
  • use that check to switch between defVar(dataset, "layer", σ, ... and the pressure level equivalent
  • use that check to include on_pressure_levels(output) && interpolate_on_pressure_levels!(grid3D_pressure_levels, grid3D, ... for every 3D variable
  • write that interpolate_on_pressure_levels! function as outlined above

that's it! how hard can it be? 😝

@weathermanbarnes
Copy link

weathermanbarnes commented Dec 5, 2024

I hope to get stuck into the code and if I can, I will definitely add to it! Congrats on the work you have done so far!! I look forward to playing around with the code more!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
output 📤 NetCDF or JLD2 output, their writers, formats, etc vertical ⬆️ Affecting the vertical dimension
Projects
None yet
Development

No branches or pull requests

2 participants