Add low-resolution dataset for DYAMOND summer initial conditions.

atmos_dyamond_summer/Project.toml

@@ -1,2 +1,7 @@
 [deps]
 ClimaArtifactsHelper = "6ffa2572-8378-4377-82eb-ea11db28b255"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"

atmos_dyamond_summer/create_artifact.jl

@@ -1,8 +1,276 @@
 using ClimaArtifactsHelper
+using Interpolations
+using NCDatasets
+using DataStructures
+import ClimaParams
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TP
 
 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 = 7000.0
+const P0 = 1e5
+const log_P0 = log(P0)
+const z_min, z_max = 10.0, 80E3
+
+Plvl(z) = P0 * exp(-z / H_EARTH)
+Plvl_inv(P) = -H_EARTH * log(P / P0)
+Plvl_inv2(log_p_by_p0) = -H_EARTH * log_p_by_p0
+
+compute_P(hyam, hybm, surfacepress) = hyam + hybm * surfacepress
+
+param_set = TP.ThermodynamicsParameters(FT)
+const grav = TP.grav(param_set)
+
+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, small=false)
+    ncin = NCDataset(infile, "r")
+    ncout = NCDataset(outfile, "c")
+
+    lonidx = 1:skip:ncin.dim["lon"]
+    latidx = 1:skip:ncin.dim["lat"]
+    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))
+        log_p_by_p0 = log.(P) .- log_P0
+        source_z = Plvl_inv2.(log_p_by_p0)
+
+        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 = interpz_3d(target_z, source_z, ncin["q"][lonidx, latidx, zidx, 1])
+        q[:, :, :] = max.(_q, FT(0))
+
+        # clwc (Specific cloud liquid water content)
+        clwc = defVar(ncout, "clwc", FT, ("lon", "lat", "z",), attrib = ncin["clwc"].attrib)
+        clwc[:, :, :] = interpz_3d(target_z, source_z, ncin["clwc"][lonidx, latidx, zidx, 1])
+
+        # ciwc (Specific cloud ice water content)
+        ciwc = defVar(ncout, "ciwc", FT, ("lon", "lat", "z",), attrib = ncin["ciwc"].attrib)
+        ciwc[:, :, :] = interpz_3d(target_z, source_z, ncin["ciwc"][lonidx, latidx, zidx, 1])
+
+        # crwc (Specific rain water content)
+        crwc = defVar(ncout, "crwc", FT, ("lon", "lat", "z",), attrib = ncin["crwc"].attrib)
+        crwc[:, :, :] = interpz_3d(target_z, source_z, ncin["crwc"][lonidx, latidx, zidx, 1])
+
+        # cswc (Specific snow water content)
+        cswc = defVar(ncout, "cswc", FT, ("lon", "lat", "z",), attrib = ncin["cswc"].attrib)
+        cswc[:, :, :] = interpz_3d(target_z, source_z, ncin["cswc"][lonidx, latidx, zidx, 1])
+
+        # rho (density)
+        phase_partition = TD.PhasePartition.(q)
+
+        rho = defVar(ncout, "rho", FT, ("lon", "lat", "z",), attrib = OrderedDict(
+            "Datatype" => string(FT),
+            "standard_name" => "air_density",
+            "long_name" => "air_density",
+            "units" => "Kg/m**3",
+            ),
+        )
+        rho[:, :, :] = TD.air_density.(param_set, t[:, :, :], P, phase_partition)
+
+        # total energy
+        e_kin = (u .* u .+ v .* v .+ w .* w) .* FT(0.5)
+        e_pot = repeat(reshape(target_z, 1, 1, source_nz), nlon, nlat, 1) .* grav
+
+        e_tot = defVar(ncout, "e_tot", FT, ("lon", "lat", "z",), attrib = OrderedDict(
+            "Datatype" => string(FT),
+            "standard_name" => "total_energy",
+            "long_name" => "total_energy",
+            "units" => "J",
+            ),
+        )
+        e_tot[:, :, :] = TD.total_energy.(param_set, e_kin, e_pot, t[:, :, :], phase_partition)
+
+        # 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")
+@time create_initial_conditions(FILE_PATH, "DYAMOND_SUMMER_ICS_p98deg.nc", skip=7)