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/1824c3a57ce212641edd4da6efe24496ae3ca2e49bafc44df65a19f353d11d31.png _images/411258c4d05ce16b00569ba5ce5f52c6c89800ec7e27de5e9156ce6ba57b3e51.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/15d56dc0c6b9742acfd4c0a0a42cdfd0d8ce7e11e12e54b7ed3ab9d9926e68bf.png _images/2bf00a0fe4afe09522f60170b11e05f8dbc1ebf76e3562c9bb0b2e091a4eef06.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/603e13e6e0bbcae441788098ea8e8dc2ee44f1fff366f077def8410a23437a07.png _images/bdddc404d11a2202303f064c7d60e65153a65804d01197d18e441975fec74c1e.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/9e197a1c1c5156c632fd2e5a889d5e1e3348f8b1ebeed40882afb9f5b39b739c.png _images/2381ae06b50fd5360aeddbc9a307245f6332ba4af7da806075e24e390566f7db.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/220e2177c8ce1028712d716030e9316a68cabd1d5acb170d0014a57b121c9008.png _images/4cd09d350a95195a0a8a6a2f62eeec1447ea81f401283926c8e1879d99f03ce8.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/debb481640985a94ac6db88151cb3288a0df2945dbe0612e501fb0b8d64727b7.png _images/3798a0afb50e8f036c8a7896e2c54501c53744a97dcad5a61ec5e8d868348435.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/7b5a717ebbb29d1cc16ad0ed84358a879cc8345d71edcc0feb62c93159829401.png _images/9172246ce46171a479c8ca31399e41563dad68758957fd4d45f62defc3c99245.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/6450fe32257007e92959db21468bd62ca4a61e10d2d5f3109bc2274bda6d3121.png _images/305c13a745c205c01e745d57c046b554be4fbfa959759763d887fc6c4902acaa.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/f8a574af50910844992b02d9a6a5cda86ee2a6db7f07b20bd30a47c44232ad9c.png _images/4fc2a255b8874f3f9f95428701d42eb3265264c70631b1f82a913dd6c27fe4c2.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/cc76d7397af85baeb2a73fb4d4b0b86e255a5989989fffd7b92bf02264e0b0f2.png _images/3e19d56a3acfd9f51e2dba289bcc7cccb7f056037f09ea2e5c55cab1299a7aa2.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/ea7feceed0b362eb3befccf56fba34ae154d5b636355439695771c9f81b1a476.png _images/dfd17cc513e6f082528caec98e26587708b8c09e139f570f5fc6aa91a64ab445.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/d9c269ce4c599964fa7bdfe8e11c03d9a7e2674445d4a605659acad3eb16beb6.png _images/e19eb1f52a96d6e67bc11e5f860202b99b946199a342bb50c94540c2445b7673.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/539ca1fcc437eb807dd0270c269022ae051a98abf1747655a3aa250ba1e3066d.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/c216cd73c83de2531489f794a1a5781a194f09f60e5fa98b0b2a113b14c30cc3.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/107e6fe253bfa5769d3fe72039e4e8411bf293f352cdb85f72ad4dbf75faf065.png _images/aff4f3081734f544741eed43cacf5b32c7402cf45b08ab301774020366c361c9.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/67ca3580e3e6cfec40be2c2f8a87d9bca1cdc571b5a7c844ebeebcafb12744aa.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/af41c7175e5dbfc0f6dafbf6abc8bccef0f472fe7088f42e120bc7e661676b2b.png _images/856448bd226a700bf38cf55b35b03192f4ac1bd6fcfea4b1a75b5031910807c0.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/0a93a3c7b6aa9187c5a8be2d8750e4aad407e78b9da0bfca9fee83f4f8f22e81.png