Heat and fresh water fluxes¶

Created by Ivan Lima on Mon Jul 18 2022 12:41:39 -0400

In [1]:
import xarray as xr
import pandas as pd
import numpy as np
import holoviews as hv
import hvplot.xarray, hvplot.pandas
import os, datetime, warnings, my_bokeh_themes
from cesm_utils import fixtime, spd
warnings.filterwarnings('ignore')
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Fri Aug 12 10:20:24 2022
In [2]:
from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851

case_list = ['cesm2.B1850.f09_g17.1dpop.011', 'cesm2.B1850.f09_g17.1dpop.013',
             'cesm2.B1850.f09_g17.1dpop.test01', 'cesm2.B1850.f09_g17.1dpop.test02', 'cesm2.B1850.f09_g17.1dpop.test03']
# case_list = ['cesm2.B1850.f09_g17.1dpop.010', 'cesm2.B1850.f09_g17.1dpop.011', 'cesm2.B1850.f09_g17.1dpop.012',
#              'cesm2.B1850.f09_g17.1dpop.013', 'cesm2.B1850.f09_g17.1dpop.014']

def global_mean(ds, vname):
    weights = ds.TAREA.where(ds[vname].isel(time=0).notnull())
    da = ds[vname].weighted(weights.fillna(0)).mean(dim=['nlat','nlon'])
    da.attrs = ds[vname].attrs
    da.time.attrs = ds.time.attrs
    return da

units = {}
latent_heat_fusion_mks = 333700
def get_data():
    varlist = ['SHF', 'SFWF', 'QFLUX', 'MELTH_F', 'MELT_F']
    case_dict = {}
    data_dict = {}
    for case in case_list:
        for vname in varlist:
            with xr.open_dataset(os.path.join(datadir, '{}.{}.annual.nc'.format(case, vname)), decode_times=False, chunks={'time':10}) as ds:
                ds = fixtime(ds, start_yr)
                ds['time'] = ds.time.to_index().to_datetimeindex().year
                units[vname] = ds[vname].units
                if case in ['cesm2.B1850.f09_g17.1dpop.011', 'cesm2.B1850.f09_g17.1dpop.013']:
                    ds = ds.isel(time=slice(None, 30)) # just first 30 years
                data_dict[vname] = global_mean(ds, vname).to_series()

            df = pd.DataFrame(data_dict)

        df['SHF+QFLUX'] = df.SHF + df.QFLUX
        df['-QFLUX/L_f'] = - df.QFLUX/latent_heat_fusion_mks
        df['SFWF-QFLUX/L_f'] = df.SFWF - df.QFLUX/latent_heat_fusion_mks
        df['MELTH_F+QFLUX'] = df.MELTH_F + df.QFLUX
        df['MELT_F-QFLUX/L_f'] = df.MELT_F - df.QFLUX/latent_heat_fusion_mks
        df['SHF-MELTH_F'] = df.SHF - df.MELTH_F
        df['SFWF-MELT_F'] = df.SFWF - df.MELT_F
        
        units['-QFLUX/L_f'] = units['SFWF']
        units['SHF+QFLUX'] = units['SHF']
        units['SFWF-QFLUX/L_f'] = units['SFWF']
        units['MELTH_F+QFLUX'] = units['MELTH_F']
        units['MELT_F-QFLUX/L_f'] = units['MELT_F']
        units['SHF-MELTH_F'] = units['SHF']
        units['SFWF-MELT_F'] = units['SFWF']
        
        case_dict[case] = df
    return case_dict

def plot_vars(varlist, ncols=2):
    case_labels = ['.'.join(case.split('.')[-2:]) for case in case_list]
    plist = hv.Layout()
    for vname in varlist:
        df = pd.concat([dd[case][vname] for case in case_list], axis=1, keys=case_labels)
        p = df.hvplot(aspect=1.1, frame_height=350, grid=True, title='Global annual mean {}'.format(vname), ylabel=units[vname])
        plist += p
    
    return plist.cols(ncols).opts(shared_axes=False)

dd = get_data()

for case in sorted(case_list):
    print('{}: {}'.format(case, case_info[case]))

cesm2.B1850.f09_g17.1dpop.011: 1D + G-terms (corrected) + SSS restoring (tau = 30 days) + EBM off + start year = 1161
cesm2.B1850.f09_g17.1dpop.013: 1D + G-terms (corrected) + SSS restoring (tau = 30 days) + EBM off + start year = 6 restoring run
cesm2.B1850.f09_g17.1dpop.test01: start = yr 1161 from PI control, gterm = 1-20 mean from 1dpop.010
cesm2.B1850.f09_g17.1dpop.test02: start = yr 1161 from PI control, gterm = 6-20 mean from 1dpop.010
cesm2.B1850.f09_g17.1dpop.test03: start = yr 6 from 1dpop.010, gterm = 6-20 mean from 1dpop.010
In [3]:
plot_vars(['SHF','SFWF','QFLUX','-QFLUX/L_f','SHF+QFLUX','SFWF-QFLUX/L_f'])
Out[3]:
In [4]:
plot_vars(['MELTH_F','MELT_F','MELTH_F+QFLUX','MELT_F-QFLUX/L_f','SHF-MELTH_F','SFWF-MELT_F'])
Out[4]:

Ice volume¶

In [7]:
def fix_ice_time(ds_in,yr0=1):
    ds_out = ds_in.assign(time = ds_in.time.values)
    ds_out['time'] = ds_out.time.assign_attrs(**ds_in.time.attrs)
    ds_out['time'].attrs['units'] = ds_in.time.units.replace('0001',str(yr0))
    if 'time_bound' in ds_out:
        ds_out = ds_out.drop('time_bound')

    return xr.decode_cf(ds_out)

def get_ice_volume(nh=True):
    data_dict = {}
    for case in case_list:
        col_name = '.'.join(case.split('.')[-2:])
        with xr.open_dataset(os.path.join(datadir, '{}.cice.annual.nc'.format(case)), decode_times=False, chunks={'time':10}) as ds:
            ds = fix_ice_time(ds, start_yr)
            ds['time'] = ds.time.to_index().to_datetimeindex().year
            if case in ['cesm2.B1850.f09_g17.1dpop.011', 'cesm2.B1850.f09_g17.1dpop.013']:
                ds = ds.isel(time=slice(None, 30)) # just first 30 years
            ice_vol = ds.aice * ds.tarea * ds.hi
            if nh: # northern hemisphere
                data_dict[col_name] = ice_vol.where(ds.TLAT >=0).sum(dim=['nj','ni']).to_series() * 1.e-9 # m^3 -> Km^2
            else: # southern hemisphere
                data_dict[col_name] = ice_vol.where(ds.TLAT <0).sum(dim=['nj','ni']).to_series() * 1.e-9 # m^3 -> Km^2

    return pd.DataFrame(data_dict)

df_ice_vol_nh = get_ice_volume(nh=True)
df_ice_vol_sh = get_ice_volume(nh=False)
In [8]:
p1 = df_ice_vol_nh.hvplot(aspect=1.1, frame_height=350, grid=True, title='Global ice volume - northern hemisphere', ylabel='Km^3')
p2 = df_ice_vol_sh.hvplot(aspect=1.1, frame_height=350, grid=True, title='Global ice volume - southern hemisphere', ylabel='Km^3')
(p1 + p2).cols(2).opts(shared_axes=False)
Out[8]: