Add low-resolution dataset for DYAMOND summer initial conditions.

sriharshakandala · sriharshakandala · commit c31beb754581 · 2024-11-25T14:23:43.000-08:00
diff --git 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 = "5ef4bb950eceb86430464c8637ede7f075a1449f"
+
+    [[DYAMOND_SUMMER_ICS_p98deg.download]]
+    sha256 = "da552c32c4e66a053992f3b2c96ebaeb78c1aded1818b4e7c2deea77d521dda1"
+    url = "https://caltech.box.com/shared/static/ad8d774iwaezctwud1g1lqbp9irbjj38.gz"
diff --git 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
@@ -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
@@ -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 = FT(0), FT(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 = reverse(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)
