SSTs in flat earth experiment

[torsten20]

When I plot the SSTs from my flat earth simulation (Pre-Defined_Experiments/Flat_Earth.ipynb) it also shows values over land. Is that just the surface temperature everywhere or are the SSTs over land to be discarded?

[milankl]

No we don't combine land surface temperature and sea surface temperature. Those should be separate variables in the netCDF output. Both should come with a mask though, do they not? You have some flexibility/options on this so it depends a little bit of what you did. Can you paste in the script?

[torsten20]

In my output file the sst has only dimensions time, lat, lon. I have a variable lsm (land sea mask) in the output. So I could use that to only plot ocean points. But the sst field itself contains data over land. I have another variable st (soil temperature) in my output. This one has dimensions time, soil layer, lat, lon. But there is no surface temperature variable for land specifically.

My model is based on Pier Luigi's climate run, just with orography removed. The code is in this script

It basically comes down to this:

using SpeedyWeather
# components
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
output = NetCDFOutput(spectral_grid, PrimitiveWetModel, path="/home/hy913279/summerschool2025/", id="Flat_Earth", output_dt=Minute(180))
add!(output, SpeedyWeather.TemperatureOutput(), SpeedyWeather.LargeScalePrecipitationOutput(), SpeedyWeather.ConvectivePrecipitationOutput())  # output also temperature and precip

# Some parameterisation choices
vertical_advection = WENOVerticalAdvection(spectral_grid)
large_scale_condensation = ImplicitCondensation(spectral_grid)
convection = SimplifiedBettsMiller(spectral_grid)
temperature=LandBucketTemperature(spectral_grid)
soil_moisture=LandBucketMoisture(spectral_grid)
land = LandModel(spectral_grid; temperature, soil_moisture)
ocean = SlabOcean(spectral_grid, mixed_layer_depth=50)

#Try to fix ocean albedo
albedo = Albedo(GlobalConstantAlbedo(spectral_grid, albedo=0.06), AlbedoClimatology(spectral_grid))

orography=NoOrography(spectral_grid)

time_stepping = Leapfrog(spectral_grid, Δt_at_T31=Minute(40)) #CPL this is the standar time step

# create model, initialize, run
model = PrimitiveWetModel(spectral_grid; orography, vertical_advection, time_stepping, large_scale_condensation, convection, land, ocean, albedo, output=output)

add!(model, SpeedyWeather.AllOutputVariables()...)
simulation = initialize!(model)
set!(albedo.land, (λ, φ) -> 0., add=true)

run!(simulation, period=Day(30*365), output=true)

[milankl]

You are right that currently soil temperatures (after some 6 months)

and sea surface temperatures are not masked

and I agree that this is confusing! You are right though that the land-sea mask is provided which in principle lets you mask in a post-processing fashion but that's ideally not something the user should bother with. Applying a mask before output I do have some constraints though that I will discuss in a SpeedyWeather issue which will be linked here in a bit!

[milankl]

Opened an issue here ☝🏼 feel free to add your opinion on this. For the time being use the land-sea mask, pick a threshold and mask SST / soil temperature manually. What you currently see is initial conditions with no changes when a grid point is fully ocean or fully land

[milankl]

On a related note, because we don't mask soil moisture either but @pierluigividale suggested to start with soil moisture of 1 as default initial conditions we also have soil moisture being 1 over the ocean, which kind of makes sense I guess???

[pierluigividale]

That is a funny one. Yes, it sort of makes sense.

Now, with the change I had to make to the soil levels, with the first soil layer having a thickness of 0.5m, it is a bit harder to get rid of moisture near the surface, but eventually it does happen anyway. Previously, with a thin top soil, the adjustment was fast, and then it was just up to vegetation to cause an adjustment lower down. All are compromises.

[milankl]

I forgot that we already have a solution for this, you need to set mask=true for all components that you want masked in the initial conditions, e.g.

temperature = LandBucketTemperature(spectral_grid, mask=true)
soil_moisture = LandBucketMoisture(spectral_grid, mask=true)
land = LandModel(spectral_grid; temperature, soil_moisture)
ocean = SlabOcean(spectral_grid, mixed_layer_depth=50, mask=true)

this will set ocean grid points with only land according to the land-sea mask to NaN and vice versa. Looks therefore like

(I think ncvis uses the lowermost colour to replace masked ... but the value is definitely set nc-wise to masked). But we seem to have an interpolation problem with NaNs on land

not sure where this is from!

[torsten20]

Thanks for all the info @milankl. I will add mask=true in my model.

[milankl]

Great! But note that this is a somewhat philosophical question: A masked array in netCDF or also in numpy e.g.

In [21]: np.ma.masked_array(np.random.rand(2,3), mask=np.random.rand(2, 3) > 0.5)
Out[21]:
masked_array(
  data=[[--, 0.9585668056536504, 0.8779267129117457],
        [0.48460071802898297, 0.1928554137579932, --]],
  mask=[[ True, False, False],
        [False, False,  True]],
  fill_value=1e+20)

has a boolean mask, not a fractional mask that we use. So how to convert a fractional mask to a boolean has a free parameter, the threshold (maybe 0.5? maybe 0, 1?), so it's up to you how to do that. With mask=true as suggested above we use 1 as a threshold for ocean (mask points that are only land) and 0 as a threshold for land (mask points that are only ocean)

[torsten20]

I ran my simulation with mask=true and got the masked SST and ST fields. However, when I try to combine both fields into one combined surface temperature field, I am left with 6 grid cells (in a T31 simulation) which are still NaN.

I combined both fields with:
T_flat=np.where(np.isnan(ds_flat.st[:,0,:,:]),ds_flat.sst,ds_flat.st[:,0,:,:])

Plotting T_flat for one specific time step, I get:

I looked at several different timesteps and the NaNs are always in the same location. When I plot st and sst separately I can see that the masks don't match in those locations.

[pierluigividale]

I can see that the joys of masking are not exclusive to the UM.

[torsten20]

I think I will use the landmask from the output file since it also allows to consider hybrid land/ocean points. And it is not really a lot of extra effort.

[milankl]

So this problem does not occur when using the original octahedral Gaussian grid

# check whether SST or land is finite (=not masked)
julia> m = isfinite.(simulation.prognostic_variables.ocean.sea_surface_temperature) .|
isfinite.(simulation.prognostic_variables.land.soil_temperature[:, 1])

julia> all(m .== true)
true

So all grid points are either SST or land or both. (Replacing .| (or) with .& (and) and sum(m)/length(m) tells me some 41% of the grid points are both ocean and land!)

Meaning that this is likely a problem of interpolating one field with NaN values onto another grid. But as you say if you want to combine them into one surface temperature field then I would indeed not mask them and weight them according to the mask! You technically still want to do that on the original grid before interpolation but that means you'd need to define a new output variable that does that for you on the fly instead of dealing with the netCDF output later!