Add CESM and RACMO smb analyses to energy model

Instead of placing these in the SMB analysis module, it makes sense to
keep all the RACMO-CESM comparison stuff together, since these will also
be used for evaluating the energy balance in RACMO

energy/energy/cesm_racmo23_smb.py

+import os
+
+import Ngl
+import numpy as np
+import numpy.ma as ma
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+
+describe = """CESM_RACMO23_smb plot."""
+
+
+def make_plot(config=None, out_path='.',
+              racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/',
+              cism_path='/lustre/atlas1/cli115/world-shared/4ue/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in TITAN ------------------------
+    f_cism = os.path.join(cism_path, 'Greenland_5km_v1.1_SacksRev_c110629.nc')
+    f_cesm_lnd_climo_jja = os.path.join(cesm_path,
+                                        'postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_JJA_196006_200508_climo.nc')
+    f_racmo_smb_jja = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.smb.1980-1999.JJA.nc')
+    f_racmo_smb_jja_remapped = os.path.join(racmo_path,
+                                            'remapped_racmo/racmo23_GRN_monthly.smb.1980-1999.remap2cesm.JJA.nc')
+    f_racmo_mask = os.path.join(racmo_path, 'RACMO23_masks_ZGRN11.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # f_cism get the following, matrixsize[301,561]
+    # usrf(0,:,:), lat(0,:,:), lon(0,:,:)
+    # Goal: plot 3 figures, radiation overlaying on elevation usrf
+
+    # read f_cism, elevation of 5km CESM
+    # adding the [:] transform the variable to type numpy.ndarray
+    ncid0 = Dataset(f_cism)
+    usrf = ncid0.variables['usrf'][0, :, :]
+    lat = ncid0.variables['lat'][0, :, :]
+    lon = ncid0.variables['lon'][0, :, :]
+
+    # read f_cesm_lnd_climo_jja, CESM variable
+    ncid1 = Dataset(f_cesm_lnd_climo_jja)
+    smb_cesm = ncid1.variables['QICE'][0, :, :]
+    smb_cesm = smb_cesm * 3.156e7  # convert from mm/s to kg/m^2/yr
+    lat1 = ncid1.variables['lat'][:]
+    lon1 = ncid1.variables['lon'][:]
+
+    # use gris as a mask to mask smb array
+    gris_mask = ncid1.variables['gris_mask'][0, :, :]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    smb_cesm_mask = ma.masked_array(smb_cesm, mask=gris_mask.mask)
+    smb_cesm = smb_cesm_mask
+
+    # read f_racmo_smb_jja_remapped, the remapped RACMO file to calculate difference
+    ncid3 = Dataset(f_racmo_smb_jja_remapped)
+    smb_remap = ncid3.variables['smb'][0, :, :]
+    smb_remap = smb_remap
+
+    smb_remap_mask = ma.masked_array(smb_remap, mask=gris_mask.mask)
+    smb_remap = smb_remap_mask
+
+    # read f_racmo_smb_jja, RACMO variable
+    ncid2 = Dataset(f_racmo_smb_jja)
+    smb = ncid2.variables['smb'][0, :, :]
+    # lat2 = ncid2.variables['lat'][:]
+    # lon2 = ncid2.variables['lon'][:]
+    smb = smb
+
+    # read f_racmo_mask, the lat/lon for RACMO data
+    ncid4 = Dataset(f_racmo_mask)
+    lat4 = ncid4.variables['lat'][:]
+    lon4 = ncid4.variables['lon'][:]
+    racmo_elev = ncid4.variables['Topography'][0, 0, :, :]
+
+    # Use gris as a mask to mask smb array.
+    gris_mask = ncid4.variables['GrIS_mask'][:]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    smb_mask = ma.masked_array(smb, mask=gris_mask.mask)
+    smb = smb_mask
+
+    diff = smb_cesm - smb_remap
+
+    # diffmin = (np.min(diff))
+    # diffmax = (np.max(diff))
+    # print("Min ANN SMB anomaly: {}".format(diffmin))
+    # print("Max ANN SMB anomaly: {}".format(diffmax))
+
+    # ------- PLOT --------
+    #  Open a workstation for drawing the plots
+    wkres = Ngl.Resources()
+    # wkres.wkColorMap = "matlab_jet"
+    wkres.wkColorMap = "BlueWhiteOrangeRed"
+    wkres.wkOrientation = "portrait"  # "portrait" or "landscape"
+    wks_type = "png"
+    wks_img = os.path.join(out_path, "CESM_RACMO23_smb_ANN")
+    wks = Ngl.open_wks(wks_type, wks_img, wkres)
+
+    # --- for the map -------
+    # Define plotting area, Greenland
+    mres = Ngl.Resources()
+    mres.nglDraw = False  # Don't draw individual plots
+    mres.nglFrame = False  # Don't advance frame.
+
+    mres.pmTickMarkDisplayMode = "Never"  # Turn off map tickmarks.
+    mres.mpGridAndLimbOn = False  # Turn off grid and limb lines.
+    mres.mpProjection = "Aitoff"
+    mres.mpLimitMode = "LatLon"  # limit map via lat/lon, to zoom in
+    mres.mpCenterLatF = 70.  # map area
+    mres.mpCenterLonF = -44.
+    mres.mpMinLatF = 57.
+    mres.mpMaxLatF = 85.
+    mres.mpMinLonF = -55.
+    mres.mpMaxLonF = -30.
+    mres.mpOutlineOn = False
+    mres.mpFillOn = False
+    mres.mpPerimOn = True  # add box around map
+
+    # --- for the smbedo contour of CESM -------
+    res1 = Ngl.Resources()
+    res1.nglDraw = False  # Don't draw individual plots
+    res1.nglFrame = False  # Don't advance frame.
+    res1.cnLineLabelsOn = False
+    res1.cnFillOn = True
+    res1.cnLinesOn = False
+    res1.cnLineLabelsOn = False
+    res1.cnFillMode = "RasterFill"
+    res1.cnLevelSelectionMode = "ExplicitLevels"
+    res1.cnLevels = [-4000, -2000, -1000, -700, -500, -300, -200, -100, -50, -20, 20, 50, 100, 200, 300, 500, 700, 1000,
+                     2000, 5000]
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+    # res1.pmLabelBarOrthogonalPosF = -0.05
+
+    res1.sfXArray = lon1
+    res1.sfYArray = lat1
+
+    # --- for the smbedo contour of RACMO -------
+    res2 = Ngl.Resources()
+    res2.nglDraw = False  # Don't draw individual plots
+    res2.nglFrame = False  # Don't advance frame.
+    res2.cnLineLabelsOn = False
+    res2.cnFillOn = True
+    res2.cnLinesOn = False
+    res2.cnLineLabelsOn = False
+    res2.cnFillMode = "RasterFill"
+    res2.cnLevelSelectionMode = "ExplicitLevels"
+    res2.cnLevels = [-4000, -2000, -1000, -700, -500, -300, -200, -100, -50, -20, 20, 50, 100, 200, 300, 500, 700, 1000,
+                     2000, 5000]
+    # res2.cnLevels = [-12800,-6400,-3200,-1600,-800,-400,-200,-50,-25,0,25,50,100,200,400,800,1600,3200,6400]
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    # res2.pmLabelBarHeightF    = 0.5         # Change height of labelbar
+    # res2.pmLabelBarWidthF     = 0.2           # Change width of labelbar
+    # res2.pmLabelBarOrthogonalPosF = -0.02      # Move labelbar closer to plot
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res2.sfXArray = lon4
+    res2.sfYArray = lat4
+
+    # --- for the smb contour of CESM-RACMO -------
+    res3 = Ngl.Resources()
+    res3.nglDraw = False  # Don't draw individual plots
+    res3.nglFrame = False  # Don't advance frame.
+    res3.cnLineLabelsOn = False
+    res3.cnFillOn = True
+    res3.cnLinesOn = False
+    res3.cnLineLabelsOn = False
+    res3.cnFillMode = "RasterFill"
+    res3.cnLevelSelectionMode = "ExplicitLevels"
+    res3.cnLevels = [-4000, -2000, -1000, -700, -500, -300, -200, -100, -50, -20, 20, 50, 100, 200, 300, 500, 700, 1000,
+                     2000, 5000]
+
+    # res3.cnLevels = np.arange(-1000, 1000, 200)
+    res3.lbLabelBarOn = True  # Turn on labelbar.
+    res3.lbOrientation = "Vertical"  # Verticle labelbar
+    res3.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res3.sfXArray = lon1
+    res3.sfYArray = lat1
+
+    # ---- for the elevation -------
+    sres = Ngl.Resources()
+    sres.nglDraw = False  # Don't draw individual plots
+    sres.nglFrame = False  # Don't advance frame.
+    sres.cnFillOn = False
+    sres.cnLinesOn = True
+    sres.cnLineLabelsOn = False
+    sres.cnLevelSelectionMode = "ExplicitLevels"
+    sres.cnLevels = [0, 1000, 2000, 3000]
+    sres.sfXArray = lon
+    sres.sfYArray = lat
+
+    # ---- for the elevation of RACMO -------
+    sres1 = Ngl.Resources()
+    sres1.nglDraw = False  # Don't draw individual plots
+    sres1.nglFrame = False  # Don't advance frame.
+    sres1.cnFillOn = False
+    sres1.cnLinesOn = True
+    sres1.cnLineLabelsOn = False
+    sres1.cnLevelSelectionMode = "ExplicitLevels"
+    sres1.cnLevels = [0, 1000, 2000, 3000]
+    sres1.sfXArray = lon4
+    sres1.sfYArray = lat4
+
+    # ---- Overlay plots, each one has its own ID
+    # overlay ice on base map, and then overlay elevation on ice
+    # usrf is elevation 5km CESM, racmo_elev is elevation of RACMO
+    usrf_plot1 = Ngl.contour(wks, usrf, sres)
+    usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1)
+    usrf_plot3 = Ngl.contour(wks, usrf, sres)
+
+    smb_cesm_plot = Ngl.contour(wks, smb_cesm, res1)
+    smb_plot = Ngl.contour(wks, smb, res2)
+    diff_plot = Ngl.contour(wks, diff, res3)
+    # diff_plot = Ngl.contour(wks,Remap_smb,res1)
+
+    # Creat multiple figures and draw, which now contains the elevation and radiation
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM, SMB", "RACMO, SMB", "CESM-RACMO, SMB"]
+
+    nmap = 3
+    plot = []
+    for i in range(nmap):
+        mres.tiMainString = map_title[i]
+        plot.append(Ngl.map(wks, mres))
+
+    # Overlay everything on the map plot.
+    Ngl.overlay(plot[0], smb_cesm_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+
+    Ngl.overlay(plot[1], smb_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+
+    Ngl.overlay(plot[2], diff_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+
+    # Using overlay, CANNOT set panel property
+    # text_ndc should be before panel
+    # textres               = Ngl.Resources()
+    # textres.txFontHeightF = 0.02   # Size of title.
+    # Ngl.text_ndc(wks,"SMB",0.48,0.85,textres)
+
+    Ngl.panel(wks, plot, [1, nmap])
+
+    Ngl.end()
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(wks_img + '.' + wks_type))
+    img_elem = el.image('CESM_RACMO23_smb',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    return img_list
+
+
+if __name__ == '__main__':
+    make_plot()

energy/energy/timeseries_acab.py

+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+describe = """TimeSeries_T2M and BoxPlot_T2M plot."""
+
+
+def make_plot(config=None, out_path='.',
+              racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in rhea, the two files save the same data ------------------------
+    f_cesm_glc_acab_aavg = os.path.join(cesm_path, 'postproc', 'glc', 'tseries',
+                                        'b.e10.BG20TRCN.f09_g16.002.cism.h.acab_aavg.nc')
+    f_racmo_ts_smb_aavg = os.path.join(racmo_path, 'tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cesm_glc_acab_aavg
+    # time from 100-255, both time and acab have 156 data
+    ncid1 = Dataset(f_cesm_glc_acab_aavg, mode='r')
+    time = ncid1.variables['time'][:]
+    acab = ncid1.variables['acab'][:]  # m/yr
+    ncid1.set_auto_scale(True)  # put the scale-factor in the time-series plot
+
+    # acab = acab #unit transfer
+    acab = acab * 1000  # unit transfer to kg/m2/yr
+    time = time + 1751  # time from 1851-2006
+
+    # read f_racmo_ts_smb_aavg
+    ncid2 = Dataset(f_racmo_ts_smb_aavg, mode='r')
+    time2 = ncid2.variables['time'][:]
+    smb = ncid2.variables['smb'][:]  # mmWE = kg/m2/yr
+    ncid2.set_auto_scale(True)  # put the scale-factor in the time-series plot
+
+    # NOTE: add months to make the full year, divide my motnhs to get years, add to get 1958
+    time2 = (time2 + 6.5) / 12 + 1957
+
+    # print(np.max(acab))
+    # print(np.min(acab))
+    # print(np.max(time))
+    # print(np.min(time))
+    # print(np.max(smb))
+    # print(np.min(smb))
+    # print(np.max(time2))
+    # print(np.min(time2))
+
+    # ------- Time Series PLOT --------
+    fig = plt.figure(figsize=(8, 4))
+    ax = fig.add_subplot(111)
+
+    ax.plot(time, acab, 'r-')
+    ax.plot(time2, smb, 'b-')
+
+    ax.set_ylim(-10, 250)
+    ax.set_xlim(1851, 2013)
+    ax.xaxis.set_ticks(np.arange(1851, 2013, 26))
+    ax.set_xticklabels(['1851', '1877', '1903', '1929', '1955', '1981', '2006'])
+    ax.set_ylabel('SMB (kg/m2/yr)', fontsize=14)
+    ax.set_xlabel('time', fontsize=14)
+    # ax.text(-12,5,'T-2m, JJA',fontsize=18)
+    # ax.legend([s1,s2],['>20%','>99%'],loc='lower right')
+
+    ts_img_path = os.path.join(out_path, 'TimeSeries_SMB.png')
+    plt.savefig(ts_img_path)
+    plt.close()
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(ts_img_path))
+    img_elem = el.image('TimeSeries_SMB',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    # ------ Box plot of statistics ---------
+    # If I have 3 sets of data, put them a matrix with 3 columns,
+    # then boxplot(matrix) in one figure data=[a1,a2,a3]
+
+    fig = plt.figure(figsize=(4.8, 4))
+
+    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
+    # ax = fig.add_subplot(121)
+    # ax.boxplot(acab,0,'')
+    # ax.set_ylim([0,230])
+    # ax.set_ylabel('acab(kg/m2/yr)',fontsize=14)
+    # ax.set_xticklabels('')
+
+    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] with outliners
+    # --- where the Q1 is 25% percentile and Q3 is 75% percential, IQR=Q3-Q1
+    # --- the outlier is data point located outside the whiskers.
+    ax = fig.add_subplot(121)
+    ax.boxplot(acab, sym='gd')
+    ax.set_ylabel('SMB(km/m2/yr)', fontsize=14)
+    ax.set_ylim(0, 230)
+    ax.set_xticklabels('')
+    ax = fig.add_subplot(122)
+    ax.boxplot(smb, sym='gd')
+    ax.set_ylabel('SMB(km/m2/yr)', fontsize=14)
+    ax.set_ylim(0, 230)
+    ax.set_xticklabels('')
+
+    bx_img_path = os.path.join(out_path, 'BoxPlot_SMB.png')
+    plt.savefig(bx_img_path)
+    plt.close()
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(bx_img_path))
+    img_elem = el.image('BoxPlot_QICE_GT',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    return img_list
+
+
+if __name__ == '__main__':
+    plt.switch_backend('agg')
+    make_plot()

energy/energy/timeseries_qice.py

+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+describe = """TimeSeries_T2M and BoxPlot_T2M plot."""
+
+
+def make_plot(config=None, out_path='.',
+              racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in rhea, the two files save the same data ------------------------
+    f_cesm_lnd_ts_qice = os.path.join(cesm_path, 'postproc', 'lnd', 'tseries',
+                                      'b.e10.BG20TRCN.f09_g16.002.clm2.aavg.h0.yrly_ai.QICE.nc')
+    f_racmo_ts_smb_aavg = os.path.join(racmo_path, 'tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cesm_lnd_ts_qice
+    # time from 100-255, both time and acab have 156 data
+    ncid1 = Dataset(f_cesm_lnd_ts_qice, mode='r')
+    time = ncid1.variables['time'][:]
+    acab = ncid1.variables['QICE'][:]  # m/yr
+
+    # unit transfer
+    acab = acab * 3.156e7  # convert from mm/s to kg/m^2/yr
+    time = (time - time[1]) / 365 + 1851  # time from 1851-2006
+
+    # read f_racmo_ts_smb_aavg
+    ncid2 = Dataset(f_racmo_ts_smb_aavg, mode='r')
+    time2 = ncid2.variables['time'][:]
+    smb = ncid2.variables['smb'][:]  # mmWE
+
+    # NOTE: add months to make the full year, divide my motnhs to get years, add to get 1958
+    time2 = (time2 + 6.5) / 12 + 1957
+
+    # maxacab = (np.max(acab))
+    # minacab = (np.min(acab))
+    # maxtime = (np.max(time))
+    # mintime = (np.min(time))
+    # maxsmb = (np.max(smb))
+    # minsmb = (np.min(smb))
+    # maxtime2 = (np.max(time2))
+    # mintime2 = (np.min(time2))
+
+    # print("Max QICE: {}".format(maxacab))
+    # print("Min QICE: {}".format(minacab))
+    # print("Max time: {}".format(maxtime))
+    # print("Min time: {}".format(mintime))
+    # print("Max sum smb: {}".format(maxsmb))
+    # print("Min sum smb: {}".format(minsmb))
+    # print("Max time2: {}".format(maxtime2))
+    # print("Min time2: {}".format(mintime2))
+
+    # ------- Time Series PLOT --------
+    fig = plt.figure(figsize=(8, 4))
+    ax = fig.add_subplot(111)
+
+    ax.plot(time, acab, 'r-')
+    ax.plot(time2, smb, 'b-')
+
+    ax.set_ylim(0, 1200)
+    ax.set_xlim(1850, 2013)
+    ax.xaxis.set_ticks(np.arange(1850, 2013, 26))
+    ax.set_xticklabels(['1850', '1876', '1902', '1928', '1954', '1980', '2006'])
+    ax.set_ylabel('SMB (kg/m2/yr)', fontsize=14)
+    ax.set_xlabel('time', fontsize=14)
+    # ax.text(-12,5,'T-2m, JJA',fontsize=18)
+    # ax.legend([s1,s2],['>20%','>99%'],loc='lower right')
+
+    ts_img_path = os.path.join(out_path, 'TimeSeries_QICE.png')
+    plt.savefig(ts_img_path)
+    plt.close()
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(ts_img_path))
+    img_elem = el.image('TimeSeries_QICE',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    # ------ Box plot of statistics ---------
+    # If I have 3 sets of data, put them a matrix with 3 columns,
+    # then boxplot(matrix) in one figure data=[a1,a2,a3]
+
+    fig = plt.figure(figsize=(8, 4))
+
+    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
+    # ax = fig.add_subplot(121)
+    # ax.boxplot(acab,0,'')
+    # ax.set_ylim([0,230])
+    # ax.set_ylabel('acab(kg/m2/yr)',fontsize=14)
+    # ax.set_xticklabels('')
+
+    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] with outliners
+    # --- where the Q1 is 25% percentile and Q3 is 75% percential, IQR=Q3-Q1
+    # --- the outlier is data point located outside the whiskers.
+    ax = fig.add_subplot(121)
+    ax.boxplot(acab, sym='gd')
+    ax.set_ylabel('SMB(km/m2/yr)', fontsize=14)
+    ax.set_ylim(-50, 230)
+    ax.set_xticklabels('')
+    ax = fig.add_subplot(122)
+    ax.boxplot(smb, sym='gd')
+    ax.set_ylabel('SMB(km/m2/yr)', fontsize=14)
+    ax.set_ylim(-50, 230)
+    ax.set_xticklabels('')
+
+    bx_img_path = os.path.join(out_path, 'BoxPlot_QICE.png')
+    plt.savefig(bx_img_path)
+    plt.close()
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(bx_img_path))
+    img_elem = el.image('BoxPlot_QICE',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    return img_list
+
+
+if __name__ == '__main__':
+    plt.switch_backend('agg')
+    make_plot()

energy/energy/timeseries_qice_gt.py

+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+describe = """TimeSeries_T2M and BoxPlot_T2M plot."""
+
+
+def make_plot(config=None, out_path='.',
+              racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in rhea, the two files save the same data ------------------------
+    f_cesm_lnd_ts_qice = os.path.join(cesm_path, 'postproc', 'lnd', 'tseries',
+                                      'b.e10.BG20TRCN.f09_g16.002.clm2.aavg.h0.yrly_ai.QICE.nc')
+    f_racmo_ts_smb_aavg = os.path.join(racmo_path, 'tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cesm_lnd_ts_qice
+    # time from 100-255, both time and acab have 156 data
+    ncid1 = Dataset(f_cesm_lnd_ts_qice, mode='r')
+    time = ncid1.variables['time'][:]
+    acab = ncid1.variables['QICE'][:]  # m/yr
+
+    # unit transfer
+    acab = acab * 3.156e7  # convert from mm/s to kg/m^2/yr
+    acab = acab * 1.7  # Convert to Gt/yr, for validation, compare to V13
+    time = (time - time[1]) / 365 + 1851  # time from 1851-2006
+
+    # read f_racmo_ts_smb_aavg
+    ncid2 = Dataset(f_racmo_ts_smb_aavg, mode='r')
+    time2 = ncid2.variables['time'][:]
+    smb = ncid2.variables['smb'][:]  # mmWE = kg/m2
+    smb = smb * 1.7  # unit Transfer to Gt from kg/m2
+
+    # NOTE: add months to make the full year, divide my motnhs to get years, add to get 1958
+    time2 = (time2 + 6.5) / 12 + 1957
+
+    # maxacab = (np.max(acab))
+    # minacab = (np.min(acab))
+    # maxtime = (np.max(time))
+    # mintime = (np.min(time))
+    # maxsmb = (np.max(smb))
+    # minsmb = (np.min(smb))
+    # maxtime2 = (np.max(time2))
+    # mintime2 = (np.min(time2))
+
+    # print("Max QICE: {}".format(maxacab))
+    # print("Min QICE: {}".format(minacab))
+    # print("Max time: {}".format(maxtime))
+    # print("Min time: {}".format(mintime))
+    # print("Max smb Gt/yr: {}".format(maxsmb))
+    # print("Min smb Gt/yr: {}".format(minsmb))
+    # print("Max time2: {}".format(maxtime2))
+    # print("Min time2: {}".format(mintime2))
+
+    # ------- Time Series PLOT --------
+    fig = plt.figure(figsize=(8, 4))
+    ax = fig.add_subplot(111)
+
+    ax.plot(time, acab, 'r-')
+    ax.plot(time2, smb, 'b-')
+
+    ax.set_ylim(0, 700)
+    ax.set_xlim(1851, 2013)
+    ax.xaxis.set_ticks(np.arange(1851, 2013, 26))
+    ax.set_xticklabels(['1851', '1877', '1903', '1929', '1955', '1981', '2006'])
+    ax.set_ylabel('SMB (kg/m2/yr)', fontsize=14)
+    ax.set_xlabel('time', fontsize=14)
+    # ax.text(-12,5,'T-2m, JJA',fontsize=18)
+    # ax.legend([s1,s2],['>20%','>99%'],loc='lower right')
+