Surface restoring

Regrid surface salinity and potential temperature from land filled WOA dataset to the ocean model grid.

import xarray as xr
import numpy as np
from datetime import datetime
import os, subprocess
import xesmf

Read tx1_4 grid

ds_out = xr.open_dataset('../mesh/tx1_4_grid.nc').rename({'tlon': 'lon','tlat': 'lat', 'qlon': 'lon_b','qlat': 'lat_b',})

Read WOA file with land fill

# WOA sfc state file with land fill, created using create_filled_sfc.py
fname = '/glade/work/gmarques/cesm/mom6_input_files/WOA_MOM6/woa18_sfc_state_filled.nc'
woa = xr.open_dataset(fname, decode_times=False)

# average between two-layers (depth = 0 and depth = 10, depth indices 0 and 2)
woa_s_an_surface_ave = woa.s_an.isel(depth=[0,2]).mean('depth')
woa_theta0_surface_ave = woa.theta0.isel(depth=[0,2]).mean('depth')

def regrid_tracer(fld, ds_in, ds_out, method='bilinear'):

    regrid = xesmf.Regridder(
        ds_in,
        ds_out,
        method=method,
        periodic=True,
    )
    fld_out = regrid(fld)
    return fld_out

ds_out_s_an = regrid_tracer(woa_s_an_surface_ave, ds_in, ds_out)

ds_out_theta0 = regrid_tracer(woa_theta0_surface_ave, ds_in, ds_out)

Create state file for MOM6

We opted to create this file via ncgen to avoid issues with FMS reading the netCDF file.

!ncgen -o state_restore_tx1_4.nc state_restore_tx1_4.cdl

state = xr.open_dataset('state_restore_tx1_4.nc', decode_times=False)

ds_out

<xarray.Dataset>
Dimensions:  (ny: 1080, nx: 1440, nxp: 1441, nyp: 1081)
Dimensions without coordinates: ny, nx, nxp, nyp
Data variables: (12/20)
    lon      (ny, nx) float64 -286.9 -286.6 -286.4 -286.1 ... 73.0 73.0 73.0
    lat      (ny, nx) float64 -80.03 -80.03 -80.03 -80.03 ... 50.03 50.0 50.0
    ulon     (ny, nxp) float64 -287.0 -286.8 -286.5 -286.2 ... 73.0 73.0 73.0
    ulat     (ny, nxp) float64 -80.03 -80.03 -80.03 -80.03 ... 50.01 50.0 50.0
    vlon     (nyp, nx) float64 -286.9 -286.6 -286.4 -286.1 ... 73.0 73.0 73.0
    vlat     (nyp, nx) float64 -80.05 -80.05 -80.05 -80.05 ... 50.03 50.0 50.0
    ...       ...
    tarea    (ny, nx) float64 2.315e+07 2.315e+07 2.315e+07 ... 3.339e+05 198.9
    tmask    (ny, nx) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    angle    (ny, nx) float64 -0.06156 -0.06156 -0.06156 ... -4.024 -0.4197
    depth    (ny, nx) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    ar       (ny, nx) float64 1.0 1.0 1.0 1.0 ... 0.1318 0.1046 0.281 1.457e+03
    egs      (ny, nx) float64 0.0004327 0.0004327 ... 5.196e-05 1.268e-06
Attributes:
    Description:  CESM MOM6 1/4 degree grid
    Author:       Frank, Fred, Gustavo (gmarques@ucar.edu)
    Created:      2022-12-27T10:20:11.076318
    type:         Glogal 1/4 degree grid file

xarray.Dataset

• Dimensions:
  • ny: 1080
  • nx: 1440
  • nxp: 1441
  • nyp: 1081
• Coordinates: (0)
• Data variables: (20)
  • lon
    (ny, nx)
    float64
    -286.9 -286.6 -286.4 ... 73.0 73.0

    long_name :

    array of t-grid longitudes

    units :

    degrees_east

    array([[-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           [-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           [-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           ...,
           [-286.99985 , -286.992491, -286.980779, ...,   72.980779,   72.992491,
              72.99985 ],
           [-286.999995, -286.997094, -286.989071, ...,   72.989071,   72.997094,
              72.999995],
           [-286.999998, -286.999988, -286.996584, ...,   72.996584,   72.999988,
              72.999998]])

  • lat
    (ny, nx)
    float64
    -80.03 -80.03 -80.03 ... 50.0 50.0

    long_name :

    array of t-grid latitudes

    units :

    degrees_north

    array([[-80.028331, -80.028331, -80.028331, ..., -80.028331, -80.028331,
            -80.028331],
           [-79.984948, -79.984948, -79.984948, ..., -79.984948, -79.984948,
            -79.984948],
           [-79.941378, -79.941378, -79.941378, ..., -79.941378, -79.941378,
            -79.941378],
           ...,
           [ 49.987955,  49.994464,  50.030566, ...,  50.030566,  49.994464,
             49.987955],
           [ 49.99279 ,  49.995509,  50.0312  , ...,  50.0312  ,  49.995509,
             49.99279 ],
           [ 49.997603,  49.997603,  50.031502, ...,  50.031502,  49.997603,
             49.997603]])

  • ulon
    (ny, nxp)
    float64
    -287.0 -286.8 -286.5 ... 73.0 73.0

    long_name :

    array of u-grid longitudes

    units :

    degrees_east

    array([[-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           [-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           [-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           ...,
           [-287.      , -286.998953, -286.986114, ...,   72.986114,   72.998953,
              73.      ],
           [-287.      , -286.999988, -286.992498, ...,   72.992498,   72.999988,
              73.      ],
           [-287.      , -286.999996, -286.997835, ...,   72.997835,   72.999996,
              73.      ]])

  • ulat
    (ny, nxp)
    float64
    -80.03 -80.03 -80.03 ... 50.0 50.0

    long_name :

    array of u-grid latitudes

    units :

    degrees_north

    array([[-80.028331, -80.028331, -80.028331, ..., -80.028331, -80.028331,
            -80.028331],
           [-79.984948, -79.984948, -79.984948, ..., -79.984948, -79.984948,
            -79.984948],
           [-79.941378, -79.941378, -79.941378, ..., -79.941378, -79.941378,
            -79.941378],
           ...,
           [ 49.987952,  49.988079,  50.010207, ...,  50.010207,  49.988079,
             49.987952],
           [ 49.99279 ,  49.99279 ,  50.010899, ...,  50.010899,  49.99279 ,
             49.99279 ],
           [ 49.997603,  49.997603,  50.011226, ...,  50.011226,  49.997603,
             49.997603]])

  • vlon
    (nyp, nx)
    float64
    -286.9 -286.6 -286.4 ... 73.0 73.0

    long_name :

    array of v-grid longitudes

    units :

    degrees_east

    array([[-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           [-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           [-286.875   , -286.625   , -286.375   , ...,   72.375   ,   72.625   ,
              72.875   ],
           ...,
           [-286.999972, -286.994926, -286.985014, ...,   72.985014,   72.994926,
              72.999972],
           [-286.999997, -286.99898 , -286.992934, ...,   72.992934,   72.99898 ,
              72.999997],
           [-287.      , -287.      , -287.      , ...,   73.      ,   73.      ,
              73.      ]])

  • vlat
    (nyp, nx)
    float64
    -80.05 -80.05 -80.05 ... 50.0 50.0

    long_name :

    array of v-grid latitudes

    units :

    degrees_north

    array([[-80.049953, -80.049953, -80.049953, ..., -80.049953, -80.049953,
            -80.049953],
           [-80.006663, -80.006663, -80.006663, ..., -80.006663, -80.006663,
            -80.006663],
           [-79.963186, -79.963186, -79.963186, ..., -79.963186, -79.963186,
            -79.963186],
           ...,
           [ 49.990374,  49.995028,  50.030925, ...,  50.030925,  49.995028,
             49.990374],
           [ 49.995199,  49.995955,  50.031392, ...,  50.031392,  49.995955,
             49.995199],
           [ 50.      ,  50.      ,  50.03154 , ...,  50.03154 ,  50.      ,
             50.      ]])

  • lon_b
    (nyp, nxp)
    float64
    -287.0 -286.8 -286.5 ... 73.0 73.0

    long_name :

    array of q-grid longitudes

    units :

    degrees_east

    array([[-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           [-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           [-287.      , -286.75    , -286.5     , ...,   72.5     ,   72.75    ,
              73.      ],
           ...,
           [-287.      , -286.999792, -286.98942 , ...,   72.98942 ,   72.999792,
              73.      ],
           [-287.      , -286.999993, -286.995317, ...,   72.995317,   72.999993,
              73.      ],
           [-287.      , -287.      , -287.      , ...,   73.      ,   73.      ,
              73.      ]])

  • lat_b
    (nyp, nxp)
    float64
    -80.05 -80.05 -80.05 ... 50.0 50.0

    long_name :

    array of q-grid latitudes

    units :

    degrees_north

    array([[-80.049953, -80.049953, -80.049953, ..., -80.049953, -80.049953,
            -80.049953],
           [-80.006663, -80.006663, -80.006663, ..., -80.006663, -80.006663,
            -80.006663],
           [-79.963186, -79.963186, -79.963186, ..., -79.963186, -79.963186,
            -79.963186],
           ...,
           [ 49.990374,  49.990382,  50.0106  , ...,  50.0106  ,  49.990382,
             49.990374],
           [ 49.995199,  49.995199,  50.011105, ...,  50.011105,  49.995199,
             49.995199],
           [ 50.      ,  50.      ,  50.011264, ...,  50.011264,  50.      ,
             50.      ]])

  • dxt
    (ny, nx)
    float64
    4.814e+03 4.814e+03 ... 0.2751

    long_name :

    x-distance between u-points, centered at t

    units :

    meters

    array([[4.813835e+03, 4.813835e+03, 4.813835e+03, ..., 4.813835e+03,
            4.813835e+03, 4.813835e+03],
           [4.834565e+03, 4.834565e+03, 4.834565e+03, ..., 4.834565e+03,
            4.834565e+03, 4.834565e+03],
           [4.855381e+03, 4.855381e+03, 4.855381e+03, ..., 4.855381e+03,
            4.855381e+03, 4.855381e+03],
           ...,
           [1.304876e+02, 2.621169e+03, 4.897064e+03, ..., 4.897064e+03,
            2.621169e+03, 1.304876e+02],
           [1.458182e+00, 2.359978e+03, 4.854614e+03, ..., 4.854614e+03,
            2.359978e+03, 1.458182e+00],
           [2.751040e-01, 1.199777e+03, 4.831483e+03, ..., 4.831483e+03,
            1.199777e+03, 2.751040e-01]])

  • dyt
    (ny, nx)
    float64
    4.814e+03 4.814e+03 ... 337.2 401.0

    long_name :

    y-distance between v-points, centered at t

    units :

    meters

    array([[4813.839928, 4813.839928, 4813.839928, ..., 4813.839928, 4813.839928,
            4813.839928],
           [4834.570092, 4834.570092, 4834.570092, ..., 4834.570092, 4834.570092,
            4834.570092],
           [4855.386732, 4855.386732, 4855.386732, ..., 4855.386732, 4855.386732,
            4855.386732],
           ...,
           [ 538.500149,  388.589022,  622.924932, ...,  622.924932,  388.589022,
             538.500149],
           [ 536.564026,  310.070157,  568.034707, ...,  568.034707,  310.070157,
             536.564026],
           [ 400.963903,  337.167707,  505.487671, ...,  505.487671,  337.167707,
             400.963903]])

  • dxCv
    (ny, nx)
    float64
    4.824e+03 4.824e+03 ... 1.0

    long_name :

    x-distance between q-points, centered at v

    units :

    meters

    array([[4.824189e+03, 4.824189e+03, 4.824189e+03, ..., 4.824189e+03,
            4.824189e+03, 4.824189e+03],
           [4.844962e+03, 4.844962e+03, 4.844962e+03, ..., 4.844962e+03,
            4.844962e+03, 4.844962e+03],
           [4.865822e+03, 4.865822e+03, 4.865822e+03, ..., 4.865822e+03,
            4.865822e+03, 4.865822e+03],
           ...,
           [2.590796e+01, 2.461434e+03, 4.874095e+03, ..., 4.874095e+03,
            2.461434e+03, 2.590796e+01],
           [5.511444e-01, 1.978699e+03, 4.838980e+03, ..., 4.838980e+03,
            1.978699e+03, 5.511444e-01],
           [1.000000e+00, 1.075999e+03, 4.827576e+03, ..., 4.827576e+03,
            1.075999e+03, 1.000000e+00]])

  • dyCu
    (ny, nx)
    float64
    4.814e+03 4.814e+03 ... 401.0 401.0

    long_name :

    y-distance between q-points, centered at u

    units :

    meters

    array([[4813.839928, 4813.839928, 4813.839928, ..., 4813.839928, 4813.839928,
            4813.839928],
           [4834.570092, 4834.570092, 4834.570092, ..., 4834.570092, 4834.570092,
            4834.570092],
           [4855.386732, 4855.386732, 4855.386732, ..., 4855.386732, 4855.386732,
            4855.386732],
           ...,
           [ 513.028972,  496.215177,  743.936135, ...,  496.215177,  513.028972,
             539.358327],
           [ 535.864512,  424.896916,  698.217343, ...,  424.896916,  535.864512,
             536.576314],
           [ 400.964049,  335.207629,  648.168424, ...,  335.207629,  400.964049,
             400.963867]])

  • dxCu
    (ny, nx)
    float64
    4.814e+03 4.814e+03 ... 0.2751

    long_name :

    x-distance between t-points, centered at u

    units :

    meters

    array([[4.813835e+03, 4.813835e+03, 4.813835e+03, ..., 4.813835e+03,
            4.813835e+03, 4.813835e+03],
           [4.834565e+03, 4.834565e+03, 4.834565e+03, ..., 4.834565e+03,
            4.834565e+03, 4.834565e+03],
           [4.855381e+03, 4.855381e+03, 4.855381e+03, ..., 4.855381e+03,
            4.855381e+03, 4.855381e+03],
           ...,
           [1.896646e+03, 4.392170e+03, 5.661919e+03, ..., 3.126064e+03,
            8.550114e+02, 1.304876e+02],
           [1.694955e+03, 4.291331e+03, 5.637108e+03, ..., 2.923262e+03,
            6.664812e+02, 1.458182e+00],
           [1.108340e+03, 3.698618e+03, 5.624055e+03, ..., 2.332642e+03,
            9.171206e+01, 2.751040e-01]])

  • dyCv
    (ny, nx)
    float64
    4.829e+03 4.829e+03 ... 2.543e+03

    long_name :

    y-distance between t-points, centered at v

    units :

    meters

    array([[4829.387556, 4829.387556, 4829.387556, ..., 4829.387556, 4829.387556,
            4829.387556],
           [4850.182577, 4850.182577, 4850.182577, ..., 4850.182577, 4850.182577,
            4850.182577],
           [4871.064317, 4871.064317, 4871.064317, ..., 4871.064317, 4871.064317,
            4871.064317],
           ...,
           [ 537.121473,  331.082776,  581.735419, ...,  581.735419,  331.082776,
             537.121473],
           [ 401.650103,  280.120349,  521.260732, ...,  521.260732,  280.120349,
             401.650103],
           [2543.213164, 2543.391414, 2654.089183, ..., 2654.089183, 2543.391414,
            2543.213164]])

  • tarea
    (ny, nx)
    float64
    2.315e+07 2.315e+07 ... 198.9

    long_name :

    area of t-cells

    units :

    meters^2

    array([[2.314817e+07, 2.314817e+07, 2.314817e+07, ..., 2.314817e+07,
            2.314817e+07, 2.314817e+07],
           [2.334797e+07, 2.334797e+07, 2.334797e+07, ..., 2.334797e+07,
            2.334797e+07, 2.334797e+07],
           [2.354947e+07, 2.354947e+07, 2.354947e+07, ..., 2.354947e+07,
            2.354947e+07, 2.354947e+07],
           ...,
           [1.179670e+05, 1.392588e+06, 3.353762e+06, ..., 2.845335e+06,
            9.742488e+05, 6.212597e+04],
           [7.163850e+03, 1.038214e+06, 3.066370e+06, ..., 2.542306e+06,
            7.305991e+05, 4.106613e+03],
           [2.175777e+02, 6.901252e+05, 2.754726e+06, ..., 2.203130e+06,
            3.338877e+05, 1.988981e+02]])

  • tmask
    (ny, nx)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    long_name :

    ocean fraction at t-cell centers

    units :

    none

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]])

  • angle
    (ny, nx)
    float64
    -0.06156 -0.06156 ... -0.4197

    long_name :

    angle grid makes with latitude line

    units :

    degrees

    array([[-6.155584e-02, -6.155584e-02, -6.155584e-02, ..., -6.155584e-02,
            -6.155584e-02, -6.155584e-02],
           [-6.154763e-02, -6.154763e-02, -6.154763e-02, ..., -6.154763e-02,
            -6.154763e-02, -6.154763e-02],
           [-6.153935e-02, -6.153935e-02, -6.153935e-02, ..., -6.153935e-02,
            -6.153935e-02, -6.153935e-02],
           ...,
           [ 1.213988e+01,  7.541037e+01,  8.263257e+01, ..., -8.044006e+01,
            -5.694835e+01, -1.546709e+00],
           [ 4.707465e-01,  7.913764e+01,  8.538624e+01, ..., -8.380810e+01,
            -5.561404e+01, -1.398931e-01],
           [ 1.370915e+00,  8.420102e+01,  8.838555e+01, ..., -8.772906e+01,
            -4.023689e+00, -4.196843e-01]])

  • depth
    (ny, nx)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    long_name :

    depth at h points

    units :

    meters

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]])

  • ar
    (ny, nx)
    float64
    1.0 1.0 1.0 ... 0.281 1.457e+03

    long_name :

    grid aspect ratio (dyt/dxt)

    units :

    none

    array([[1.000001e+00, 1.000001e+00, 1.000001e+00, ..., 1.000001e+00,
            1.000001e+00, 1.000001e+00],
           [1.000001e+00, 1.000001e+00, 1.000001e+00, ..., 1.000001e+00,
            1.000001e+00, 1.000001e+00],
           [1.000001e+00, 1.000001e+00, 1.000001e+00, ..., 1.000001e+00,
            1.000001e+00, 1.000001e+00],
           ...,
           [4.126830e+00, 1.482502e-01, 1.272037e-01, ..., 1.272037e-01,
            1.482502e-01, 4.126830e+00],
           [3.679678e+02, 1.313869e-01, 1.170092e-01, ..., 1.170092e-01,
            1.313869e-01, 3.679678e+02],
           [1.457499e+03, 2.810253e-01, 1.046237e-01, ..., 1.046237e-01,
            2.810253e-01, 1.457499e+03]])

  • egs
    (ny, nx)
    float64
    0.0004327 0.0004327 ... 1.268e-06

    long_name :

    grid effective grid spacing

    units :

    degrees

    array([[4.326716e-04, 4.326716e-04, 4.326716e-04, ..., 4.326716e-04,
            4.326716e-04, 4.326716e-04],
           [4.345349e-04, 4.345349e-04, 4.345349e-04, ..., 4.345349e-04,
            4.345349e-04, 4.345349e-04],
           [4.364059e-04, 4.364059e-04, 4.364059e-04, ..., 4.364059e-04,
            4.364059e-04, 4.364059e-04],
           ...,
           [3.088732e-05, 1.061235e-04, 1.646896e-04, ..., 1.516934e-04,
            8.876361e-05, 2.241489e-05],
           [7.611553e-06, 9.163120e-05, 1.574753e-04, ..., 1.433884e-04,
            7.686694e-05, 5.762911e-06],
           [1.326500e-06, 7.470746e-05, 1.492585e-04, ..., 1.334812e-04,
            5.196372e-05, 1.268280e-06]])

• Attributes: (4)

  Description :

  CESM MOM6 1/4 degree grid

  Author :

  Frank, Fred, Gustavo (gmarques@ucar.edu)

  Created :

  2022-12-27T10:20:11.076318

  type :

  Glogal 1/4 degree grid file

jm, im = ds_out.lat.shape
state['LAT'] = ds_out.lat[:,int(im/2)].values
state['LON'] = ds_out.lon[int(jm/2),:].values

Overwrite salt and thetao

state['salt'].values[:] = ds_out_s_an.values[:]

state['theta0'].values[:] = ds_out_theta0.values[:]

TODO: Compare original and remapped fields

# Global attrs
state.attrs['title'] = 'surface salinity and potential temperature from WOA filled over continents'
state.attrs['src_file'] = fname
state.attrs['dst_grid_name'] = 'tx1_4'
state.attrs['author'] = 'Gustavo Marques (gmarques@ucar.edu)'
state.attrs['date'] = datetime.now().isoformat()
state.attrs['created_using'] = 'https://github.com/NCAR/tx1_4/state_restoring/state_restoring.ipynb'
# save 20200616
fname1 = 'state_restore_tx1_4_{}{}{}.nc'.format(datetime.now().isoformat()[0:4],datetime.now().isoformat()[5:7],
       datetime.now().isoformat()[8:10])
state.to_netcdf(fname1)