From b07ebc67ceb396bb1328fd509f14673c7455b015 Mon Sep 17 00:00:00 2001 From: sriharshakandala Date: Sat, 23 Nov 2024 13:53:00 -0800 Subject: [PATCH] Add low-resolution dataset for DYAMOND summer initial conditions. --- atmos_dyamond_summer/OutputArtifacts.toml | 6 + atmos_dyamond_summer/Project.toml | 3 + atmos_dyamond_summer/README.md | 8 +- atmos_dyamond_summer/create_artifact.jl | 249 +++++++++++++++++++++- 4 files changed, 264 insertions(+), 2 deletions(-) diff --git a/atmos_dyamond_summer/OutputArtifacts.toml b/atmos_dyamond_summer/OutputArtifacts.toml index 1e6a8ff..3d057dd 100644 --- a/atmos_dyamond_summer/OutputArtifacts.toml +++ b/atmos_dyamond_summer/OutputArtifacts.toml @@ -1,2 +1,8 @@ [DYAMOND_summer_initial_conditions] git-tree-sha1 = "931ce29dded2271fddd537a2492f3aef260624fd" +[DYAMOND_SUMMER_ICS_p98deg] +git-tree-sha1 = "e6827d642305864aa2a142789805150ec04c9fdb" + + [[DYAMOND_SUMMER_ICS_p98deg.download]] + sha256 = "1f933b995be6548080cbc4e73fc5c9f0fe82e9f6ec566ea69b2a414c763c7529" + url = "https://caltech.box.com/shared/static/1nqueb39ro33b2ps7zm2df9hpbkps6wy.gz" diff --git a/atmos_dyamond_summer/Project.toml b/atmos_dyamond_summer/Project.toml index 8e4179f..95bdcbc 100644 --- a/atmos_dyamond_summer/Project.toml +++ b/atmos_dyamond_summer/Project.toml @@ -1,2 +1,5 @@ [deps] ClimaArtifactsHelper = "6ffa2572-8378-4377-82eb-ea11db28b255" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" diff --git a/atmos_dyamond_summer/README.md b/atmos_dyamond_summer/README.md index 1371f2d..9bade36 100644 --- a/atmos_dyamond_summer/README.md +++ b/atmos_dyamond_summer/README.md @@ -2,11 +2,17 @@ ## Overview The DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains (DYAMOND) protocol provides a framework for the intercomparison of global storm-resolving models. -This artifact contains the initial conditions needed for the DYAMOND summer simulations, which simulate 40 days and nights, begining on 1 August 2016. +This artifact contains the initial conditions needed for the DYAMOND summer simulations, which simulate 40 days and nights, begining on 1 August 2016. + +The downloaded artifact has a resolution of about 0.14 degrees. A low-resolution version of this data, with a resolution of about 1 degree (0.98 degrees) and a custom vertical grid that is uniform across all columns, is generated. The following procedure is followed for generating the low-resolution data: + +- Since the horizontal grid is uniformly spaced, a low-resolution horizontal grid is generated by skipping some horizontal points. +- The input file uses hybrid sigma pressure coordinates in the vertical direction. First, we compute the pressure at the vertical levels using $P = ap(k) + P_S b(k)$. In the file, this is $P_S$ = `exp(lnsp)`, b = `hybm` and a = `hyam`. Next, we assume a hydrostatic profile and convert pressure to altitude over the mean sea level using $P = P0*exp(-z / H)$ with scale height $H$. We assume $P^* = 1e5$ (Pa) and $H = 7000$(m). Now, we select our target z, from 10m to 80km, and interpolate all source data onto this grid. We use "Flat" boundary conditions for those points that are outside the range of definition. ## File summmary - `ifs_oper_T1279_2016080100.nc`: Restart file for DYAMOND summer simulation (~18 GB) +- `DYAMOND_SUMMER_ICS_p98deg.nc`: Initial conditions for DYAMOND summer simulation at 0.98 degree resolution (~436 MB) ## References diff --git a/atmos_dyamond_summer/create_artifact.jl b/atmos_dyamond_summer/create_artifact.jl index 5932557..0eab918 100644 --- a/atmos_dyamond_summer/create_artifact.jl +++ b/atmos_dyamond_summer/create_artifact.jl @@ -1,8 +1,255 @@ using ClimaArtifactsHelper +using Interpolations +using NCDatasets +using DataStructures const FILE_URL = "https://swift.dkrz.de/v1/dkrz_ab6243f85fe24767bb1508712d1eb504/SAPPHIRE/DYAMOND/ifs_oper_T1279_2016080100.nc" const FILE_PATH = "ifs_oper_T1279_2016080100.nc" artifact_name = "DYAMOND_summer_initial_conditions" - create_artifact_guided_one_file(FILE_PATH; artifact_name = artifact_name, file_url = FILE_URL) + +const FT = Float32 +const H_EARTH = FT(7000.0) +const P0 = FT(1e5) +const z_min, z_max = 10.0, 80E3 + +Plvl(z) = P0 * exp(-z / H_EARTH) +Plvl_inv(P) = -H_EARTH * log(P / P0) + +compute_P(hyam, hybm, surfacepress) = hyam + hybm * surfacepress + +function interpz_3d(target_z, z3d, var3d) + nx, ny, nz = size(z3d) + target_var3d = similar(var3d, nx, ny, length(target_z)) + @inbounds begin + for yi in 1:ny, xi in 1:nx + source_z = view(z3d, xi, yi, :) + itp = extrapolate( + interpolate( + (source_z,), + (view(var3d, xi, yi, :)), + Gridded(Linear()), + ), + Flat(), + ) + target_var3d[xi, yi, :] .= itp.(target_z) + end + end + return target_var3d +end + +function create_initial_conditions(infile, outfile; skip = 1) + ncin = NCDataset(infile, "r") + ncout = NCDataset(outfile, "c", attrib = copy(ncin.attrib)) + ncout.attrib["history"] *= "; Modified by CliMA (see atmos_dyamond_summer in ClimaArtifacts)" + + lonidx = 1:skip:ncin.dim["lon"] + latidx = 1:skip:ncin.dim["lat"] + # In this dataset vertical levels are ordered in the direction of + # increasing pressure (decreasing elevation). This is being reordered + # to stay consistent with our code which orders vertical levels in the + # direction of increasing elevation + zidx = ncin.dim["lev"]:-1:1 + + nlon = length(lonidx) + nlat = length(latidx) + source_nz = length(zidx) + + @inbounds begin + # get source z + sp = repeat(exp.(ncin["lnsp"][lonidx, latidx, 1, 1]), 1, 1, source_nz) # surface pressure + hyam = repeat(reshape(ncin["hyam"][zidx], 1, 1, source_nz), nlon, nlat, 1) + hybm = repeat(reshape(ncin["hybm"][zidx], 1, 1, source_nz), nlon, nlat, 1) + + P = FT.(compute_P.(hyam, hybm, sp)) + source_z = Plvl_inv.(P) + + nz = source_nz #length(zidx) + + # create a target z grid with nz points + target_z = FT.(Plvl_inv.(range(Plvl(z_min), Plvl(z_max), nz))) + # defining dimensions + defDim(ncout, "lon", nlon) + defDim(ncout, "lat", nlat) + defDim(ncout, "z", nz) + # longitude + lon = defVar(ncout, "lon", FT, ("lon",), attrib = ncin["lon"].attrib) + lon[:] = Array{FT}(ncin["lon"][lonidx]) + + # latitude + lat = defVar(ncout, "lat", FT, ("lat",), attrib = ncin["lat"].attrib) + lat[:] = Array{FT}(ncin["lat"][latidx]) + + z = defVar(ncout, "z", FT, ( "z",), attrib = OrderedDict( + "Datatype" => string(FT), + "standard_name" => "altitude", + "long_name" => "altitude", + "units" => "m", + ), + ) + z[:] = target_z + + # p (pressure) + p = defVar(ncout, "p", FT, ("lon", "lat", "z",), attrib = OrderedDict( + "Datatype" => string(FT), + "standard_name" => "pressure", + "long_name" => "pressure", + "units" => "Pa", + ), + ) + p[:, :, :] = P + + + # u (eastward_wind) + u = defVar(ncout, "u", FT, ("lon", "lat", "z",), attrib = ncin["u"].attrib) + u[:, :, :] = interpz_3d(target_z, source_z, ncin["u"][lonidx, latidx, zidx, 1]) + + # v (northward_wind) + v = defVar(ncout, "v", FT, ("lon", "lat", "z",), attrib = ncin["v"].attrib) + v[:, :, :] = interpz_3d(target_z, source_z, ncin["v"][lonidx, latidx, zidx, 1]) + + # w (vertical velocity) + w = defVar(ncout, "w", FT, ("lon", "lat", "z",), attrib = ncin["w"].attrib) + w[:, :, :] = interpz_3d(target_z, source_z, ncin["w"][lonidx, latidx, zidx, 1]) + + # t (air_temperature) + t = defVar(ncout, "t", FT, ("lon", "lat", "z",), attrib = ncin["t"].attrib) + t[:, :, :] = interpz_3d(target_z, source_z, ncin["t"][lonidx, latidx, zidx, 1]) + + # q (specific_humidity) + q = defVar(ncout, "q", FT, ("lon", "lat", "z",), attrib = ncin["q"].attrib) + q[:, :, :] = max.(interpz_3d(target_z, source_z, ncin["q"][lonidx, latidx, zidx, 1]), FT(0)) + + # clwc (Specific cloud liquid water content) + clwc = defVar(ncout, "clwc", FT, ("lon", "lat", "z",), attrib = ncin["clwc"].attrib) + clwc[:, :, :] = max.(interpz_3d(target_z, source_z, ncin["clwc"][lonidx, latidx, zidx, 1]), FT(0)) + + # ciwc (Specific cloud ice water content) + ciwc = defVar(ncout, "ciwc", FT, ("lon", "lat", "z",), attrib = ncin["ciwc"].attrib) + ciwc[:, :, :] = max.(interpz_3d(target_z, source_z, ncin["ciwc"][lonidx, latidx, zidx, 1]), FT(0)) + + # crwc (Specific rain water content) + crwc = defVar(ncout, "crwc", FT, ("lon", "lat", "z",), attrib = ncin["crwc"].attrib) + crwc[:, :, :] = max.(interpz_3d(target_z, source_z, ncin["crwc"][lonidx, latidx, zidx, 1]), FT(0)) + + # cswc (Specific snow water content) + cswc = defVar(ncout, "cswc", FT, ("lon", "lat", "z",), attrib = ncin["cswc"].attrib) + cswc[:, :, :] = max.(interpz_3d(target_z, source_z, ncin["cswc"][lonidx, latidx, zidx, 1]), FT(0)) + + # tsn (temperature_in_surface_snow) + tsn = defVar(ncout, "tsn", FT, ("lon", "lat",), attrib = ncin["tsn"].attrib) + tsn[:, :] = ncin["tsn"][lonidx, latidx, 1] + + # skt (skin temperature) + skt = defVar(ncout, "skt", FT, ("lon", "lat",), attrib = ncin["skt"].attrib) + skt[:, :] = ncin["skt"][lonidx, latidx, 1] + + # sst (sea surface temperature) + sst = defVar(ncout, "sst", FT, ("lon", "lat",), attrib = ncin["sst"].attrib) + sst[:, :] = ncin["sst"][lonidx, latidx, 1] + + # stl1 (soil temperature level 1) + stl1 = defVar(ncout, "stl1", FT, ("lon", "lat",), attrib = ncin["stl1"].attrib) + stl1[:, :] = ncin["stl1"][lonidx, latidx, 1] + + # stl2 (soil temperature level 2) + stl2 = defVar(ncout, "stl2", FT, ("lon", "lat",), attrib = ncin["stl2"].attrib) + stl2[:, :] = ncin["stl2"][lonidx, latidx, 1] + + # stl3 (soil temperature level 3) + stl3 = defVar(ncout, "stl3", FT, ("lon", "lat",), attrib = ncin["stl3"].attrib) + stl3[:, :] = ncin["stl3"][lonidx, latidx, 1] + + # stl4 (soil temperature level 4) + stl4 = defVar(ncout, "stl4", FT, ("lon", "lat",), attrib = ncin["stl4"].attrib) + stl4[:, :] = ncin["stl4"][lonidx, latidx, 1] + + # sd (lwe_thickness_of_surface_snow_amount) + sd = defVar(ncout, "sd", FT, ("lon", "lat",), attrib = ncin["sd"].attrib) + sd[:, :] = ncin["sd"][lonidx, latidx, 1] + + # rsn (snow density) + rsn = defVar(ncout, "rsn", FT, ("lon", "lat",), attrib = ncin["rsn"].attrib) + rsn[:, :] = ncin["rsn"][lonidx, latidx, 1] + + # asn (snow albedo) + asn = defVar(ncout, "asn", FT, ("lon", "lat",), attrib = ncin["asn"].attrib) + asn[:, :] = ncin["asn"][lonidx, latidx, 1] + + # src (skin reservoir content) + src = defVar(ncout, "src", FT, ("lon", "lat",), attrib = ncin["src"].attrib) + src[:, :] = ncin["src"][lonidx, latidx, 1] + + # ci (Sec-ice cover) + ci = defVar(ncout, "ci", FT, ("lon", "lat",), attrib = ncin["ci"].attrib) + ci[:, :] = ncin["ci"][lonidx, latidx, 1] + + # swvl1 (volumetric soil water layer 1) + swvl1 = defVar(ncout, "swvl1", FT, ("lon", "lat",), attrib = ncin["swvl1"].attrib) + swvl1[:, :] = ncin["swvl1"][lonidx, latidx, 1] + + # swvl2 (volumetric soil water layer 2) + swvl2 = defVar(ncout, "swvl2", FT, ("lon", "lat",), attrib = ncin["swvl2"].attrib) + swvl2[:, :] = ncin["swvl2"][lonidx, latidx, 1] + + # swvl3 (volumetric soil water layer 3) + swvl3 = defVar(ncout, "swvl3", FT, ("lon", "lat",), attrib = ncin["swvl3"].attrib) + swvl3[:, :] = ncin["swvl3"][lonidx, latidx, 1] + + # swvl4 (volumetric soil water layer 4) + swvl4 = defVar(ncout, "swvl4", FT, ("lon", "lat",), attrib = ncin["swvl4"].attrib) + swvl4[:, :] = ncin["swvl4"][lonidx, latidx, 1] + + # slt (soil type) + slt = defVar(ncout, "slt", FT, ("lon", "lat",), attrib = ncin["slt"].attrib) + slt[:, :] = ncin["slt"][lonidx, latidx, 1] + + # lsm (land-sea mask) + lsm = defVar(ncout, "lsm", FT, ("lon", "lat",), attrib = ncin["lsm"].attrib) + lsm[:, :] = ncin["lsm"][lonidx, latidx, 1] + + # sr (surface roughness length) + sr = defVar(ncout, "sr", FT, ("lon", "lat",), attrib = ncin["sr"].attrib) + sr[:, :] = ncin["sr"][lonidx, latidx, 1] + + # cvl (low vegetation cover) + cvl = defVar(ncout, "cvl", FT, ("lon", "lat",), attrib = ncin["cvl"].attrib) + cvl[:, :] = ncin["cvl"][lonidx, latidx, 1] + + # cvh (high vegetation cover) + cvh = defVar(ncout, "cvh", FT, ("lon", "lat",), attrib = ncin["cvh"].attrib) + cvh[:, :] = ncin["cvh"][lonidx, latidx, 1] + + # sdor (standard deviation or orography) + sdor = defVar(ncout, "sdor", FT, ("lon", "lat",), attrib = ncin["sdor"].attrib) + sdor[:, :] = ncin["sdor"][lonidx, latidx, 1] + + # isor (anisotropy of sub-grid scale orography) + isor = defVar(ncout, "isor", FT, ("lon", "lat",), attrib = ncin["isor"].attrib) + isor[:, :] = ncin["isor"][lonidx, latidx, 1] + + # anor (angle of sub-grid scale orography) + anor = defVar(ncout, "anor", FT, ("lon", "lat",), attrib = ncin["anor"].attrib) + anor[:, :] = ncin["anor"][lonidx, latidx, 1] + + # slor (slope of sub-grid scale orography) + slor = defVar(ncout, "slor", FT, ("lon", "lat",), attrib = ncin["slor"].attrib) + slor[:, :] = ncin["slor"][lonidx, latidx, 1] + end + close(ncout) + return nothing +end + +println("creating initial conditions for 0.98⁰ resolution") + +artifact_dir_p98deg = artifact_name_p98deg = "DYAMOND_SUMMER_ICS_p98deg" + +if !isdir(artifact_dir_p98deg) + mkdir(artifact_dir_p98deg) +end +file_path_p98deg = joinpath(@__DIR__, artifact_dir_p98deg, artifact_name_p98deg * ".nc") + +@time create_initial_conditions(FILE_PATH, file_path_p98deg, skip=7) + +create_artifact_guided(artifact_dir_p98deg, artifact_name = artifact_name_p98deg, append = true)