Mixed layer depth#

%load_ext autoreload
%autoreload 2
%%capture 
# comment above line to see details about the run(s) displayed
from misc import *
from mom6_tools.m6plot import myStats, annotateStats, xycompare
import cartopy.crs as ccrs
import cartopy.feature
import intake
%matplotlib inline
# load obs-based mld from oce-catalog
obs = 'mld-deboyer-tx2_3v2'
catalog = intake.open_catalog(diag_config_yml['oce_cat'])
print('\n Reading climatology from: ', obs)
mld_obs = catalog[obs].to_dask().where(grd_xr[0].wet > 0.)
# rename coords and and fix time
mld_obs = mld_obs.rename({'lon' : 'longitude', 'lat' : 'latitude', 'time' : 'month'})
mld_obs["month"] = np.arange(1,len(mld_obs.month)+1)
months = [0,1,2]
obs_JFM = np.ma.masked_invalid(mld_obs.mld.isel(month=months).mean('month').values)
months = [6,7,8]
obs_JAS = np.ma.masked_invalid(mld_obs.mld.isel(month=months).mean('month').values)
obs_winter = obs_JAS.copy(); obs_summer = obs_JAS.copy()
j = np.abs( grd[0].geolat[:,0] - 0. ).argmin()
obs_winter[j::,:] = obs_JFM[j::,:]
obs_summer[0:j,:] = obs_JFM[0:j,:]
 Reading climatology from:  mld-deboyer-tx2_3v2

Monthly#

months = ['January', 'February', 'March', 'April', 
          'May', 'June', 'July', 'August', 'September', 
          'October', 'November', 'December']

ds = []
for path, case, i in zip(ocn_path, casename, range(len(label))):
  ds.append(xr.open_dataset(path+case+'_MLD_monthly_clima.nc'))

January#

m=0
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/fbfffb5a3f0bd66376b7a22f782f155f16855def62eddf331ca9f63dc8a87ec8.png _images/10bcc784ac55d8a741536560a3a22236b051f9fe243ac6417f2a86f479e5a782.png

February#

m=1
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/04e724638e6a1aea3055b645877008ac43d8b0fc65b0a8e84a794cd82ccf9ac2.png _images/54a3f43d005d9b7562eb378892515d3a1445aedd76df730fe6e67c2a0be6e0f4.png

March#

m=2
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/d4ce3e2fa80ba01814e21a46b034b30a1d7ad7d1e50e87e1ca4c8b10feb5f71f.png _images/f77bb242ceabba93a76622d36e1bdd9c1de6409aaf51389c2d4369d1c887c8d3.png

April#

m=3
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/d20cdb9c7eaaf461cf8f979e5f47298e7754419ee600a62b070c0f626a4c6a22.png _images/a40442b85d754a821233af819916fae0331012ab299dd864d1bde553ddd491c1.png

May#

m=4
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/14d4fa9a5da0459821312d332066555ce887b8a48d58a17b0809a188b0b1cd7d.png _images/d84138aacaa563eb81281e2ef89694ec88313a2bea345f398b39ad2e3d0f15c0.png

June#

m=5
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).mean('xh').plot(y="yh", ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).mean('xh').plot(y="yh", ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/59cb3e70b187f6989af72e5d3e3aa61bfb1e4682f8e2ff8bd751bce0da96dbe6.png _images/05ebc9c68d03036fe15f1c77412ce55c53cc348c43a952d442b603080b1319eb.png

July#

m=6
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/4fc7c3f212fc5e1b39147068f746b15dd6bec7fd58fc88d34faa3ae62a45c6e2.png _images/11aa49c1154e6f83f1b116b20291458131e4b76aaabc6a9e65fb886dc0590af4.png

August#

m=7
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/7e6c61e7a9d9beb13998a5294e9c4fb3321873e8c8cf634b069bb5f66e7af2e0.png _images/63f90fb39597e6fa1c0d44ad3cd00a4a0e17c943b19b2891e53b189ea6896729.png

September#

m=8
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/483cd1bb86cbedbb1d9cf08540600ff272521348ce9783d57c5919707f652ddd.png _images/a16192f29c2aa2a59e2306b0a5dd78da436ea43bac8c1a7bec35baa3de7098cb.png

October#

m=9
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/276829b3a5e181429f26c2ac27974a84ffb2cd76ed4245eac11cfa0e485fb189.png _images/5cd7df139628bc58608e7cd0d1b0eb45f47362ac1d8669e9d5a26ac29379d84d.png

November#

m=10
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/6457c3f960de0c23fcde3ce1f6bb2779f62bf61cda6cb3a062d21158596c0c66.png _images/9e5ac1979442472cb33e055581fd4d8a33f95e69558571432fc2296578339f0e.png

December#

m=11
for i in range(len(label)):
    model = np.ma.masked_invalid(ds[i].mlotst.isel(month=m).values)
    obs   = np.ma.masked_where(grd[i].wet == 0, mld_obs.mld.isel(month=m).values)
    xycompare(model, 
            obs, 
            grd[i].geolon, grd[i].geolat, grd[0].areacello,
            title1 = label[i], 
            title2 = 'obs (deBoyer)', 
            suptitle=str(months[m]) +' climatology, '+ str(start_date) + ' to ' + str(end_date),
            colormap=plt.cm.nipy_spectral, dcolormap=plt.cm.bwr,
            clim = (0,1000), extend='max',
            dlim=(-200,200))
      
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    ds[i].mlotst.isel(month=m).weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])
    
mld_obs.mld.isel(month=m).weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, '+str(months[m]))
ax.grid()
ax.legend();
_images/2eed8888328c4f7a2fab19024a6beca66ab20c78d8a4834627f43de1fa5609e1.png _images/3b990dd32ab71d1614d79fc284e38494ee546a57ba7413ec4664165716112b7d.png

All months, obs#

cmap = plt.cm.nipy_spectral.copy()
cmap.set_bad(color='gray') 

g = mld_obs.mld.plot(x="longitude", y="latitude", col='month', 
            col_wrap=3,
            robust=True,
            figsize=(14,12),
            cmap=cmap,
            vmin=0., vmax=1000.,
            cbar_kwargs={"orientation": "horizontal", "pad": 0.05, "label": 'MLD (m), 0.03 density criteria'},
           )
# Add a suptitle
g.fig.suptitle('deBoyer, 2004', fontsize=16, y=1.02);
_images/7e1deaef0873bfbac375bc8f3f7174529316d41e925b197849eb82db395eaaf4.png

All months, model#

for path, case, i in zip(ocn_path, casename, range(len(casename))):
  ds =  xr.open_dataset(path+case+'_MLD_monthly_clima.nc')
  g = ds.mlotst.plot(x="longitude", y="latitude", col='month', 
                col_wrap=3,
                robust=True,
                figsize=(14,12),
                cmap=cmap,
                vmin=0., vmax=1000.,
                cbar_kwargs={"orientation": "horizontal", "pad": 0.05, "label": 'MLD (m), 0.03 density criteria'},
               )
  # Add a suptitle
  g.fig.suptitle('Run {}, mean between {} to {}'.format(label[i], start_date, end_date), fontsize=16, y=1.02);
_images/16590cf853d003fd46b479b703f1f5ea4d20b13c8a189533b9b83f7f72118efe.png

All months, model - obs#

cmap = plt.cm.bwr.copy()
cmap.set_bad(color='gray')  

for path, case, i in zip(ocn_path, casename, range(len(casename))):
  ds =  xr.open_dataset(path+case+'_MLD_monthly_clima.nc')
  diff = ds.mlotst - mld_obs.mld
  g = diff.plot(x="longitude", y="latitude", col='month', 
                col_wrap=3,
                robust=True,
                figsize=(14,12),
                cmap=cmap,
                vmin=-200, vmax=200.,
                cbar_kwargs={"orientation": "horizontal", "pad": 0.05, "label": 'MLD bias (m)'},
               )
  # Add a suptitle
  g.fig.suptitle('Run {} - Deboyer (2004)'.format(label[i]), fontsize=16, y=1.02);
_images/c0be9e9441e022268a628bac973b572bf024bc6c0b82d4fa7ec24d8990a2aec6.png

Winter#

%matplotlib inline

title = 'Mean Winter MLD, JFM(NH), JAS(SH)'
try:
  area= grd[0].area_t
except:
  area= grd[0].areacello
plot_map(obs_winter, area, grd[0], 'deBoyer', 
           vmin=0, vmax=1500, suptitle=title)
_images/7e3af9a930af38ac84fc68d3c1258acb61ddba8810106ebfa60c9d734deafc5c.png
title1 = 'Mean Winter MLD, JFM(NH), JAS(SH)'
title2 = 'Mean Winter MLD (model - obs), JFM(NH), JAS(SH)'

for path, case, i in zip(ocn_path, casename, range(len(casename))):
  da = xr.open_dataset(path+case+'_MLD_winter.nc')
  try:
    area= grd[i].area_t
  except:
    area= grd[i].areacello

  # model
  plot_map(da.MLD_winter.values, area, grd[i], label[i], 
           vmin=0, vmax=1500, suptitle=title1) 
  # model - obs
  diff = (da.MLD_winter.values - obs_winter)
  plot_map(diff, area, grd[i], label[i] + ' - deBoyer', 
           vmin=-500, vmax=500, suptitle=title2, cmap='bwr')
_images/7621070a6b128c8fa759adecd3bb1e6bd6222d69ad40313573486ff8a2ce5758.png _images/e4cf69e9a09ac2e254eca0d53829bda418178162f907f78db654bd71205f1735.png
obs_winter
masked_array(
  data=[[--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        ...,
        [--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --]],
  mask=[[ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        ...,
        [ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True]],
  fill_value=1e+20,
  dtype=float32)
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    da = xr.open_dataset(path+case+'_MLD_winter.nc')
    da.MLD_winter.weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])

obs_winter_da = xr.DataArray(
    data=obs_winter,  # New data (all zeros)
    coords=da.coords,
    dims=da.dims
)

obs_winter_da.weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, Winter JFM(NH), JAS(SH)')
ax.grid()
ax.legend();
_images/1a3305f1f5d23d3602120239ee1e0c629c4dbdf8b92ad3e3302f491e74c993aa.png

Summer#

title = 'Mean Summer MLD, JFM(SH), JAS(NH)'
try:
  area= grd[0].area_t
except:
  area= grd[0].areacello
plot_map(obs_summer, area, grd[0], 'deBoyer', vmin=0, vmax=200, 
           suptitle=title, nh='JAS', sh='JFM',)
_images/4e1f3cb323b288ba6d862d6819d10bf4425240c5b07d8d6f0e9a559bf360c561.png
title1 = 'Mean Summer MLD, JFM(SH), JAS(NH)'
title2 = 'Mean Summer MLD (model - obs), JFM(SH), JAS(NH)'

for path, case, i in zip(ocn_path, casename, range(len(casename))):
  da = xr.open_dataset(path+case+'_MLD_summer.nc')
  try:
    area= grd[i].area_t
  except:
    area= grd[i].areacello

  # model
  plot_map(da.MLD_summer.values, area, grd[i], label[i], 
           vmin=0, vmax=200, suptitle=title1, nh='JAS', sh='JFM') 
  # model - obs
  diff = (da.MLD_summer.values - obs_summer)
  plot_map(diff, area, grd[i], label[i] + ' - deBoyer', 
           vmin=-100, vmax=100, suptitle=title2, nh='JAS', 
           sh='JFM',cmap='bwr')
_images/bec278bfc254c7259c73a5e4e3d07f582260553ac1ac321d4463ef6979dae5b6.png _images/71d7aa29b3b034ddd0afe81ea61f0de3960737add18ba7057fbcbedd8d07239c.png
fig, ax = plt.subplots(figsize=(8,10))
for i in range(len(label)):
    da = xr.open_dataset(path+case+'_MLD_summer.nc')
    da.MLD_summer.weighted(grd_xr[i].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label=label[i])

obs_summer_da = xr.DataArray(
    data=obs_summer,  # New data (all zeros)
    coords=da.coords,
    dims=da.dims
)

obs_summer_da.weighted(grd_xr[0].areacello.fillna(0)).mean('xh').plot(y="yh", 
                                        ax=ax, label='obs (deBoyer)')
ax.set_title('Zonally averaged MLD, Summer JFM(SH), JAS(NH)')
ax.grid()
ax.legend();
_images/049e7a334000eeb688e2bda0ef52a115692654a2ef05a73ab392af7b8cbc7c62.png