Sea Surface Temperature Update

Updates SSTs in model dump files during an AUS2200 run. Useful for creating runs in which SST evolves over time according to known historical data.

import iris
import mule
import xarray as xr
import numpy as np

dump_fn  = '/scratch/ly62/dr4292/cylc-run/u-cs142-20160603T0000-continuous-140E-evolvingSST/share/cycle/20160606T1800Z/aus2200/d0198/RA3/ics/umnsaa_da006'
dump_out = '/scratch/ly62/dr4292/cylc-run/u-cs142-20160603T0000-continuous-140E-evolvingSST/share/cycle/20160606T1800Z/aus2200/d0198/RA3/ics/umnsaa_da006_newsst'
era5_fn  = '/scratch/ly62/dr4292/cylc-run/u-cs142-20160603T0000-continuous-140E-evolvingSST/share/cycle/20160606T1800Z/ec/um/ec_cb006'

cubes = iris.load(dump_fn)

/g/data/hh5/public/apps/cms_conda/envs/analysis3-24.01/lib/python3.10/site-packages/iris/fileformats/_ff.py:722: _WarnComboLoadingDefaulting: The stash code m01s00i060 is on a grid 22 which has not been explicitly handled by the fieldsfile loader. Assuming the data is on a P grid.
  warnings.warn(

era5 = iris.load(era5_fn)

/g/data/hh5/public/apps/cms_conda/envs/analysis3-24.01/lib/python3.10/site-packages/iris/fileformats/_ff.py:722: _WarnComboLoadingDefaulting: The stash code m01s00i060 is on a grid 22 which has not been explicitly handled by the fieldsfile loader. Assuming the data is on a P grid.
  warnings.warn(

for f in cubes:
    s=f.attributes['STASH']
    if s.model==1 and s.section==0 and s.item == 24:
        um_sst = xr.DataArray.from_iris(f)
        break

for f in era5:
    s=f.attributes['STASH']
    if s.model==1 and s.section==0 and s.item == 24:
        era5_sst = xr.DataArray.from_iris(f)
        break

era5_sst

<xarray.DataArray 'surface_temperature' (latitude: 2170, longitude: 2650)> Size: 23MB
dask.array<filled, shape=(2170, 2650), dtype=float32, chunksize=(2170, 2650), chunktype=numpy.ndarray>
Coordinates:
  * latitude      (latitude) float32 9kB -49.28 -49.26 -49.24 ... -6.357 -6.337
  * longitude     (longitude) float32 11kB 113.8 113.8 113.8 ... 166.2 166.2
    pseudo_level  int64 8B ...
    realization   int64 8B ...
    time          datetime64[ns] 8B ...
Attributes:
    standard_name:  surface_temperature
    units:          K
    source:         Data from Met Office Unified Model
    um_version:     12.2
    STASH:          m01s00i024

xarray.DataArray

'surface_temperature'

• latitude: 2170
• longitude: 2650

• dask.array<chunksize=(2170, 2650), meta=np.ndarray>

  Array Chunk Bytes 21.94 MiB 21.94 MiB Shape (2170, 2650) (2170, 2650) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

• Coordinates: (5)
  • latitude
    (latitude)
    float32
    -49.28 -49.26 ... -6.357 -6.337

    standard_name :

    latitude

    units :

    degrees

    array([-49.283   , -49.263203, -49.2434  , ...,  -6.3764  ,  -6.356602,
            -6.3368  ], dtype=float32)

  • longitude
    (longitude)
    float32
    113.8 113.8 113.8 ... 166.2 166.2

    standard_name :

    longitude

    units :

    degrees

    array([113.765   , 113.7848  , 113.804596, ..., 166.1756  , 166.1954  ,
           166.2152  ], dtype=float32)

  • pseudo_level
    ()
    int64
    ...

    long_name :

    pseudo_level

    [1 values with dtype=int64]

  • realization
    ()
    int64
    ...

    standard_name :

    realization

    [1 values with dtype=int64]

  • time
    ()
    datetime64[ns]
    ...

    standard_name :

    time

    [1 values with dtype=datetime64[ns]]

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Index([ -49.28300094604492,  -49.26320266723633,  -49.24340057373047,
           -49.223602294921875, -49.203800201416016,  -49.18400192260742,
            -49.16419982910156,  -49.14440155029297,  -49.12459945678711,
           -49.104801177978516,
           ...
           -6.5149993896484375,  -6.495201110839844,   -6.47540283203125,
            -6.455600738525391,  -6.435802459716797, -6.4160003662109375,
            -6.396202087402344,  -6.376399993896484,  -6.356601715087891,
            -6.336799621582031],
          dtype='float32', name='latitude', length=2170))

  • longitude
    PandasIndex

    PandasIndex(Index([113.76499938964844, 113.78479766845703, 113.80459594726562,
           113.82440185546875, 113.84420013427734, 113.86399841308594,
           113.88379669189453, 113.90360260009766, 113.92340087890625,
           113.94319915771484,
           ...
           166.03700256347656, 166.05679321289062, 166.07659912109375,
           166.09640502929688, 166.11619567871094, 166.13600158691406,
           166.15579223632812, 166.17559814453125, 166.19540405273438,
           166.21519470214844],
          dtype='float32', name='longitude', length=2650))

• Attributes: (5)

  standard_name :

  surface_temperature

  units :

  K

  source :

  Data from Met Office Unified Model

  um_version :

  12.2

  STASH :

  m01s00i024

um_sst

<xarray.DataArray 'surface_temperature' (latitude: 2120, longitude: 2600)> Size: 22MB
dask.array<filled, shape=(2120, 2600), dtype=float32, chunksize=(2120, 2600), chunktype=numpy.ndarray>
Coordinates:
  * latitude      (latitude) float32 8kB -48.79 -48.77 -48.75 ... -6.852 -6.832
  * longitude     (longitude) float32 10kB 114.3 114.3 114.3 ... 165.7 165.7
    pseudo_level  int64 8B ...
    realization   int64 8B ...
    time          datetime64[ns] 8B ...
Attributes:
    standard_name:  surface_temperature
    units:          K
    source:         Data from Met Office Unified Model
    um_version:     12.2
    STASH:          m01s00i024

xarray.DataArray

'surface_temperature'

• latitude: 2120
• longitude: 2600

• dask.array<chunksize=(2120, 2600), meta=np.ndarray>

  Array Chunk Bytes 21.03 MiB 21.03 MiB Shape (2120, 2600) (2120, 2600) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

• Coordinates: (5)
  • latitude
    (latitude)
    float32
    -48.79 -48.77 ... -6.852 -6.832

    standard_name :

    latitude

    units :

    degrees

    array([-48.788   , -48.7682  , -48.748398, ...,  -6.871399,  -6.851597,
            -6.831799], dtype=float32)

  • longitude
    (longitude)
    float32
    114.3 114.3 114.3 ... 165.7 165.7

    standard_name :

    longitude

    units :

    degrees

    array([114.26   , 114.2798 , 114.2996 , ..., 165.6806 , 165.70041, 165.7202 ],
          dtype=float32)

  • pseudo_level
    ()
    int64
    ...

    long_name :

    pseudo_level

    [1 values with dtype=int64]

  • realization
    ()
    int64
    ...

    standard_name :

    realization

    [1 values with dtype=int64]

  • time
    ()
    datetime64[ns]
    ...

    standard_name :

    time

    [1 values with dtype=datetime64[ns]]

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Index([ -48.78799819946289,   -48.7681999206543,  -48.74839782714844,
           -48.728599548339844, -48.708797454833984,  -48.68899917602539,
            -48.66919708251953,  -48.64939880371094,  -48.62959671020508,
           -48.609798431396484,
           ...
            -7.009998321533203,  -6.990200042724609,   -6.97039794921875,
            -6.950599670410156,  -6.930797576904297,  -6.910999298095703,
            -6.891197204589844,   -6.87139892578125,  -6.851596832275391,
            -6.831798553466797],
          dtype='float32', name='latitude', length=2120))

  • longitude
    PandasIndex

    PandasIndex(Index([114.26000213623047, 114.27980041503906, 114.29959869384766,
           114.31940460205078, 114.33920288085938, 114.35900115966797,
           114.37879943847656, 114.39860534667969, 114.41840362548828,
           114.43820190429688,
           ...
           165.54200744628906, 165.56179809570312, 165.58160400390625,
           165.60140991210938, 165.62120056152344, 165.64100646972656,
           165.66079711914062, 165.68060302734375, 165.70040893554688,
           165.72019958496094],
          dtype='float32', name='longitude', length=2600))

• Attributes: (5)

  standard_name :

  surface_temperature

  units :

  K

  source :

  Data from Met Office Unified Model

  um_version :

  12.2

  STASH :

  m01s00i024

era5_sst_on_um = era5_sst.isel(longitude=slice(25,-25),latitude=slice(25,-25))

era5_sst_on_um = era5_sst_on_um.assign_coords(latitude=um_sst.latitude,longitude=um_sst.longitude)

for f in cubes:
    s=f.attributes['STASH']
    if s.model==1 and s.section==0 and s.item == 30:
        um_lsm = xr.DataArray.from_iris(f)
        break

um_lsm

<xarray.DataArray 'land_binary_mask' (latitude: 2120, longitude: 2600)> Size: 44MB
dask.array<filled, shape=(2120, 2600), dtype=int64, chunksize=(2120, 2600), chunktype=numpy.ndarray>
Coordinates:
  * latitude      (latitude) float32 8kB -48.79 -48.77 -48.75 ... -6.852 -6.832
  * longitude     (longitude) float32 10kB 114.3 114.3 114.3 ... 165.7 165.7
    pseudo_level  int64 8B ...
    realization   int64 8B ...
    time          datetime64[ns] 8B ...
Attributes:
    standard_name:  land_binary_mask
    source:         Data from Met Office Unified Model
    um_version:     12.2
    STASH:          m01s00i030

xarray.DataArray

'land_binary_mask'

• latitude: 2120
• longitude: 2600

• dask.array<chunksize=(2120, 2600), meta=np.ndarray>

  Array Chunk Bytes 42.05 MiB 42.05 MiB Shape (2120, 2600) (2120, 2600) Dask graph 1 chunks in 3 graph layers Data type int64 numpy.ndarray

• Coordinates: (5)
  • latitude
    (latitude)
    float32
    -48.79 -48.77 ... -6.852 -6.832

    standard_name :

    latitude

    units :

    degrees

    array([-48.788   , -48.7682  , -48.748398, ...,  -6.871399,  -6.851597,
            -6.831799], dtype=float32)

  • longitude
    (longitude)
    float32
    114.3 114.3 114.3 ... 165.7 165.7

    standard_name :

    longitude

    units :

    degrees

    array([114.26   , 114.2798 , 114.2996 , ..., 165.6806 , 165.70041, 165.7202 ],
          dtype=float32)

  • pseudo_level
    ()
    int64
    ...

    long_name :

    pseudo_level

    [1 values with dtype=int64]

  • realization
    ()
    int64
    ...

    standard_name :

    realization

    [1 values with dtype=int64]

  • time
    ()
    datetime64[ns]
    ...

    standard_name :

    time

    [1 values with dtype=datetime64[ns]]

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Index([ -48.78799819946289,   -48.7681999206543,  -48.74839782714844,
           -48.728599548339844, -48.708797454833984,  -48.68899917602539,
            -48.66919708251953,  -48.64939880371094,  -48.62959671020508,
           -48.609798431396484,
           ...
            -7.009998321533203,  -6.990200042724609,   -6.97039794921875,
            -6.950599670410156,  -6.930797576904297,  -6.910999298095703,
            -6.891197204589844,   -6.87139892578125,  -6.851596832275391,
            -6.831798553466797],
          dtype='float32', name='latitude', length=2120))

  • longitude
    PandasIndex

    PandasIndex(Index([114.26000213623047, 114.27980041503906, 114.29959869384766,
           114.31940460205078, 114.33920288085938, 114.35900115966797,
           114.37879943847656, 114.39860534667969, 114.41840362548828,
           114.43820190429688,
           ...
           165.54200744628906, 165.56179809570312, 165.58160400390625,
           165.60140991210938, 165.62120056152344, 165.64100646972656,
           165.66079711914062, 165.68060302734375, 165.70040893554688,
           165.72019958496094],
          dtype='float32', name='longitude', length=2600))

• Attributes: (4)

  standard_name :

  land_binary_mask

  source :

  Data from Met Office Unified Model

  um_version :

  12.2

  STASH :

  m01s00i030

updated_sst=(((1-um_lsm)*era5_sst_on_um) + um_lsm*um_sst)

### Sanity check - identical over land
(um_lsm*(um_sst-updated_sst)).plot()

<matplotlib.collections.QuadMesh at 0x1534a9516f20>

### Sanity check - different over ocean
((1-um_lsm)*(um_sst-updated_sst)).plot()

<matplotlib.collections.QuadMesh at 0x1534a48cd5d0>

d=updated_sst.data.compute()

dump_mule = mule.load_umfile(dump_fn)

/g/data/hh5/public/apps/cms_conda/envs/analysis3-24.01/lib/python3.10/site-packages/mule/stashmaster.py:259: UserWarning:
Unable to load STASHmaster from version string, path does not exist
Path: $UMDIR/vn12.2/ctldata/STASHmaster/STASHmaster_A
Please check that the value of mule.stashmaster.STASHMASTER_PATH_PATTERN is correct for your site/configuration
  warnings.warn(msg)

mf_out = dump_mule.copy()

array_provider = mule.ArrayDataProvider(d)

for f in dump_mule.fields:
    if f.lbuser4 != 24:
        mf_out.fields.append(f)
    else:
        new_field = f.copy()
        new_field.set_data_provider(array_provider)
        mf_out.fields.append(new_field)

mf_out.to_file(dump_out)

new_cubes = iris.load(dump_out)

/g/data/hh5/public/apps/cms_conda/envs/analysis3-24.01/lib/python3.10/site-packages/iris/fileformats/_ff.py:722: _WarnComboLoadingDefaulting: The stash code m01s00i060 is on a grid 22 which has not been explicitly handled by the fieldsfile loader. Assuming the data is on a P grid.
  warnings.warn(

for f in new_cubes:
    s=f.attributes['STASH']
    if s.model==1 and s.section==0 and s.item == 24:
        new_um_sst = xr.DataArray.from_iris(f)
        break

### Sanity check - identical over land
(um_lsm*(um_sst-new_um_sst)).plot()

<matplotlib.collections.QuadMesh at 0x1534a8f61de0>

### Sanity check - different over ocean
((1-um_lsm)*(um_sst-new_um_sst)).plot()

<matplotlib.collections.QuadMesh at 0x1534a8d5ecb0>