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/09fc50880a98f43ca3a810bfed499dfc08a0cc36ac9b01829ae542ae4575c890.png _images/d73d5cea6625cd4e374d6bbf689950d4cf5490595284245e6f669852f49fca43.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/ad3a0314dfe3104a809df6277ada3285d45ce5fb2ab9eb0d94d3febd68bdf0a2.png _images/78ee31affe350b065e1bf536113953b6c6b75f3494e2a1fd108d37f9599a6bb6.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/05a594439b79144e7859f54970646d32c9426fd3a476a76ce3769a562a815950.png _images/7a1859ef413ba87dd73bf1e41e6247e8f8776c9d8376634924d8b36221b3b1d9.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/82dff0e804fa83c3addbc587ab42c1db4d9dda1af44d4439f998c3f8e4fb2474.png _images/b2801bbef222a9979436f866e08df86872610d7c8bc8e9b0c2bb390f23de35b7.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/d87fb14ad8fcafd08036560f8fbe710ae45b5c3484ded3efde793581336a69f0.png _images/b53c4ac93762a7fbee8512314f88c2255111bace87dd8e7ad4f21c2cdf177ff4.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/4f66bf9fe4e5362f66ff8978b6fa7191e870bdf4e0c84742f7619722fadd3629.png _images/f1cfac8080509f54fe16d77e659d26069b4672144d5dc64c991da580e2189c3e.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/2f8263a160ae3db8c5eff9ad3755642ea5e1b49aa6979fb76b6126cf6d4f0dcc.png _images/0e733bdfa341aa8d95e1a9ceb47baf2744bf54ea6fc767071ae7045f70434f92.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/a6093d77ad6fae1c57d5865962276215bc0b9611443011d92fc4470298090001.png _images/3b361c894d73ace0c557e45e98edb7d38e9ec6fd490448b6aece300e20b267f1.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/ed238b97481e0471ac28ae51786a0758ce1bd5a53885334767c820e84bcf0e4f.png _images/f2ccc13deeee9de8a1b568f6bcb798701e30c439b9e33bcf865195bada182b68.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/67803606592c4eeec407128700020165d8b92e562c10cda077739808054ad011.png _images/42b7b0baaebb199e704bb10309d2755acfb5e76f6bc3d2af4a08ff88865b8fbf.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/df1fbb987e65df91d617a970d1ef5ee88087ba7a2018089a9cda9581e2d636eb.png _images/bb32f476f59d565576cd8e8f4ade3d7948d28474b44b21b22d86174609c77bce.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/5a42cc9dd5e2e9c3c267c129853f196ac8b38b810ebf2ea228ab7e2152a94b3e.png _images/74dc5752b00c3a73230832b709d3325a3cbf5b4093e8cb0a27cc17dfc8b39e4a.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/04829ee08fd85a4dfa0f3c5e7aa37c2140238b00aaf28f49f2ef91dd46176e49.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/c6841731c462987992da8b1fdb3ac4471f86d945a74a3b1a8c73459dd5c6d1bb.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/51280eedfbe14678148c06b506e2a345293c11885454c6478f9aed3ad6926b7b.png _images/b06430a7f3633060a9877d95bddf461a1e0a9cd614353e6bf463fd608b49894a.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/f7690b0880f0382808d6d198c2a6edc8f291b139e70d714dfe2d6d117662265b.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/2820afd0400dab6925db8c3d25e4b70f6fdb93f665b0d133174bd8cf3a67a3a5.png _images/09e78152b40570c2e705447561c304d2cd2e244e42656136869fd633640da355.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/7cc034afcc437242dd7178fbcc751194bdcca17883c13b8f20f23f18eb2f35be.png