Coupled Simple Land Simulations¶
Marysa M Laguë¶
2016.09.13¶
Comparison of a coupled Simple Land/CAM5 to a coupled CLM45/CAM5 simulation¶
The spun-up coupled model we're comparing to is an old control running CLM4.5/CAM5/Slab Ocean.
The Simple Land model is built from cesm svn version: cesm1_5_beta05_mml_land
This run of the simple land model was updated to have brighter glaciers (the resulting albedo from my initial glacier albedo parameters - 0.9 vis, 0.4 nir, from the literature - ended up with an actual albedo that was considerably darker than both ceres and CLM4.5 (though CERES is mostly just measuring the SW albedo). I'm now prescribing visible albedo to be 0.95, and nir to be 0.6, which is on the high end of observed. Its still too dark/warm at high latitdues, so I've also (a) increased the snow albedo to 0.8/0.5 vis/nir from 0.7/0.4, but that run isn't finished.
Data atmosphere forced simulations of the land model¶
I've also forced the Simple Land model and a (development) version of CLM5 (r164) with the same data atmophsere to compare how just the land model behaves, without allowing any atmospheric feedbacks. Naturally it is much closer to CLM than when we couple to CAM, since the reference height temperatures and humidities are the same in the DATM forcing.
For these offline simulations, the land models are being forced with the same atmosphere, so that (in part) is why their surface temperatures are in better agreeance (than with the coupled versions). On the whole, however, SL is cooler than CLM. Also, since we're shutting off the land-atm feedback effect, nothing can get THAT out of hand with either model...
Major caveats:¶
The SL coupled simulation is run with an 1850 atmosphere, while the CLM45 coupled simulation was run at 2000 conditions - ie almost 100 ppm CO2 more, so it isn't entirely the land model's fault that CLM45 is almost a degree warmer than SL. The obvious solution to this would be to run an SL simulation with a year 2000 atmosphere, and see how it goes.
Prior to doing that, though, I'd like to sort out the glacier albedo stuff a bit more. I had initially given the glaciers vis and nir albedos from ground observations in the literature; these resulted in too dark + too warm glaciers in the simulations. I've been ramping up the albedo towards the upper end of observations (vis especially) to get closer to the ceres albedo (which really is just shortwave...) and the clm45 albedo.
# For interactive in-line plots:
#%matplotlib nbagg
# For inline plots:
%matplotlib inline
import matplotlib
import numpy as np
import os
import datetime
import netCDF4 as nc
import xarray as xr
from scipy import interpolate
import numpy.matlib
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import brewer2mpl as cbrew
import scipy.io as sio
from IPython.display import display
from IPython.display import HTML
import IPython.core.display as di # Example: di.display_html('<h3>%s:</h3>' % str, raw=True)
# This line will hide code by default when the notebook is exported as HTML
di.display_html('<script>jQuery(function() {if (jQuery("body.notebook_app").length == 0) { jQuery(".input_area").toggle(); jQuery(".prompt").toggle();}});</script>', raw=True)
# This line will add a button to toggle visibility of code blocks, for use with the HTML export version
di.display_html('''<button onclick="jQuery('.input_area').toggle(); jQuery('.prompt').toggle();">Toggle code</button>''', raw=True)
def mml_fig(LN,LT,mapdata,ds,myvar,proj,title=None,clim=None,colmap=None,units=None):
#x_data, y_data, x_label, y_label, title=None, xlim=None, ylim=None
''' Desctiption of function goes here. Also, single and double quotes are the same...'''
fig = plt.figure(figsize=(6,6))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
mp = Basemap(projection='robin',lon_0=180.,lat_0 = 0) # can't make it start anywhere other than 180???
mp.drawcoastlines()
mp.drawmapboundary(fill_color='1.') # make map background white
parallels = np.arange(-90.,90,20.)
mp.drawparallels(parallels,labels=[1,0,0,0],fontsize=10,linewidth=0.5,dashes=[1,2])
meridians = np.arange(0.,360.,20.)
mp.drawmeridians(meridians,linewidth=0.5,dashes=[1,2])
#(x, y) = m(LONXY, LATXY)
cs = mp.pcolormesh(LN,LT,mapdata,cmap=plt.cm.inferno,latlon=True)
if cm:
cs.cmap = colmap
else:
cs.cmap = plt.cm.inferno
cbar = mp.colorbar(cs,location='bottom',pad="5%")
cbar.set_label('units: '+ds[myvar].units)
if title:
plt.title(title,fontsize=12)
else:
plt.title(ds[myvar].long_name,fontsize=12)
if clim:
cbar.set_clim(clim[0],clim[1])
cs.set_clim(clim[0],clim[1])
#plt.suptitle('units?')
#plt.show()
#plt.show()
return fig, mp, ax, cbar, cs
def mml_map(LN,LT,mapdata,ds,myvar,proj,title=None,clim=None,colmap=None,cb_ttl=None):
# need to have already opened a figure/axis
#plt.sca(ax)
mp = Basemap(projection='robin',lon_0=180.,lat_0 = 0) # can't make it start anywhere other than 180???
mp.drawcoastlines()
mp.drawmapboundary(fill_color='1.') # make map background white
parallels = np.arange(-90.,90,20.)
mp.drawparallels(parallels,labels=[1,0,0,0],fontsize=10,linewidth=0.5,dashes=[1,2])
meridians = np.arange(0.,360.,20.)
mp.drawmeridians(meridians,linewidth=0.5,dashes=[1,2])
#(x, y) = m(LONXY, LATXY)
cs = mp.pcolormesh(LN,LT,mapdata,cmap=plt.cm.inferno,latlon=True)
if cm:
cs.cmap = colmap
else:
cs.cmap = plt.cm.inferno
cbar = mp.colorbar(cs,location='bottom',pad="5%")
if cb_ttl:
cbar.set_label(cb_ttl,fontsize=12)
else:
cbar.set_label('units: '+ds[myvar].units,fontsize=12)
if title:
plt.title(title,fontsize=12)
else:
plt.title(ds[myvar].long_name,fontsize=12)
if clim:
cbar.set_clim(clim[0],clim[1])
cs.set_clim(clim[0],clim[1])
#plt.suptitle('units?')
#plt.show()
#plt.show()
return mp, cbar, cs
def sl_clm_fig(myvar,myds,mapdata_1,mapdata_2,mapdata_3,ttl_1,ttl_2,ttl_3,ttl_main,clim_abs,clim_diff,units,cmap_abs,cmap_diff):
### Make a 3x3 subplot
fig, axes = plt.subplots(3, 3, figsize=(20,12))
# hide the plots below the diagonal
axes[1, 0].axis('off')
axes[2, 0].axis('off')
axes[2, 1].axis('off')
# CLM 4.5:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_1
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
#cbar0.title('units = unitless')
# SL_clm:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_2
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL_ceres:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_3
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units )
# SL_clm - CLM45:
ax4 = axes.flatten()[4]
plt.sca(ax4)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata12 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = 'Simple Land (clm '+r'$\alpha$'+') - CLM4.5'
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata12,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
# SL_cer - CLM45:
ax5 = axes.flatten()[5]
plt.sca(ax5)
mapdata_nan = mapdata_3 - mapdata_1
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata13 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = 'Simple Land (ceres '+r'$\alpha$'+') - CLM4.5'
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata13,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
# SL_cer - SL_clm:
ax8 = axes.flatten()[8]
plt.sca(ax8)
mapdata_nan = mapdata_3 - mapdata_2
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata23 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = 'Simple Land ceres '+r'$\alpha$'+' - clm '+r'$\alpha$'
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata23,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
# Figure title
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
return fig,axes
# Path to save figures:
figpath = '/Users/mlague/GoogleDrive/School/Models/SimpleLand/Analysis/figures/'
# Point at the data sets
ext_dir = '/Users/mlague/GoogleDrive/School/Models/SimpleLand/Analysis/netcdfs/'
# Coupled simulations:
clm45_clm_file = ext_dir + 'cesm13_cam5clm45_2000_ctrl.clm2.h0.20-50_year_avg.nc'
sl_clm_file = ext_dir + 'cam_albedo_ceres_ann2.clm2.h0.20-25_year_avg.nc'
clm45_cam_file = ext_dir + 'cesm13_cam5clm45_2000_ctrl.cam.h0.20-50_year_avg.nc'
sl_cam_file = ext_dir + 'cam_albedo_ceres_ann2.cam.h0.20-25_year_avg.nc'
clm45_clm_offline = ext_dir + 'datm_clm45_default.clm2.h0.10-20_year_avg.nc'
sl_clm_offline = ext_dir + 'datm_sl_albedo_ceres_ann3.clm2.h0.10-20_year_avg.nc'
# Open the coupled data sets in xarray
ds_clm45_clm = xr.open_dataset(clm45_clm_file)
ds_sl_clm = xr.open_dataset(sl_clm_file)
ds_clm45_cam = xr.open_dataset(clm45_cam_file)
ds_sl_cam = xr.open_dataset(sl_cam_file)
# Open the offline data sets (land only) in xarray
ds_clm45_off = xr.open_dataset(clm45_clm_offline)
ds_sl_off = xr.open_dataset(sl_clm_offline)
# open a cam area file produced in matlab using an EarthEllipsoid from a cam5 f19 lat/lon data set
area_f19_mat = sio.loadmat('/Users/mlague/GoogleDrive/School/Models/SimpleLand/Analysis/f19_area.mat')
area_f19 = area_f19_mat['AreaGrid']
Load some standard variables¶
lat, lon, landmask
lat = ds_sl_clm['lat'].values
lon = ds_sl_clm['lon'].values
landmask = ds_sl_clm['landmask'].values
LN,LT = np.meshgrid(lon,lat)
#print(np.shape(LN))
#print(np.shape(landmask))
Analysis¶
Make some plots of results, see how the land models compare with the same atmospheric forcing!
Coupled Simulations¶
Albedo = FSR/FSDS¶
Albedo calculated from the land surface datasets
# Plots
myvar = 'MML_fsr'
myds = ds_sl_clm
mapdata_nan = ds_clm45_clm.mean('time')['FSR'].values.squeeze()/ds_clm45_clm.mean('time')['FSDS'].values.squeeze()
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_clm.mean('time')['MML_fsr'].values.squeeze()/ds_sl_clm.mean('time')['MML_fsds'].values.squeeze()
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'unitless'
ttl_main = 'Coupled: Annual Mean Albedo = FSR/FSDS'
clim_abs = [0,1]
clim_diff = [-.25,.25]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/albedo_ceres_m_clm45.png', bbox_inches='tight')
We can compare this to what the clearsky albedo coming out of CAM is in these coupled runs:
# Plots
myvar = 'FSNSC'
myds = ds_sl_cam
mapdata_nan = (ds_clm45_cam.mean('time')['FSDSC'].values.squeeze() - ds_clm45_cam.mean('time')['FSNSC'].values.squeeze())/ds_clm45_cam.mean('time')['FSDSC'].values.squeeze()
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = (ds_sl_cam.mean('time')['FSDSC'].values.squeeze() - ds_sl_cam.mean('time')['FSNSC'].values.squeeze())/ds_sl_cam.mean('time')['FSDSC'].values.squeeze()
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'unitless'
ttl_main = 'Coupled: Annual Mean Albedo from CAM = $(FSDSC-FSNSC)/FSDSC$'
clim_abs = [0,1]
clim_diff = [-.25,.25]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/albedo_ceres_m_clm45.png', bbox_inches='tight')
As we would expect, the land albedos look pretty much the same calculated from either cam or clm (as they should - something would be very broken if they didn't!) The differences we're seeing over the oceans are a result of the active CICE model running in the coupled simulations.
The albedo over the ice sheets is still too low (too dark) in the simple land model - I think a lot of this is due to me using too low a value for nir albedo. The snow piling up on the ice sheets shouldn't have too big an effect on the albedo as I think I masked the snow albedo to have glacier albedo over the ice sheets, anyhow.
The too-low albedo at high latitues off the ice sheet is likely due to not having an appropriately tuned snow masking depth (I started with 100 kg/m2, which is roughly 25 cm of snow, then switched it to 30 kg/m2, or about 7.5 cm of snow. The results being shown here are of the first 10 years after the switch, and it definitely looks better than when it was just set to 100 kg/m2... if anything, now it is too low over hte boreal region (but still too high for the Arctic islands?) Likely snowmasking depth will be a useful tuning knob if I try to reproduce "real worl" or CLM-like results.
Temperature¶
We'll now look at the difference in 2m reference height temperature, and in surface temperature, between the two models.
First, maps of 2m air temperature:
# Plots
myvar = 'TREFHT'
myds = ds_sl_cam
mapdata_nan = (ds_clm45_cam.mean('time')['TREFHT'].values.squeeze())
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
mapdata_nan = (ds_sl_cam.mean('time')['TREFHT'].values.squeeze())
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'K'
ttl_main = 'Coupled: Annual Mean 2m Air T from CAM (TREFHT)'
clim_abs = [235,300]
clim_diff = [-7.5,7.5]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/trefht_ceres_m_clm45.png', bbox_inches='tight')
# Plots
myvar = 'TS'
myds = ds_sl_cam
mapdata_nan = (ds_clm45_cam.mean('time')['TS'].values.squeeze())
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
mapdata_nan = (ds_sl_cam.mean('time')['TS'].values.squeeze())
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'K'
ttl_main = 'Coupled: Annual Mean Surface Temperature from CAM (TS)'
clim_abs = [235,300]
clim_diff = [-7.5,7.5]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/ts_ceres_m_clm45.png', bbox_inches='tight')
We've got some pretty big temperature differences in the coupled simulation; lower down I show the results of the two land models couped to the same data atmosphere. Without the feedback between the land and the atmosphere, the differences aren't as drastic, but the simple model is being forced to behave more like clm by the prescribed reference height temperature, humidity, etc. Mostly, the simple land model is just too cold. (But not nearly as cold as when I ran it using top-of-atmosphere albedo from ceres!)
If we ignore the spatially inhomogenous temperature differences between the two simulations and just consider global average temperature, things don't look quite so concerning.
Solid lines show the surface temperature, while dashed lines show the 2m air temperature. CLM is shown in green, while the simple land model is shown in blue.
months = range(1,13)
myvar = 'MML_ts'
fig = plt.figure(figsize=(16,5))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
data1a = np.nansum(np.nansum(ds_clm45_clm['TSA'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
data2a = np.nansum(np.nansum(ds_sl_clm['MML_l2a_tref2m'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
plt.plot(months,data1a,'g--',linewidth=2.0,label='CLM45 T2m')
plt.plot(months,data2a,'b--',linewidth=2.0,label='SL T2m')
data1 = np.nansum(np.nansum(ds_clm45_clm['TV'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
data2 = np.nansum(np.nansum(ds_sl_clm[myvar].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
plt.plot(months,data1,'g',linewidth=2.0,label='CLM45 TV')
plt.plot(months,data2,'b',linewidth=2.0,label='SL TS')
plt.legend(fontsize=18,loc='best')
plt.xlim([1,12])
plt.xticks(months,fontsize=18)
plt.yticks(fontsize=18)
plt.grid()
plt.xlabel('Month',fontsize=18)
plt.ylabel('Global Avg T [K]',fontsize=18)
plt.title('Global Avg Temperature (Coupled Simulations, Land Only)',
fontsize=24, y=1.08)
#plt.ylim([98,104.5])
plt.savefig(figpath+'/cam_runs/lnd_TS_T2m_line.pdf', bbox_inches='tight')
Weird thing to note here: $\Delta$ TS-T2m in the Simple Land model is much smaller than that in CLM4.5. This would suggest the MO coupling is behaving differently in the two setups. The Simple Land is setup right now with a globally uniform "roughness" of 10m, which I would have thought was unrealistically high for most places... a high roughness could lead to lots of turbulent mixing which you could argue then leads to similar temperatures between the sfc and 2m (from a physics point of view - I should go check the MO equations and refresh my memory on how the roughness length $z_{0m}$ comes into play with the logarithmic temperature profile).
months = range(1,13)
myvar = 'MML_ts'
fig = plt.figure(figsize=(16,5))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
data1a = np.nansum(np.nansum(ds_clm45_cam['TREFHT'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19)
data2a = np.nansum(np.nansum(ds_sl_cam['TREFHT'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19)
plt.plot(months,data1a,'g--',linewidth=2.0,label='CLM45 T2m')
plt.plot(months,data2a,'b--',linewidth=2.0,label='SL T2m')
data1 = np.nansum(np.nansum(ds_clm45_cam['TS'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19)
data2 = np.nansum(np.nansum(ds_sl_cam['TS'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19)
plt.plot(months,data1,'g',linewidth=2.0,label='CLM45 TS')
plt.plot(months,data2,'b',linewidth=2.0,label='SL TS')
plt.legend(fontsize=18,loc='best')
plt.xlim([1,12])
plt.xticks(months,fontsize=18)
plt.yticks(fontsize=18)
plt.grid()
plt.xlabel('Month',fontsize=18)
plt.ylabel('Global Avg T [K]',fontsize=18)
plt.title('Global Avg Temperature (Coupled Simulations, Whole Globe)',
fontsize=24, y=1.08)
#plt.ylim([98,104.5])
plt.savefig(figpath+'/cam_runs/atm_TS_T2m_line.pdf', bbox_inches='tight')
The simple land model is a bit colder than CLM, even though the ice sheets are still a bit too dark. The differnce isn't as large over land, but its still there.
Top of Atmosphere Energy Budget?¶
Here I'll check the globally integrated TOA energy imbalance for both models, calculated as
$\text{TOA_imbalance} = \text{area} \times (\text{FSNT} - \text{FLNT}) \text{ PetaWatts}$
or
$\text{TOA_imbalance} = \text{area} \times (\text{FSNT} - \text{FLNT})/\text{sum(area)} \text{ W/m}^2$
TOA_imb_clm45_pw = 10**-15*np.sum(area_f19*(ds_clm45_cam.mean('time')['FSNT'].values.squeeze()-ds_clm45_cam.mean('time')['FLNT'].values.squeeze()))
TOA_imb_sl_pw = 10**-15*np.sum(area_f19*(ds_sl_cam.mean('time')['FSNT'].values.squeeze()-ds_sl_cam.mean('time')['FLNT'].values.squeeze()))
TOA_imb_clm45_wm2 = np.sum(area_f19*(ds_clm45_cam.mean('time')['FSNT'].values.squeeze()-ds_clm45_cam.mean('time')['FLNT'].values.squeeze()))/np.sum(area_f19)
TOA_imb_sl_wm2 = np.sum(area_f19*(ds_sl_cam.mean('time')['FSNT'].values.squeeze()-ds_sl_cam.mean('time')['FLNT'].values.squeeze()))/np.sum(area_f19)
print('TOA Energy Imbalance [PetaWatts]:')
print('CLM45 = ',TOA_imb_clm45_pw)
print('Simple Land = ',TOA_imb_sl_pw)
print('\n TOA Energy Imbalance [W/m2]:')
print('CLM45 = ',TOA_imb_clm45_wm2)
print('Simple Land = ',TOA_imb_sl_wm2)
Okay, so the TOA energy balance is off, about 2-2.5 times worse than with CLM. On the bright side, its in the same ballpark as the CLM imbalance, so we didn't break things too badly. Of course, it would be nice if the energy budget were just closed regardless of what the surface looked like (ie that the surface and atm would just adjust to close the system).
FSNT is positive down, FLNT is positive up, so if TOA < 0, FLNT > FSNT => Earth is losing heat over the long term ... interesting... when yellowstone comes back up, make a time series of the global mean temperature and see if its trending downwards. Also check what my CO$_2$ concentration is - I think its built off an 1850 run, so it'd be at 287.5, which would right there explain why it is colder (coupled) than the CLM4.5 simulation, which was using 355 ppm, aiming for a ~ yr 2000 simulation. So these aren't direct comparisons, just something to give an idea of where the biggest issues are likely to be.
Turbulent fluxes¶
We can also take a look at the senible and latent heat fluxes. To keep things easier to look at, I'm just going to plot them over land, for now.
# Plots
myvar = 'MML_lhflx'
myds = ds_sl_clm
#mapdata_1 = ds_clm45.mean('time')['TSA'].values
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
#mapdata_3 = ds_sl_cer.mean('time')[myvar].values
mapdata_nan = ds_clm45_clm.mean('time')['EFLX_LH_TOT'].values
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_clm.mean('time')[myvar].values
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'W/m2'
ttl_main = 'Annual Mean Latent Heat Flux'
clim_abs = [0, 100]
clim_diff = [-30, 30]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/lhflx_ceres_m_clm45.png', bbox_inches='tight')
The general pattern of latent heat flux looks approximately right. Looking at the difference plot, we see that the latent heat fluxes are the most off over Indonesia & Australia, with Australia not fluxing enough water (I think its bucket is empty... see lower down... probably it was too easy to evaporate water from it, then it never gets filled with rain?), and Indonesia fluxing way too much.
I should note, however, that I'm still running with globally uniform surface parameters for everything except albedo in these simulations. The roughness height = 10m everywhere, the snowmasking depth is ~25 cm (100 kg/m2, which is probably too high), the bucket capacity = 200 kg/m2, and the evaporative resistance = 20 s/m.
# Plots
myvar = 'MML_shflx'
myds = ds_sl_clm
#mapdata_1 = ds_clm45.mean('time')['TSA'].values
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
#mapdata_3 = ds_sl_cer.mean('time')[myvar].values
mapdata_nan = ds_clm45_clm.mean('time')['FSH'].values
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_clm.mean('time')[myvar].values
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'W/m2'
ttl_main = 'Annual Mean Sensible Heat Flux'
clim_abs = [0, 100]
clim_diff = [-30,30]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/shflx_ceres_m_clm45.png', bbox_inches='tight')
The sensible heat fluxes also have some pretty drastic differences, particularly in the mid-latitudes and tropics (thoug the Amazon isn't so bad...) The simple land model is getting particularly high sensible heat fluxes over northern Africa and south/central Asia.
Water¶
Here I take a quick look at how much snow and water are in their respective buckets in the annual mean.
# Plots
myvar = 'MML_snow'
myds = ds_sl_clm
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_clm.mean('time')[myvar].values
mapdata_3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_3 = 'Simple Land'
units = 'W/m2'
ttl_main = 'Annual Mean Snow Depth (mm = kg/m2)'
clim_abs = [0, 80]
clim_diff = [-30, 30]
cmap_abs = plt.cm.Blues
cmap_diff = plt.cm.RdBu_r
#sl_clm_fig(myvar,myds,mapdata_1,mapdata_2,mapdata_3,ttl_1,ttl_2,ttl_3,ttl_main,clim_abs,clim_diff,units,cmap_abs,cmap_diff)
#plt.savefig(figpath+'/cam_runs/T2m_9maps.pdf', bbox_inches='tight')
#plt.savefig(figpath+'/cam_runs/T2m_9maps.png', bbox_inches='tight')
fig, ax = plt.subplots(1, 1, figsize=(8,5))
mml_map(LN,LT,mapdata_3,myds,myvar,'moll',title='Simple Land',clim=clim_abs,colmap=cmap_abs,cb_ttl='units = '+units)
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/snow_sl.png', bbox_inches='tight')
Okay, so a snowmasking depth of 100 kg/m2 is probably way too big (about 25 cm, see below), which is what's keeping my high latitude albedo too dark. I've modified it to 30 kg/m2, or 7.5 cm of snow to go to full snow albedo... which will probably be to bright, but it could be a useful tuning tool. Looking more closely at how much snow ends up in the high-latitude snow bucket during winter months, and looking at that along with winter albedo, should give some indication of if the snow masking depth is the main cause of our high-latitude albedo problems.
Lots of snow piling up over the ice sheets is fine, as I haven't got a cap/limit on the snow bucket capacity. If we ever couple this to a dynamic ocean, we'd need to impose something like that, otherwise the ocean will get super salty because all the fresh water will just get stuck in snow on the ice sheets. In the current implementation though, the water budget isn't closed, and once snow exceeds the snow masking depth, adding more snow has no additional effect on albedo (or anything else, except that its there and available if you ever tried to melt it).
Snow density/depth calculations: Assuming the snow has a density of 400 kg/m3 (X kg/m2 / 400 kg/m3 = X/400 m = X/4 cm). I got that guess at snow density (which ranges from 100 kg/m3 for hoar (dusty light snow) to 800 kg/m3 for wet snow / packed firn) from http://sciencelearn.org.nz/Contexts/Icy-Ecosystems/Looking-closer/Snow-and-ice-density, which in turn got it from [Paterson, W.S.B. 1994. The Physics of Glaciers.] (which I still need to look up).
# Plots
myvar = 'MML_water'
myds = ds_sl_clm
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_clm.mean('time')[myvar].values
mapdata_3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_3 = 'Simple Land'
units = 'W/m2'
ttl_main = 'Annual Mean Bucket Water (mm = kg/m2)'
clim_abs = [0, 200]
clim_diff = [-30, 30]
cmap_abs = plt.cm.Blues
cmap_diff = plt.cm.RdBu_r
#sl_clm_fig(myvar,myds,mapdata_1,mapdata_2,mapdata_3,ttl_1,ttl_2,ttl_3,ttl_main,clim_abs,clim_diff,units,cmap_abs,cmap_diff)
#plt.savefig(figpath+'/cam_runs/T2m_9maps.pdf', bbox_inches='tight')
#plt.savefig(figpath+'/cam_runs/T2m_9maps.png', bbox_inches='tight')
fig, ax = plt.subplots(1, 1, figsize=(8,5))
mml_map(LN,LT,mapdata_3,ds_sl_clm,myvar,'moll',title='Simple Land',clim=clim_abs,colmap=cmap_abs,cb_ttl='units = '+units)
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/water_sl.png', bbox_inches='tight')
The water bucket looks pretty full, esp at high latitudes... I might want to consider initializing it with less water (eg 10 kg in 200 kg bucket, instead of 150)? Would it fill up to these levels, or is the bucket just super full here because its got startup-water in it that never evaporated? On the bright side, at least dry places like the Australian Outback and the Sahara Desert are looking nice and empty. The buckets are "full" looking under the ice sheets because I initialized them with 150 kg/m2 of water in them, but since they're frozen they haven't been able to evaporate much of that water out (though sublimation is allowed... but it would come from snow first, and snow is piling up in large quantities).
Atmospheric variables¶
In particular, how do the cloud fields look? Thanks to Bill Sacks, I've got the coupler data from the CLM45-CAM5 run, which means the simple land model is now giving the atmosphere dust fluxes equal each month to the monthly mean dust flux that came out of the CLM45-CAM5 run.
# Plots
myvar = 'CLDLOW'
myds = ds_sl_cam
mapdata_nan = (ds_clm45_cam.mean('time')['CLDLOW'].values.squeeze())
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
mapdata_nan = (ds_sl_cam.mean('time')['CLDLOW'].values.squeeze())
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'Cloud Fraction'
ttl_main = 'Coupled: Annual Mean Low Cloud Fraction from CAM (CLDLOW)'
clim_abs = [0,1]
clim_diff = [-0.25,.25]
cmap_abs = plt.cm.Blues
cmap_diff = plt.cm.RdBu
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/cldlow_ceres_m_clm45.png', bbox_inches='tight')
# Plots
myvar = 'CLDTOT'
myds = ds_sl_cam
mapdata_nan = (ds_clm45_cam.mean('time')['CLDTOT'].values.squeeze())
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
mapdata_nan = (ds_sl_cam.mean('time')['CLDTOT'].values.squeeze())
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'Cloud Fraction'
ttl_main = 'Coupled: Annual Mean Total Cloud Fraction from CAM (CLDLOW)'
clim_abs = [0,1]
clim_diff = [-0.25,.25]
cmap_abs = plt.cm.Blues
cmap_diff = plt.cm.RdBu
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/cldtot_ceres_m_clm45.png', bbox_inches='tight')
Over land, the clouds look particularly dense over the southeastern USA, and lacking over Australia.
The cloud decks over the East Pacific are weaker in the simple land model... over South America this would suggest the air coming off the continent is too dry (hot? or just lacking water? check!); I need to refresh my memory over what circulation drives the one off Mexico/California, but it should be roughly the same idea (trade winds).
I am pretty sure the version of CAM that the Simple Land model is coupled to was also under developement, so its possible that the cloud-deck changes are a result of changes in CAM. To validate this, I should run a "control" simulation building CLM4.5 (its still in the source code) with CAM5 (from my checked out version of cesm that the simple land is built around), and compare that to the SL-CAM simulation.
Offline Simulations¶
(Forced with the same data atmosphere)
Here I'll just repeat a couple of the plots to show how the simple land surface acts when its forced with a data atmosphere identical to that forcing a (development version of) CLM5.
# Plots
myvar = 'TS'
myds = ds_sl_off
mapdata_nan = (ds_clm45_off.mean('time')['TV'].values.squeeze())
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
mapdata_nan = (ds_sl_off.mean('time')['MML_ts'].values.squeeze())
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land (ts)'
ttl_2 = 'CLM4.5 (tv)'
ttl_3 = 'Simple Land - CLM4.5'
units = 'K'
ttl_main = 'Offline: Annual Mean Surface Temperature (from land model)'
clim_abs = [235,300]
clim_diff = [-7.5,7.5]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/ts_ceres_m_clm45_offline.png', bbox_inches='tight')
When forced with the same atmosphere, the simple land surface temperature is cooler than that in clm4.5. Caveat: I'm looking at vegetation temperature in CLM4.5, since a "surface skin temperature" isn't actually output - I should go take a peak to see what clm hands the coupler as for ts. Of course, the simple land may simply get lower temperatures than CLM - its probably too easy to evaporate with SL, especially in the tropics.
months = range(1,13)
myvar = 'MML_ts'
fig = plt.figure(figsize=(16,5))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
data1a = np.nansum(np.nansum(ds_clm45_off['TSA'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
data2a = np.nansum(np.nansum(ds_sl_off['MML_l2a_tref2m'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
plt.plot(months,data1a,'g--',linewidth=2.0,label='CLM45 T2m')
plt.plot(months,data2a,'b--',linewidth=2.0,label='SL T2m')
data1 = np.nansum(np.nansum(ds_clm45_off['TV'].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
data2 = np.nansum(np.nansum(ds_sl_off[myvar].values*area_f19,axis=1,keepdims=True),axis=2,keepdims=True).squeeze()/np.nansum(area_f19*landmask)
plt.plot(months,data1,'g',linewidth=2.0,label='CLM45 TV')
plt.plot(months,data2,'b',linewidth=2.0,label='SL TS')
plt.legend(fontsize=18,loc='best')
plt.xlim([1,12])
plt.xticks(months,fontsize=18)
plt.yticks(fontsize=18)
plt.grid()
plt.xlabel('Month',fontsize=18)
plt.ylabel('Global Avg T [K]',fontsize=18)
plt.title('Global Avg Temperature (Offline Simulations, Land Only)',
fontsize=24, y=1.08)
#plt.ylim([98,104.5])
plt.savefig(figpath+'/cam_runs/lnd_TS_T2m_line_offline.pdf', bbox_inches='tight')
The 2m air temperatures in the offline versions are almost identical, given they're using MO theory and the same refernece height temperature (at 30m, I think) to get the 2m T. The surface T is still about 0.5-1 K cooler in the simple model than in CLM.
Weirdest thing here, as in the coupld simulations, is how close the simple land skin temperature is to the 2mT. I've used a large roughness length everywhere, which could be the cause of this, but I need to think myself through that more carefully to make sure I didn't just convince myself of something that should actually be acting in the opposite direction... (argument: high roughness = lots of mixing = sfc and 2m temperatures more similar).
# Plots
myvar = 'MML_lhflx'
myds = ds_sl_off
#mapdata_1 = ds_clm45.mean('time')['TSA'].values
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
#mapdata_3 = ds_sl_cer.mean('time')[myvar].values
mapdata_nan = ds_clm45_off.mean('time')['EFLX_LH_TOT'].values
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_off.mean('time')[myvar].values
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'W/m2'
ttl_main = 'Annual Mean Latent Heat Flux, Offline'
clim_abs = [0, 100]
clim_diff = [-30, 30]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/lhflx_ceres_m_clm45_offline.png', bbox_inches='tight')
Latent heat fluxes are MUCH stronger in the simple land model, when forced with the data atmosphere, suggesting that it is too easy to get water out of the ground in the current model configuration.
# Plots
myvar = 'MML_shflx'
myds = ds_sl_off
#mapdata_1 = ds_clm45.mean('time')['TSA'].values
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
#mapdata_3 = ds_sl_cer.mean('time')[myvar].values
mapdata_nan = ds_clm45_off.mean('time')['FSH'].values
mapdata_clm = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
#mapdata_2 = np.nan*np.zeros(np.shape(albedo_clm45))
mapdata_nan = ds_sl_off.mean('time')[myvar].values
mapdata_sl = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl_1 = 'Simple Land'
ttl_2 = 'CLM4.5'
ttl_3 = 'Simple Land - CLM4.5'
units = 'W/m2'
ttl_main = 'Annual Mean Sensible Heat Flux, Offline'
clim_abs = [0, 100]
clim_diff = [-30,30]
cmap_abs = plt.cm.inferno
cmap_diff = plt.cm.RdBu_r
fig, axes = plt.subplots(1, 3, figsize=(20,6))
# SL:
ax0 = axes.flatten()[0]
plt.sca(ax0)
mapdata_nan = mapdata_sl
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata1 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_1
mp0, cbar0, cs0 = mml_map(LN,LT,mapdata1,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# CLM:
ax1 = axes.flatten()[1]
plt.sca(ax1)
mapdata_nan = mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata2 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_2
mp2, cbar2, cs2 = mml_map(LN,LT,mapdata2,myds,myvar,'moll',title=ttl,clim=clim_abs,colmap=cmap_abs, cb_ttl='units: '+units ) #plt.cm.BuPu_r
# SL - CLM:
ax2 = axes.flatten()[2]
plt.sca(ax2)
mapdata_nan = mapdata_sl - mapdata_clm
#create masked array where nan=mask. pcolormesh does not like NaNs.
mapdata3 = np.ma.masked_where(np.isnan(mapdata_nan),mapdata_nan)
ttl = ttl_3
mp3, cbar3, cs3 = mml_map(LN,LT,mapdata3,myds,myvar,'moll',title=ttl,clim=clim_diff,colmap=cmap_diff, cb_ttl='units: '+units )
fig.subplots_adjust(top=0.95)
fig.suptitle(ttl_main,
fontsize=18)
plt.savefig(figpath+'/cam_runs/shflx_ceres_m_clm45_offline.png', bbox_inches='tight')
Sensible heat fluxes are a lot weaker in the simple land model, but I expect this is largely related to how easy it is for turbulent fluxes to go through the latent heat pathway. Definitely something to dig in to!