<jupyter_start><jupyter_text>POP 0.1deg OMIP2 Simulation: Compute UOHC/UOSC Bias on 1x1 Grid<jupyter_code>%matplotlib inline
import os
import shutil
from glob import glob
import cftime
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import intake
import cf_xarray
import pop_tools
import xesmf as xe
#import esmlab
#import util
import warnings
warnings.filterwarnings('ignore')
xr.__version__
import dask
dask.__version__<jupyter_output><empty_output><jupyter_text>Create Dask Cluster<jupyter_code># Close out Dask Cluster and release workers:
# NOTE:  only run this cell to terminate Dask Cluster!
cluster.close()
client.close()
def get_ClusterClient():
    import dask
    from dask_jobqueue import PBSCluster
    from dask.distributed import Client
    cluster = PBSCluster(
        cores=1,
        memory='20GB',
        processes=1,
        queue='casper',
        resource_spec='select=1:ncpus=1:mem=20GB',
        project='NCGD0011',
        walltime='03:00:00',
        interface='ib0',)

    dask.config.set({
        'distributed.dashboard.link':
        'https://jupyterhub.hpc.ucar.edu/stable/user/{USER}/proxy/{port}/status'
    })
    client = Client(cluster)
    return cluster, client
cluster, client = get_ClusterClient()
cluster.scale(15) 
cluster<jupyter_output><empty_output><jupyter_text>Functions to compute depth-average<jupyter_code>def limit_depth_of_field(level_bounds, field, depth_limit):
    level_bounds_limited = level_bounds.where(level_bounds < depth_limit, depth_limit)
    delta_level = abs(level_bounds_limited[:, 1] - level_bounds_limited[:, 0])
    delta_level_limited = delta_level.where(delta_level != 0, drop=True)
    field_limited = field.where(delta_level != 0, drop=True)
    return delta_level_limited, field_limited

def compute_depthavg_ds(ds,field,depth_coord,depth=100):
    """
    Compute mean over depth for field.
    """
    newvarname = field+'_'+str(depth)+'m_mean'
    dep_bnds = depth_coord+'_bnds'
    dz,data = limit_depth_of_field(ds[dep_bnds],ds[field],depth)
    data_mean = (data*dz).sum(depth_coord,skipna=False)/dz.sum(depth_coord)
    ds[newvarname] = data_mean
    ds[newvarname].attrs['units'] = data.attrs['units']
    return ds<jupyter_output><empty_output><jupyter_text>EN4 Observations<jupyter_code>fin = '/glade/p/cgd/oce/people/yeager/OBS/MetOffice.EN4/v4.2.1/monthly/EN.4.2.1.f.analysis.g10.????.nc'
ds_en4 = xr.open_mfdataset(fin,chunks={'time':12},combine='by_coords',data_vars='minimal',parallel=True)
# Create 30-year monthly climatology
ds_en4 = ds_en4[['temperature','salinity','time_bnds','depth_bnds']].sel(time=slice('1978','2007'))

# Following line is needed to circumvent this xarray issue:
#  https://github.com/pydata/xarray/issues/2145#issuecomment-391114129
ds_en4 = ds_en4.set_coords('depth_bnds')

ds_en4 = ds_en4.groupby('time.month',squeeze=True).mean('time',keep_attrs=True)
ds_en4 = ds_en4.rename({'temperature':'thetao','salinity':'so'})
ds_en4 = ds_en4.load()
# Add depth-average fields to Dataset
ds_en4['thetao'] = ds_en4['thetao'] - 273.15
ds_en4['thetao'].attrs['units'] = 'degC'
ds_en4 = compute_depthavg_ds(ds_en4,'thetao','depth',200)
ds_en4 = compute_depthavg_ds(ds_en4,'so','depth',200)
# Write to netcdf
fout = '/glade/scratch/yeager/OMDP/EN4_1978-2007clim_200m.nc'
ds_en4 = ds_en4[['thetao_200m_mean','so_200m_mean']]
ds_en4.fillna(-9999).to_netcdf(fout,encoding={'thetao_200m_mean': {'dtype': 'float', '_FillValue': -9999},'so_200m_mean': {'dtype': 'float', '_FillValue': -9999}})
ds_en4['thetao_200m_mean'].isel(month=0).plot()<jupyter_output><empty_output><jupyter_text>Regrid EN4 data to standard 1x1 grid<jupyter_code># Define target grid
woa1x1 = xr.Dataset({'lat': (['lat'], np.arange(-89.5, 90.5, 1.0)),
                     'lon': (['lon'], np.arange(0.5, 360.5, 1.0))})
regrid_en4 = xe.Regridder(ds_en4, woa1x1, method='bilinear',periodic=True)
regrid_en4  # print basic regridder information.
ds_en4_1x1 = regrid_en4(ds_en4,keep_attrs=True)
# Write to netcdf
fout = '/glade/scratch/yeager/OMDP/EN4_1978-2007clim_200m_1x1.nc'
ds_en4_1x1.to_netcdf(fout)
ds_en4_1x1<jupyter_output><empty_output><jupyter_text>Process POP 0.1deg OMIP2 data<jupyter_code>def time_set_midmonth(ds, time_name):
    year = ds[time_name].dt.year
    month = ds[time_name].dt.month
    year = xr.where(month==1,year-1,year)
    month = xr.where(month==1,12,month-1)
    nmonths = len(month)
    newtime = [cftime.DatetimeNoLeap(year[i], month[i], 15) for i in range(nmonths)]
    ds[time_name] = newtime
    return ds
%%time
dir1 = '/glade/campaign/cesm/development/omwg/g.e20.G.TL319_t13.control.001/ocn/hist/'
filetemp = dir1+'g.e20.G.TL319_t13.control.001.pop.h.????-??.nc'
infiles = sorted(glob(filetemp))
len(infiles)
def open_POPfile(file,chunks):
    ds = time_set_midmonth(xr.open_dataset(file,chunks=chunks),'time')
    ds = ds[['TEMP','SALT','z_t','z_w']].isel(z_t=slice(0,21)).isel(z_w=slice(0,22))
    return ds
infiles[-1].split('/')[-1]<jupyter_output><empty_output><jupyter_text>First, define POP0.1 to 1x1 Regridder<jupyter_code>%%time
ds = open_POPfile(infiles[0],{'z_t':4})
ds['z_t'] = ds['z_t']/100.
ds['z_t'].attrs['units'] = 'm'
ds['z_w'] = ds['z_w']/100.
ds['z_w'].attrs['units'] = 'm'
ztbnds = np.transpose(np.array([ds['z_w'].values[0:-1],ds['z_w'].values[1:]]))
tmp = xr.DataArray(ztbnds,dims=['z_t','nbnds'],coords={'z_t':ds['z_t']})
ds['z_t_bnds'] = tmp
ds = ds.rename({'TEMP':'thetao','SALT':'so'})
ds = compute_depthavg_ds(ds,'thetao','z_t',200)
ds = compute_depthavg_ds(ds,'so','z_t',200)
ds = ds.drop(['ULONG','ULAT','thetao','so'])
ds = ds.rename({'TLAT':'latitude','TLONG':'longitude'})
ds
%%time
HRPOPregridder = xe.Regridder(ds, woa1x1, method='bilinear')
ds1x1 = HRPOPregridder(ds.load())
%%time
for infile in infiles:
    cdffile = infile.split('/')[-1]
    outfile='/glade/scratch/yeager/iHESP/'+cdffile
    if (not os.path.exists(outfile)):
        print(cdffile)
        ds = open_POPfile(infile,{'z_t':4})
        ds['z_t'] = ds['z_t']/100.
        ds['z_t'].attrs['units'] = 'm'
        ds['z_w'] = ds['z_w']/100.
        ds['z_w'].attrs['units'] = 'm'
        ztbnds = np.transpose(np.array([ds['z_w'].values[0:-1],ds['z_w'].values[1:]]))
        tmp = xr.DataArray(ztbnds,dims=['z_t','nbnds'],coords={'z_t':ds['z_t']})
        ds['z_t_bnds'] = tmp
        ds = ds.rename({'TEMP':'thetao','SALT':'so'})
        ds = compute_depthavg_ds(ds,'thetao','z_t',200)
        ds = compute_depthavg_ds(ds,'so','z_t',200)
        ds = ds.drop(['ULONG','ULAT','thetao','so'])
        ds = ds.rename({'TLAT':'latitude','TLONG':'longitude'})
        ds1x1 = HRPOPregridder(ds.load())
        ds1x1.to_netcdf(outfile)<jupyter_output>g.e20.G.TL319_t13.control.001.pop.h.0001-01.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-02.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-03.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-04.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-05.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-06.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-07.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-08.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-09.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-10.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-11.nc
g.e20.G.TL319_t13.control.001.pop.h.0001-12.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-01.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-02.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-03.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-04.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-05.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-06.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-07.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-08.nc
g.e20.G.TL319_t13.control.001.pop.h.0002-09.nc
g.e20.G.TL319[...]<jupyter_text>Compute OMIP UOHC & UOSC anomalies  Read in regridded UOHC/UOSC from FOSI-HR<jupyter_code>filetemp = '/glade/scratch/yeager/iHESP/g.e20.G.TL319_t13.control.001.pop.h.????-??.nc'
ds_hr1x1 = xr.open_mfdataset(filetemp,parallel=True,combine='by_coords')
ds_hr1x1
hrt200 = ds_hr1x1.thetao_200m_mean.load()
hrs200 = ds_hr1x1.so_200m_mean.load()
hrt200 = xr.where(hrt200 > 1.e20,np.nan,hrt200)
hrs200 = xr.where(hrs200 > 1.e20,np.nan,hrs200)
%%time
ncyc=4
nyear=61
yearrange=(1958,2018)
nyr = yearrange[1]-yearrange[0]+1
climrange = (1978,2007)
entryid = 'OMIP.iHESP.CESM13.omip2'
varname = entryid+'_thetao_0-200m_bias'
attrs = {'forcing_cycle':yearrange,'ncycle':1,'nominal_resolution':'10 km'}
nmon=732
cyc_time = [cftime.DatetimeNoLeap(yearrange[0]+year, 1+month, 15) for year in range(nyr) for month in range(12)]
hrt200['time'] = cyc_time
hrs200['time'] = cyc_time
hrt200 = hrt200.sel(time=slice(str(climrange[0]),str(climrange[1])))
hrt200_clim = hrt200.groupby('time.month').mean('time',keep_attrs=True)
hrt200_climbias = hrt200_clim - ds_en4_1x1['thetao_200m_mean']
hrs200 = hrs200.sel(time=slice(str(climrange[0]),str(climrange[1])))
hrs200_clim = hrs200.groupby('time.month').mean('time',keep_attrs=True)
hrs200_climbias = hrs200_clim - ds_en4_1x1['so_200m_mean']
hrt200.groupby('time.month').mean('time',skipna=True).isel(month=1).plot()
hrt200_climbias.name = entryid+'_thetao_0-200m_bias'
hrt200_climbias = hrt200_climbias.assign_coords({'cycle':1})
hrt200_climbias = hrt200_climbias.expand_dims({'member_id':1},axis=1)
hrt200_climbias.attrs = attrs

hrs200_climbias.name = entryid+'_so_0-200m_bias'
hrs200_climbias = hrs200_climbias.assign_coords({'cycle':1})
hrs200_climbias = hrs200_climbias.expand_dims({'member_id':1},axis=1)
hrs200_climbias.attrs = attrs
ds_out = xr.Dataset({hrt200_climbias.name:hrt200_climbias,hrs200_climbias.name:hrs200_climbias,
                    hrt200_clim.name:hrt200_clim,hrs200_clim.name:hrs200_clim})
ds_out.to_netcdf('/glade/scratch/yeager/OMDP/POP0.1.OMIP2.T200.S200.climbias.nc')
ds_out
hrt200_clim<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1754226136248608, 'fraction_numerical': 0.08735701994128961, 'mean_word_length': 5.465968586387435, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '311507', 'license': 'MIT', 'n_tokens_mistral': 4806, 'n_tokens_neox': 4279, 'n_words': 505, 'path': 'NAC/POP0.1_OMIP2.thetao_so.0-200m.1x1.ipynb', 'repo_name': 'clivar-omdp/OMIP-analysis', 'size': '9879'}