Adding the energy balance extension, which should now work on Rhea (only)

energy/energy/cesm_racmo23_albedo.py

+import os
+import Ngl
+import numpy.ma as ma
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+
+describe = """CESM_RACMO23_albedo 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_alb_jja = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.alb.1980-1999.JJA.nc')
+    f_racmo_alb_jja_remapped = os.path.join(racmo_path,
+                                            'remapped_racmo/racmo23_GRN_monthly.alb.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, the elevation data for 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, variables of CESM, lat/lon in file1 are 1D vector
+    ncid1 = Dataset(f_cesm_lnd_climo_jja)
+    fsds = ncid1.variables['FSDS'][0, :, :]
+    fsr = ncid1.variables['FSR'][0, :, :]
+    lat1 = ncid1.variables['lat'][:]
+    lon1 = ncid1.variables['lon'][:]
+
+    albedo = fsr / fsds
+
+    # use gris as a mask to mask albedo array
+    gris_mask = ncid1.variables['gris_mask'][0, :, :]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    albedo_mask = ma.masked_array(albedo, mask=gris_mask.mask)
+    albedo = albedo_mask
+
+    # read f_racmo_alb_jja, variables of RACMO
+    ncid2 = Dataset(f_racmo_alb_jja)
+    alb = ncid2.variables['alb'][0, :, :]
+
+    # 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 albedo array.
+    gris_mask = ncid4.variables['GrIS_mask'][:]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    alb_mask = ma.masked_array(alb, mask=gris_mask.mask)
+    alb = alb_mask
+
+    # read input_file3, the remapped RACMO file to calculate difference
+    ncid3 = Dataset(f_racmo_alb_jja_remapped)
+    remap_alb = ncid3.variables['alb'][0, :, :]
+
+    diff = albedo - remap_alb
+
+    # print(np.max(diff))
+    # print(np.min(diff))
+
+    # ------- 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_albedo_JJA")
+    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 albedo 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 = [0.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+
+    res1.sfXArray = lon1
+    res1.sfYArray = lat1
+
+    # --- for the albedo 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 = [0.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res2.sfXArray = lon4
+    res2.sfYArray = lat4
+
+    # --- for the albedo 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 = [-0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
+    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 for 5km CESM; racmo_elev is for RACMO
+    usrf_plot1 = Ngl.contour(wks, usrf, sres)
+    usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1)
+    usrf_plot3 = Ngl.contour(wks, usrf, sres)
+
+    albedo_plot = Ngl.contour(wks, albedo, res1)
+    alb_plot = Ngl.contour(wks, alb, res2)
+    # diff_plot = Ngl.contour(wks,Remap_alb,res1) # Remapped RACMO
+    diff_plot = Ngl.contour(wks, diff, res3)
+
+    # Creat multiple figures and draw, which now contains the elevation and radiation
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM, Albedo", "RACMO, Albedo", "CESM-RACMO, Albedo"]
+    # map_title = ["CESM, Albedo","RACMO, Albedo"]
+
+    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], albedo_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+
+    Ngl.overlay(plot[1], alb_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+
+    Ngl.overlay(plot[2], diff_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+
+    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_albedo',
+                        ' '.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/cesm_racmo23_latf.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_laft 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_latf_jja = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.latf.1980-1999.JJA.nc')
+    f_racmo_latf_jja_remapped = os.path.join(racmo_path,
+                                             'remapped_racmo/racmo23_GRN_monthly.latf.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 5km CESM
+    # adding the [:] transforms 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 variables
+    ncid1 = Dataset(f_cesm_lnd_climo_jja)
+    qsoil = ncid1.variables['QSOIL'][0, :, :]
+    lat1 = ncid1.variables['lat'][:]
+    lon1 = ncid1.variables['lon'][:]
+    
+    qsoil = -1 * 2.835 * 1e6 * qsoil
+    
+    # use gris as a mask to mask fsa array
+    gris_mask = ncid1.variables['gris_mask'][0, :, :]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+    
+    qsoil_mask = ma.masked_array(qsoil, mask=gris_mask.mask)
+    qsoil = qsoil_mask
+    
+    # read f_racmo_latf_jja, RACMO variables
+    ncid2 = Dataset(f_racmo_latf_jja)
+    latf = ncid2.variables['latf'][0, :, :]
+    
+    # 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 fsa array.
+    gris_mask = ncid4.variables['GrIS_mask'][:]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+    
+    latf_mask = ma.masked_array(latf, mask=gris_mask.mask)
+    latf = latf_mask
+    
+    # read f_racmo_latf_jja_remapped, the remapped RACMO file to calculate difference
+    ncid3 = Dataset(f_racmo_latf_jja_remapped)
+    remap_latf = ncid3.variables['latf'][0, :, :]
+    
+    diff = qsoil - remap_latf
+    
+    # print(np.max(diff))
+    # print(np.min(diff))
+    
+    # ------- 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_latf_JJA")
+    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 latf 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 = np.arange(-30, 10, 5)
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+    # res1.pmLabelBarOrthogonalPosF = -0.05
+    
+    res1.sfXArray = lon1
+    res1.sfYArray = lat1
+    
+    # --- for the latf 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 = np.arange(-30, 10, 5)
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+    
+    res2.sfXArray = lon4
+    res2.sfYArray = lat4
+    
+    # --- for the latf 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 = np.arange(-30, 20, 5)
+    res3.lbLabelBarOn = True  # Turn on labelbar.
+    res3.lbOrientation = "Vertical"  # Verticle labelbar
+    res3.lbLabelFontHeightF = 0.03  # 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 of 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)
+    
+    qsoil_plot = Ngl.contour(wks, qsoil, res1)
+    latf_plot = Ngl.contour(wks, latf, res2)
+    diff_plot = Ngl.contour(wks, diff, res3)
+    
+    # Creat multiple figures and draw, which now contains the elevation and radiation
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM, LHF", "RACMO, LHF", "CESM-RACMO, LHF"]
+    
+    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], qsoil_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+    
+    Ngl.overlay(plot[1], latf_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+    
+    Ngl.overlay(plot[2], diff_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+    
+    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_laft',
+                        ' '.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/cesm_racmo23_lwsd.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_lwsd 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_lwsd_jja = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.lwsd.1980-1999.JJA.nc')
+    f_racmo_lwsd_jja_remapped = os.path.join(racmo_path,
+                                             'remapped_racmo/racmo23_GRN_monthly.lwsd.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, variable of CESM
+    ncid1 = Dataset(f_cesm_lnd_climo_jja)
+    flds = ncid1.variables['FLDS'][0, :, :]
+    lat1 = ncid1.variables['lat'][:]
+    lon1 = ncid1.variables['lon'][:]
+
+    # use gris as a mask to mask flds array
+    gris_mask = ncid1.variables['gris_mask'][0, :, :]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    flds_mask = ma.masked_array(flds, mask=gris_mask.mask)
+    flds = flds_mask
+
+    # read f_racmo_lwsd_jja, variable of original RACMO
+    ncid2 = Dataset(f_racmo_lwsd_jja)
+    lwsd = ncid2.variables['lwsd'][0, :, :]
+
+    # 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 flds array.
+    gris_mask = ncid4.variables['GrIS_mask'][:]
+    gris_mask = ma.masked_equal(gris_mask, 0)
+
+    lwsd_mask = ma.masked_array(lwsd, mask=gris_mask.mask)
+    lwsd = lwsd_mask
+
+    # read f_racmo_lwsd_jja_remapped, the remapped RACMO file to calculate difference
+    ncid3 = Dataset(f_racmo_lwsd_jja_remapped)
+    remap_lwsd = ncid3.variables['lwsd'][0, :, :]
+
+    diff = flds - remap_lwsd
+
+    # print(np.max(diff))
+    # print(np.min(diff))
+
+    # ------- 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_lwsd_JJA")
+    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 lwsd contour of CESM -------
+    res1 = Ngl.Resources()
+    res1.nglDraw = False  # Don't draw individual plots
+    res1.nglFrame = False  # Don't advance frame.
+    res1.cnLineLabelsOn = False