Adding an ECCO prescribed atmosphere

src/DataWrangling/ECCO/ECCO.jl

@@ -42,7 +42,8 @@ import ClimaOcean.DataWrangling:
     is_three_dimensional,
     inpainted_metadata_path,
     reversed_vertical_axis,
-    default_mask_value
+    default_mask_value,
+    reversed_sign
 
 download_ECCO_cache::String = ""
 function __init__()
@@ -77,6 +78,8 @@ Base.size(::ECCO4Monthly, variable) = (720,  360, 50)
 temperature_units(::SomeECCODataset) = Celsius()
 default_mask_value(::ECCO4Monthly) = 0
 reversed_vertical_axis(::SomeECCODataset) = true
+reversed_sign(::SomeECCODataset, ::Val{:downwelling_longwave}) = true
+reversed_sign(::SomeECCODataset, ::Val{:downwelling_shortwave}) = true
 
 const ECCO2_url = "https://ecco.jpl.nasa.gov/drive/files/ECCO2/cube92_latlon_quart_90S90N/"
 const ECCO4_url = "https://ecco.jpl.nasa.gov/drive/files/Version4/Release4/interp_monthly/"
@@ -157,35 +160,47 @@ ECCO4_dataset_variable_names = Dict(
     :latent_heat_flux      => "EXFhl",
     :net_longwave          => "EXFlwnet",
     :downwelling_shortwave => "oceQsw",
-    :downwelling_longwave  => "EXFlwdn",
+    :downwelling_longwave  => "EXFlwdn",    
+    :air_temperature       => "EXFatemp",
+    :air_specific_humidity => "EXFaqh",
+    :sea_level_pressure    => "EXFpress",
+    :eastward_wind         => "EXFewind",
+    :northward_wind        => "EXFnwind",
+    :rain_freshwater_flux  => "EXFpreci",
 )
 
 ECCO2_dataset_variable_names = Dict(
     :temperature           => "THETA",
     :salinity              => "SALT",
-    :u_velocity            => "UVEL",
-    :v_velocity            => "VVEL",
+    :u_velocity            => "EVEL",
+    :v_velocity            => "NVEL",
     :free_surface          => "SSH",
     :sea_ice_thickness     => "SIheff",
     :sea_ice_concentration => "SIarea",
-    :net_heat_flux         => "oceQnet"
+    :net_heat_flux         => "oceQnet",
 )
 
 ECCO_location = Dict(
     :temperature           => (Center, Center, Center),
     :salinity              => (Center, Center, Center),
+    :u_velocity            => (Face,   Center, Center),
+    :v_velocity            => (Center, Face,   Center),
     :free_surface          => (Center, Center, Nothing),
     :sea_ice_thickness     => (Center, Center, Nothing),
     :sea_ice_concentration => (Center, Center, Nothing),
     :net_heat_flux         => (Center, Center, Nothing),
-    :u_velocity            => (Face,   Center, Center),
-    :v_velocity            => (Center, Face,   Center),
     :sensible_heat_flux    => (Center, Center, Nothing),
     :latent_heat_flux      => (Center, Center, Nothing),
     :net_longwave          => (Center, Center, Nothing),
-    :downwelling_shortwave => (Center, Center, Nothing),
     :downwelling_longwave  => (Center, Center, Nothing),
-)
+    :downwelling_shortwave => (Center, Center, Nothing),
+    :air_temperature       => (Center, Center, Nothing),
+    :air_specific_humidity => (Center, Center, Nothing),
+    :sea_level_pressure    => (Center, Center, Nothing),
+    :eastward_wind         => (Center, Center, Nothing),
+    :northward_wind        => (Center, Center, Nothing),
+    :rain_freshwater_flux  => (Center, Center, Nothing),
+)    
 
 const ECCOMetadata{D} = Metadata{<:SomeECCODataset, D}
 const ECCOMetadatum   = Metadatum{<:SomeECCODataset}
@@ -297,4 +312,6 @@ end
 
 inpainted_metadata_path(metadata::ECCOMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
 
+include("ECCO_atmosphere.jl")
+
 end # Module

src/DataWrangling/ECCO/ECCO_atmosphere.jl

@@ -0,0 +1,75 @@
+using ClimaOcean.DataWrangling: DatasetBackend
+using Oceananigans.OutputReaders
+using ClimaOcean.OceanSeaIceModels: TwoBandDownwellingRadiation
+
+"""
+    ECCOPrescribedAtmosphere([architecture = CPU(), FT = Float32];
+                              dataset = ECCO4Montly(),
+                              start_date = first_date(dataset, :temperature),
+                              end_date = last_date(dataset, :temperature),
+                              backend = DatasetBackend(10),
+                              time_indexing = Cyclical(),
+                              surface_layer_height = 10,  # meters
+                              include_rivers_and_icebergs = false,
+                              other_kw...)
+
+Return a [`PrescribedAtmosphere`](@ref) representing JRA55 reanalysis data.
+The atmospheric data will be held in `JRA55FieldTimeSeries` objects containing.
+For a detailed description of the keyword arguments, see the [`JRA55FieldTimeSeries`](@ref) constructor.
+"""
+function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32;
+                                  dataset = ECCO4Monthly(),
+                                  start_date = first_date(dataset, :air_temperature),
+                                  end_date = last_date(dataset, :air_temperature),
+                                  time_indexing = Cyclical(),
+                                  time_indices_in_memory = 10,
+                                  surface_layer_height = 2,  # meters
+                                  other_kw...)
+
+    ua_meta = Metadata(:eastward_wind;         dataset, start_date, end_date)    
+    va_meta = Metadata(:northward_wind;        dataset, start_date, end_date)    
+    Ta_meta = Metadata(:air_temperature;       dataset, start_date, end_date)    
+    qa_meta = Metadata(:air_specific_humidity; dataset, start_date, end_date)    
+    pa_meta = Metadata(:sea_level_pressure;    dataset, start_date, end_date)    
+    Ql_meta = Metadata(:downwelling_longwave;  dataset, start_date, end_date)
+    Qs_meta = Metadata(:downwelling_shortwave; dataset, start_date, end_date)
+    Fr_meta = Metadata(:rain_freshwater_flux;  dataset, start_date, end_date)
+
+    kw = (; time_indices_in_memory, time_indexing)
+    kw = merge(kw, other_kw)
+
+    ua = FieldTimeSeries(ua_meta, architecture; kw...)
+    va = FieldTimeSeries(va_meta, architecture; kw...)
+    Ta = FieldTimeSeries(Ta_meta, architecture; kw...)
+    qa = FieldTimeSeries(qa_meta, architecture; kw...)
+    pa = FieldTimeSeries(pa_meta, architecture; kw...)
+    Ql = FieldTimeSeries(Ql_meta, architecture; kw...)
+    Qs = FieldTimeSeries(Qs_meta, architecture; kw...)
+    Fr = FieldTimeSeries(Fr_meta, architecture; kw...)
+
+    auxiliary_freshwater_flux = nothing
+    freshwater_flux = (; rain = Fr)
+
+    times = ua.times
+    grid  = ua.grid
+
+    velocities = (u = ua, v = va)
+    tracers = (T = Ta, q = qa)
+    pressure = pa
+
+    downwelling_radiation = TwoBandDownwellingRadiation(shortwave=Qs, longwave=Ql)
+
+    FT = eltype(ua)
+    surface_layer_height = convert(FT, surface_layer_height)
+
+    atmosphere = PrescribedAtmosphere(grid, times;
+                                      velocities,
+                                      freshwater_flux,
+                                      auxiliary_freshwater_flux,
+                                      tracers,
+                                      downwelling_radiation,
+                                      surface_layer_height,
+                                      pressure)
+
+    return atmosphere
+end

src/DataWrangling/metadata_field.jl

@@ -150,6 +150,8 @@ end
 struct ShiftSouth end
 struct AverageNorthSouth end
 
+reversed_sign(dataset, val_name) = false
+
 @inline mangle(i, j, data, ::Nothing) = @inbounds data[i, j]
 @inline mangle(i, j, data, ::ShiftSouth) = @inbounds data[i, j-1]
 @inline mangle(i, j, data, ::AverageNorthSouth) = @inbounds (data[i, j+1] + data[i, j]) / 2
@@ -158,6 +160,8 @@ struct AverageNorthSouth end
 @inline mangle(i, j, k, data, ::ShiftSouth) = @inbounds data[i, j-1, k]
 @inline mangle(i, j, k, data, ::AverageNorthSouth) = @inbounds (data[i, j+1, k] + data[i, j, k]) / 2
 
+@inline maybe_reverse_sign(datum, reverse::Bool) = ifelse(reverse, - datum, datum)
+
 function set_metadata_field!(field, data, metadatum)
     grid = field.grid
     arch = architecture(grid)
@@ -177,6 +181,8 @@ function set_metadata_field!(field, data, metadatum)
         nothing
     end
 
+    reverse_sign = reversed_sign(metadatum.dataset, Val(metadatum.name))
+
     if ndims(data) == 2
         _kernel = _set_2d_metadata_field!
         spec = :xy
@@ -186,29 +192,29 @@ function set_metadata_field!(field, data, metadatum)
     end
 
     data = on_architecture(arch, data)
-    Oceananigans.Utils.launch!(arch, grid, spec, _kernel, field, data, mangling, temp_units)
+    Oceananigans.Utils.launch!(arch, grid, spec, _kernel, field, data, mangling, reverse_sign, temp_units)
 
     return nothing
 end
 
-@kernel function _set_2d_metadata_field!(field, data, mangling, temp_units)
+@kernel function _set_2d_metadata_field!(field, data, mangling, reverse_sign, temp_units)
     i, j = @index(Global, NTuple)
     d = mangle(i, j, data, mangling)
     d = convert_temperature(d, temp_units)
+    d = maybe_reverse_sign(d, reverse_sign)
     @inbounds field[i, j, 1] = d
 end
 
 @inline nan_convert_missing(FT, ::Missing) = convert(FT, NaN)
 @inline nan_convert_missing(FT, d::Number) = convert(FT, d)
 
-@kernel function _set_3d_metadata_field!(field, data, mangling, temp_units)
+@kernel function _set_3d_metadata_field!(field, data, mangling, reverse_sign, temp_units)
     i, j, k = @index(Global, NTuple)
-    d = mangle(i, j, k, data, mangling)
-
     FT = eltype(field)
+    d = mangle(i, j, k, data, mangling)
     d = nan_convert_missing(FT, d)
-
     d = convert_temperature(d, temp_units)
+    d = maybe_reverse_sign(d, reverse_sign)
     @inbounds field[i, j, k] = d
 end
 

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -159,6 +159,7 @@ struct DegreesCelsius end
 struct DegreesKelvin end
 
 const celsius_to_kelvin = 273.15
+
 @inline convert_to_kelvin(::DegreesCelsius, T::FT) where FT = T + convert(FT, celsius_to_kelvin)
 @inline convert_to_kelvin(::DegreesKelvin, T) = T
 

src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl

@@ -222,8 +222,7 @@ end
     interp_atmos_time_series(values(ΣJ), args...)
 
 @inline interp_atmos_time_series(ΣJ::Tuple{<:Any}, args...) =
-    interp_atmos_time_series(ΣJ[1], args...) +
-    interp_atmos_time_series(ΣJ[2], args...)
+    interp_atmos_time_series(ΣJ[1], args...) 
 
 @inline interp_atmos_time_series(ΣJ::Tuple{<:Any, <:Any}, args...) =
     interp_atmos_time_series(ΣJ[1], args...) +