odd behaviour with ice_concentration

[francispoulin]

Below is a sample code I put together for a 1d sea ice model that I put together for testing.

What is odd is that in this case I initialize the sea ice thicknesses and concentrations to 100 and 10, respectively. Then, after the first output, the thickness increases to 1000 and the concentration reduces to 1. Is there a normalization going on here?

using ClimaSeaIce
using ClimaSeaIce.SeaIceMomentumEquations
using ClimaSeaIce.Rheologies
using Oceananigans
using Oceananigans.Units 
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Printf

grid = RectilinearGrid(CPU();
                       size = (100), 
                          x = (0, 100kilometers),  
                       halo = (7),
                   topology = (Periodic, Flat, Flat))

# Atmosphere - sea ice stress
Va = YFaceField(grid)
τᵥ = Field(Va)
set!(Va, 1e-4*900*cos.(2*pi*grid.xᶠᵃᵃ[1:grid.Nx]/grid.Lx))
Oceananigans.BoundaryConditions.fill_halo_regions!(τᵥ)
compute!(τᵥ)

dynamics = SeaIceMomentumEquation(grid; 
                                  top_momentum_stress = (v=τᵥ),
                                  coriolis     = FPlane(f=1e-4),
                                  )

model = SeaIceModel(grid; 
                    dynamics,
                    ice_thermodynamics = nothing, # No ice_thermodynamics here
                    advection = WENO(order=7))

set!(model, h = 100)
set!(model, ℵ = 10)

simulation = Simulation(model, Δt = 2minutes, stop_time = hours)

h    = model.ice_thickness
ℵ    = model.ice_concentration
u, v = model.velocities

wall_time = [time_ns()]

function progress(sim) 
    h = sim.model.ice_thickness
    ℵ = sim.model.ice_concentration
    u = sim.model.velocities.u
    v = sim.model.velocities.v

    hmax = maximum(interior(h))
    ℵmin = minimum(interior(ℵ))
    umax = maximum(interior(u)), maximum(interior(v))
    step_time = 1e-9 * (time_ns() - wall_time[1])

    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e), max(h): %.2f, min(ℵ): %.2f, wtime: %s \n",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   prettytime(sim.Δt),
                   umax..., hmax, ℵmin, prettytime(step_time))

     wall_time[1] = time_ns()
end

simulation.callbacks[:progress] = Callback(progress, TimeInterval(1hours))

run!(simulation)

[simone-silvestri]

The concentration must be between 0 and 1, so to keep the volume constant (V = h * conc) the algorithm increases the thickness whenever concentration increases beyond 1 (and caps concentration)