some corrections

Author: simone-silvestri

src/DataWrangling/JRA55/JRA55_field_time_series.jl

@@ -189,6 +189,7 @@ function set!(fts::JRA55NetCDFFTSMultipleYears, backend=fts.backend)
             LX, LY, LZ = location(fts)
             i₁, i₂, j₁, j₂, TX = compute_bounding_indices(nothing, nothing, fts.grid, LX, LY, λc, φc)
 
+
             if issorted(nn)
                 data = ds[name][i₁:i₂, j₁:j₂, nn]
             else
@@ -457,7 +458,7 @@ function JRA55FieldTimeSeries(metadata::JRA55Metadata, architecture=CPU(), FT=Fl
         ds = Dataset(filepath)
         data = ds[shortname][i₁:i₂, j₁:j₂, time_indices_in_memory]
         close(ds)
-        
+
         copyto!(interior(fts, :, :, 1, :), data)
         fill_halo_regions!(fts)
 

src/InitialConditions/InitialConditions.jl

@@ -18,7 +18,7 @@ using JLD2
 # Implementation of 3-dimensional regridding
 # TODO: move all the following to Oceananigans!
 
-using Oceananigans.Fields: interpolate!
+using Oceananigans.Fields: regrid!, interpolate!
 using Oceananigans.Grids: cpu_face_constructor_x,
                           cpu_face_constructor_y,
                           cpu_face_constructor_z,
@@ -31,6 +31,44 @@ construct_grid(::Type{<:RectilinearGrid}, arch, size, extent, topology) =
 construct_grid(::Type{<:LatitudeLongitudeGrid}, arch, size, extent, topology) =
     LatitudeLongitudeGrid(arch; size, longitude = extent[1], latitude = extent[2], z = extent[3], topology)
 
+# Regrid a field in three dimensions
+function three_dimensional_regrid!(a, b)
+    target_grid = a.grid isa ImmersedBoundaryGrid ? a.grid.underlying_grid : a.grid
+    source_grid = b.grid isa ImmersedBoundaryGrid ? b.grid.underlying_grid : b.grid
+
+    topo = topology(target_grid)
+    arch = architecture(target_grid)
+    arch = child_architecture(arch)
+
+    target_y = yt = cpu_face_constructor_y(target_grid)
+    target_z = zt = cpu_face_constructor_z(target_grid)
+
+    target_size = Nt = size(target_grid)
+
+    source_x = xs = cpu_face_constructor_x(source_grid)
+    source_y = ys = cpu_face_constructor_y(source_grid)
+
+    source_size = Ns = size(source_grid)
+
+    # Start by regridding in z
+    @debug "Regridding in z"
+    zgrid   = construct_grid(typeof(target_grid), arch, (Ns[1], Ns[2], Nt[3]), (xs, ys, zt), topo)
+    field_z = Field(location(b), zgrid)
+    regrid!(field_z, zgrid, source_grid, b)
+
+    # regrid in y
+    @debug "Regridding in y"
+    ygrid   = construct_grid(typeof(target_grid), arch, (Ns[1], Nt[2], Nt[3]), (xs, yt, zt), topo)
+    field_y = Field(location(b), ygrid);
+    regrid!(field_y, ygrid, zgrid, field_z);
+
+    # Finally regrid in x
+    @debug "Regridding in x"
+    regrid!(a, target_grid, ygrid, field_y)
+
+    return a
+end
+
 include("diffuse_tracers.jl")
 
 end # module

src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl

@@ -14,7 +14,6 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model)
                       v = ZeroField(),
                       h = sea_ice.model.ice_thickness,
                       ℵ = sea_ice.model.ice_concentration,
-                      hc = sea_ice.model.ice_consolidation_thickness,
                       Tₒ = ocean.model.tracers.T,
                       Sₒ = ocean.model.tracers.S)
 
@@ -98,7 +97,6 @@ end
         uᵢ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u)
         vᵢ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v)
         hᵢ = interior_state.h[i, j, 1]
-        hc = interior_state.hc[i, j, 1]
         ℵᵢ = interior_state.ℵ[i, j, 1]
         Tₛ = interface_temperature[i, j, 1]
         Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ)
@@ -119,7 +117,7 @@ end
                               h_bℓ = atmosphere_state.h_bℓ)
 
     downwelling_radiation = (; Qs, Qℓ)
-    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, hc=hc)
+    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ)
 
     # Estimate initial interface state
     u★ = convert(FT, 1e-4)

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -291,7 +291,7 @@ function default_ai_temperature(sea_ice::SeaIceSimulation)
 end
 
 function default_ao_specific_humidity(ocean)
-    FT = Float64
+    FT = eltype(ocean)
     phase = AtmosphericThermodynamics.Liquid()
     x_H₂O = convert(FT, 0.98)
     return ImpureSaturationSpecificHumidity(phase, x_H₂O)
@@ -345,7 +345,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                              sea_ice_heat_capacity = heat_capacity(sea_ice),
                              gravitational_acceleration = g_Earth)
 
-    FT = Float64
+    FT = eltype(ocean)
 
     ocean_reference_density    = convert(FT, ocean_reference_density)
     ocean_heat_capacity        = convert(FT, ocean_heat_capacity)

src/OceanSeaIceModels/OceanSeaIceModels.jl

@@ -43,6 +43,8 @@ const default_freshwater_density = 1000
 const SeaIceSimulation = Simulation{<:SeaIceModel}
 const OceananigansSimulation = Simulation{<:HydrostaticFreeSurfaceModel}
 
+Base.eltype(ocean::OceananigansSimulation) = eltype(ocean.model.grid)
+
 sea_ice_thickness(::Nothing) = ZeroField()
 sea_ice_thickness(sea_ice::SeaIceSimulation) = sea_ice.model.ice_thickness
 

src/OceanSeaIceModels/ocean_sea_ice_model.jl

@@ -19,6 +19,8 @@ import Oceananigans.TimeSteppers: time_step!, update_state!, time
 import Oceananigans.Utils: prettytime
 import Oceananigans.Models: timestepper, NaNChecker, default_nan_checker, initialization_update_state!
 
+import Base
+
 mutable struct OceanSeaIceModel{I, A, O, F, C, Arch} <: AbstractModel{Nothing, Arch}
     architecture :: Arch
     clock :: C