Setting up a coupled model at a new resolution

Setting up a coupled model at a new resolution#

Scott Wales, CLEX CMS

We’ve been setting up an ACCESS coupled model for a new atmosphere resolution. The process is quite involved, here’s some notes on how we’ve done it.

ACCESS consists of three models - the UM atmosphere, MOM ocean and CICE sea ice models, tied together with the Oasis coupler. The coupler handles sending fields between all of the models, to do so it needs to regrid the data fields from the grid of one model to the grid of another.

Regridding weights#

Oasis requires regridding weights in SCRIP format to convert a field from one model to another. Oasis can generate these itself from SCRIP descriptions of the source and target grid, or they can be generated using ESMF_RegridWeightGen and the ESMF weights converted to SCRIP format.

The most important thing here is that the masks are consistent between all the models. If they are not you will see discontinuity in the surface fields around coastlines.

UM Grids#

The UM mask is defined by two files - the land mask, which is 1 everywhere there is greater than 0% land in a grid cell, and the land fraction, which is the fraction of land in a grid cell. These fields are only defined on the T grid - the UM uses an Arakawa C grid, with the scalar component termed the ‘T’ grid and the vector components termed the ‘U’ and ‘V’ grids.

The coupling masks define the regions on the source and target grids which are involved in the coupling. Since we’re coupling ocean and sea ice models we want to mask out land points where the ocean fields are not defined. The mask for the UM T grid should be 1 everywhere the land fraction field is less than 1, in other words everywhere with greater than 0% ocean in a cell.

For the U and V fields I used no mask, since we don’t want to conservatively interpolate vector fields and masking the offset grids is more difficult.

Since the UM uses regular lat-lon grids you can create a GRIDSPEC grid description to feed into ESMF.

MOM/CICE Grids#

The MOM and CICE masks should be equal to each other. The CICE mask is found in $ICE/INPUT/kmt.nc as variable kmt, the MOM mask in $OCN/INPUT/grid_spec.nc as variable wet.

The coupling mask for the MOM/CICE T grid should be equal to wet from grid_spec.nc.

MOM uses a tripolar grid, all of the fields needed to create a SCRIP description of a MOM grid may be found in grid_spec.nc

Masks and Conservative Regridding#

Whereever possible we want to be using conservative regridding - we don’t want our model to loose energy each timestep otherwise fields like the model temperature will drift. Conservative regridding ensures that the total field value sent by the source model equals the total field value received by the destination model.

For conservative regridding to work, the masks of the two models must be consistent as well - there can’t be areas where the ocean model thinks it’s land and the atmosphere model thinks it’s ocean.

An example point where this can happen is the ice shelves around Antarctica. The ocean model’s grid doesn’t go all the way to the south pole, so it misses the southernmost part of the Ross ice shelf. The atmosphere grid does reach the pole however, and its default land mask assumes the ice shelf is ocean. There can be other areas where this causes conflicts as well - for instance the black sea is not included in the MOM default mask.

To make sure that the masks are consistent you need to do a conservative regridding of the mask from one grid to the other. The highest resolution grid should be the source, in my case I used the ocean mask. I then conservatively regridded the wet field from MOM to the UM grid with ESMF_RegridWeightGen. It’s important to not use masks for this interpolation, since the mask itself is what we’re interpolating and we want to conserve the land area. Also since the MOM grid doesn’t cover the entire globe I turned on the ESMF ‘–ignore_unmapped’ flag.

This operation gave me an ocean fraction field on the UM grid. The ocean fraction can then be straightforwardly converted to a land fraction and land mask field to be input to the UM, although note that the interpolation can introduce a bit of noise, so the land fraction won’t be exactly 1.0 on pure land points - I used all grid cells where the land fraction was greater than 0.99 as the land mask.

With a new land mask and land fraction the UM ancillaries needed to be modified

Regridding Weights#

To calculate regridding weights I used ESMF_RegridWeightGen, a command-line tool included with ESMF. ESMF names the fields differently to SCRIP, so some modification of the output is required in order for OASIS to be able to read the files.

Conservative regridding is used for scalar fields, for vector fields I used patch regridding to for smooth interpolation without introducing artefacts from the different grid resolutions.

UM Ancillary interpolation#

Since the UM is using a new resolution and mask the various ancillary datasets also need to be changed. For the most part files from the Met Office archive at the target resolution were used. For fields like vegetation fraction which are only defined on land points (and missing data over the ocean) I ‘de-masked’ the ancillary field, by performing a nearest grid point interpolation from the masked source field to the same grid without a mask to remove the missing data areas. It’s important to use nearest grid point so that fields like the vegetation fractions correctly add up - some interpolation types also can result in negative values at some points which crashes the model.

The aerosol climatology files had to be vertically interpolated since we had also reduced the number of vertical levels in the model run. For this I used the stratify python library. The vertical levels used by the UM are hybrid height levels - they vary depending on location, following the terrain up to a certain height, then above that level they are constant height levels.

To convert from the source set of hybrid heights to the target hybrid heights I first converted both to true heights, using the formula from UMDB F03. I then used stratify to convert the aerosol fields on true source heights to the true target heights.

To modify the UM ancillary files I used the mule library as well as pandas to group the individual 2d lat-lon field slices that are contained in UM format files into the full 4d field

def categorise_fields(m):
    """
    Sorts all of the 2d field slices in mule file ``m`` into a
    pandas.DataFrame so that they can be grouped together into full
    fields
    """
    df = pandas.DataFrame({'field': m.fields})
    df['year'] = df['field'].apply(lambda f: f.lbyr)
    df['month'] = df['field'].apply(lambda f: f.lbmon)
    df['day'] = df['field'].apply(lambda f: f.lbdat)
    df['hour'] = df['field'].apply(lambda f: f.lbhr)
    df['minute'] = df['field'].apply(lambda f: f.lbmin)
    df['second'] = df['field'].apply(lambda f: f.lbsec)
    df['stash'] = df['field'].apply(lambda f: f.lbuser4)
    df['vertical_type'] =df['field'].apply(lambda f: f.lbvc)
    df['level'] = df['field'].apply(lambda f: f.lblev)
    df['pseudo'] = df['field'].apply(lambda f: f.lbuser5)

    df['blev'] = df['field'].apply(lambda f: f.blev)
    df['brlev'] = df['field'].apply(lambda f: f.brlev)

    df['bhlev'] = df['field'].apply(lambda f: f.bhlev)
    df['bhrlev'] = df['field'].apply(lambda f: f.bhrlev)

    return df

# Group the 2d slices with the same field and time value together
ancil = mule.AncilFile.from_file(ancil_path)
df = categorise_fields(ancil)
for name, g in df.groupby(['year','month','day','hour','minute','second', 'stash']):
    print("%04d%02d%02dT%02d:%02d:%02d STASH %d"%name)

    # Stack the slices of this field/time into a 3d array
    cube = numpy.stack(g['field'].apply(lambda f: f.get_data()))

def true_theta_height(level_file, orography):
    """
    Create a true height field from the UM levels namelist at ``level_file``
    and the orography field
    """
    target_levels = f90nml.read(vertlevs)['VERTLEVS']
    
    eta = numpy.array(target_levels['eta_theta'])
    const_lev = target_levels['first_constant_r_rho_level']-1

    Zsea = target_levels['z_top_of_model'] * eta
    C = (1 - eta/eta[const_lev])**2
    C[const_lev:] = 0
    Z = target_Zsea[:, numpy.newaxis, numpy.newaxis] + numpy.multiply.outer(target_C,orog)

    return Z

Running the model#

With all of this set up the next thing to do is configure the model with all of the new files then do a short test run. A single model day should tell you if the coupling fields are being sent correctly and all of the ancillary files are at the correct resolution. With that done you can move on to longer runs to check the stability of the model.