updates to process E3SM v1 data

postproc/e3sm1p0/e3sm_rad_aavg.py

+##On my laptop I installed PyNGL and PyNIO to use Nio and Ngl module for reading *.nc file and for contour plot
+##On titan, I did not install PyNGL and PyNIO, so import Dataset from netCDF4 to read *.nc file
+##probably, I will the following command for contourplot in titan python "from mpl_toolkits.basemap import Basemap"
+
+
+#!/usr/bin/env python3
+
+import numpy as np
+import numpy.ma as ma
+import os
+from netCDF4 import Dataset
+
+
+# ---------------- Data sources ------------------------
+input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/climo_regrid/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
+input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/climo_regrid/20180215.DECKv1b_H1.ne30_oEC.edison_ANN_climo.nc'
+
+# --------------------------------------------------------------
+# --------------------------------------------------------------
+
+# input_file1 get FLDS (mean downwelling longwave rad),FIRA (mean net IR longwave rad, absorbed)
+# FLDS(0,:,:), FIRA(0,:,:)
+# Matrix gris_area(Greenland ice area) has values 0 and none
+
+# Goal: print out area-weighted average of SW_d(268), SW_net(61), LW_d(235), LW_net(-46), R_net(15), SHF(7), LHF(-8),
+# these variables have size[192,288]
+# FSDS: atmospheric incident solar radiation
+# FSA: absorbed solar radiation
+# FLDS: atmospheric longwave radiation
+# FIRA: net infrared (longwave) radiation
+# FIRE: emitted infrared (longwave) radiation
+# FSH: sensible heat
+# QSOIL: ground evaporation over ice landunits
+
+
+# ---read input_file2, probably get area with ice sheet with nonzero values
+#ncid2 = Dataset(input_file2, mode='r')
+#ices = ncid2.variables['PCT_GLC_ICESHEET'][:] 
+#ncid2.close()
+
+# ---read input_files
+ncid1 = Dataset(input_file1, mode='r')
+#ncid11 = Dataset(input_file11, mode='r')
+ncid3 = Dataset(input_file3, mode='r')
+
+FSDS = ncid1.variables['FSDS'][0,:,:] 
+FSA = ncid1.variables['FSA'][0,:,:] 
+FLDS = ncid1.variables['FLDS'][0,:,:] 
+FIRA = ncid1.variables['FIRA'][0,:,:]
+FSH = ncid1.variables['FSH'][0,:,:] 
+QSOIL = ncid1.variables['QSOIL'][0,:,:] 
+RAIN = ncid3.variables['RAIN'][0,:,:] 
+SNOW = ncid3.variables['SNOW'][0,:,:] 
+QRUNOFF = ncid3.variables['QRUNOFF'][0,:,:] 
+#T2m = ncid11.variables['TREFHT'][0,:,:] 
+T2m = ncid1.variables['TSA'][0,:,:] 
+
+#E3SM
+GrisMask = ncid1.variables['gris_mask'][:,:] 
+
+SW_d  = FSDS 
+SW_net = FSA
+LW_d = FLDS
+LW_net = -1.0 * FIRA
+R_net = FSA-FIRA
+SHF = -1.0 * FSH
+LHF = -1.0 * 2.835e6 * QSOIL
+#SMB = QICE * 3.1567e7
+SMB = RAIN + SNOW - QRUNOFF - QSOIL * 3.1567e7
+T2mC = T2m - 273.15
+
+# ---grid cell area on Greenland. For other regions, the value is 9.999e35.
+area = ncid1.variables['area'][:] 
+
+ncid1.close()
+#ncid11.close()
+ncid3.close()
+
+# # # --- Option I: use greenland mask to mask the area where the ice sheet fraction is zero
+# # # Note: some non-greenland region has non-zero ice sheet fraction.
+Gris_mask = ma.masked_equal(GrisMask,0)
+area_maskice = ma.masked_array(area, mask=Gris_mask.mask)
+
+
+# # # --- Option II: use icesheet fraction mask to mask the area where the ice sheet fraction is zero
+# ices_mask = ma.masked_equal(ices,0)
+# area_maskice = ma.masked_array(area, mask=ices_mask.mask)
+
+
+# # # --- Option III: first multiply area with ices, then use greenland mask to mask the area where the ice sheet fraction is zero
+## -- this option gives the closest results to the Table 1 in the paper.
+#areaices = np.multiply(area,ices)
+#Gris_mask = ma.masked_equal(GrisMask,0)
+#area_maskice = ma.masked_array(areaices, mask=Gris_mask.mask)
+
+
+swdout=(ma.average(SW_d,weights=area_maskice))
+swnetout=(ma.average(SW_net,weights=area_maskice))
+lwdout=(ma.average(LW_d,weights=area_maskice))
+lwnetout=(ma.average(LW_net,weights=area_maskice))
+rnetout=(ma.average(R_net,weights=area_maskice))
+shfout=(ma.average(SHF,weights=area_maskice))
+lhfout=(ma.average(LHF,weights=area_maskice))
+smbout=(ma.average(SMB,weights=area_maskice))
+t2cout=(ma.average(T2mC,weights=area_maskice))
+
+print "Sw_d model:", swdout
+print "Sw_net model:", swnetout
+print "Lw_d model:", lwdout
+print "Lw_net model:", lwnetout
+print "Rnet model:", rnetout
+print "SHF model:", shfout
+print "LHF model:", lhfout
+print "SMB model:", smbout
+#print "T2M C atm model:", t2cout
+print "T2M C land model:", t2cout
+
+

postproc/e3sm1p0/main_e3sm.sh

@@ -24,17 +24,10 @@ area_file=/lustre/atlas1/$project/world-shared/$USER/surfdata_0.9x1.25_simyr1850
 map_file=/ccs/home/zender/data/maps/map_ne30np4_to_fv129x256_aave.20150901.nc
 mask_file=/lustre/atlas1/$project/world-shared/4ue/$case/postproc/mask_ne30np4_notime.nc
 version=e3sm
-if [[ $version = cesm1 ]]; then
-atm_vars=("PSL" "CLDLOW" "CLDHGH" "CLDMED" "CLDTOT" "TREFHT" "TS" "TSMX" "TSMN" "PRECT" "PRECC" "PRECL")
-lnd_h0_vars=("QICE" "TREFMNAV" "TREFMXAV" "TSA" "RAIN" "SNOW")
-lnd_h1_vars=("QICE_FRZ" "QICE_MELT" "QSOIL_ICE" "QSNOMELT_ICE" "RAIN" "SNOW")
-elif [[ $version = e3sm ]]; then
 atm_vars=("PSL" "CLDLOW" "CLDHGH" "CLDMED" "CLDTOT" "TREFHT" "TS" "TSMX" "TSMN" "PRECC" "PRECL")
 lnd_h0_vars=("FSDS" "QRUNOFF" "QSOIL" "TSA" "RAIN" "SNOW" "FIRA" "FLDS" "FSH" "QSNOMELT")
-#glc_vars=("acab" "thk")
-fi
 
-postprocess=true
+postprocess=false
 postprocess_obs=false
 make_plots=true
 
@@ -43,18 +36,9 @@ mkdir -p $outpath
 echo "$version"
 
 #postprocess all components of coupled model 
-if [ "$version" == "cesm1" ]; then
-	glc=cism
-	lnd=clm2 
-	atm=cam2 
-elif [ "$version" == "e3sm" ]; then
 	glc=cism
 	lnd=clm2 
 	atm=cam 
-else
-  	echo "incorrect model version detected"
-	exit	
-fi
 
 #run a check, is the requested data directory there?
 if test -d $path; then
@@ -96,21 +80,7 @@ else
 	exit
 fi
 
-# this is a hack- all files from E3SM runs we care about should have a gris_mask in the data
-if [ "$version" == "cesm1" ]; then
-
-if test -a $outpath/masks/${case}.$lnd.gris_mask.nc; then
-        echo "already created gris_mask and area file, moving on"
-else
-        echo "gris_mask file does not exist, creating"
-
-	mkdir -p $outpath/masks
-       	ncks -O -v gris_mask,area $path/lnd/hist/${case}.${lnd}.h0.${year_start}-01.nc $outpath/masks/${case}.${lnd}.gris_mask_time.nc 
-       	ncwa -a time $outpath/masks/${case}.${lnd}.gris_mask_time.nc $outpath/masks/${case}.${lnd}.gris_mask.nc 
-
-fi # if mask file needs to be created
-
-fi  # CESM vs E3SM
+# this is a hack for cesm runs, all files from E3SM runs we care about should have a gris_mask in the data
 
 if [ "$postprocess_obs" = "true" ]; then
 
@@ -119,12 +89,7 @@ echo "process and remap cloud obs"
 echo "process RACMO"
 . lnd_obs.sh 
 
-if [ "$version" = "cesm1" ]; then
-echo "remap RACMO to CESM"
-./remap_RACMO_to_CESM_inline.sh
-elif [ "$version" = "e3sm" ]; then
 ./remap_RACMO_to_CESM2_inline.sh
-fi #CESM version
 
 fi
 
@@ -134,33 +99,11 @@ echo "executing post processing of atm data for validation"
 echo "executing post processing of lnd data for validation"
 . lnd.sh  
 
-if [ "$version" = "cesm1" ]; then
-echo "executing post processing of glc data for validation"
-. glc.sh
-
-fi
 
 if [ "$make_plots" = "true" ]; then
 echo "make the atmosphere plots"
 mkdir -p "$outpath/atm/plots"
 
-if [ "$version" = "cesm1" ]; then
-
-python plot_files/CESM_cldlow_ANN.py 
-python plot_files/CESM_cldhgh_ANN.py
-python plot_files/CESM_cldtot_ANN.py 
-python plot_files/yearly_cycle_CLDTOT.py 
-python plot_files/yearly_cycle_CLDHGH.py 
-python plot_files/yearly_cycle_CLDLOW.py 
-
-python plot_files/TimeSeries_TREFHT.py
-
-echo "Areal annual avgd values of radiation vars: CESM"
-python model_rad_aavg.py
-
-elif [ "$version" = "e3sm" ]; then
-
-python plot_files/E3SM_cldlow_ANN.py 
 python plot_files/E3SM_cldhgh_ANN.py
 python plot_files/E3SM_cldtot_ANN.py 
 python plot_files/yearly_cycle_CLDTOT_e3sm.py 
@@ -170,8 +113,6 @@ python plot_files/yearly_cycle_CLDLOW_e3sm.py
 echo "Areal annual avgd values of radiation vars: E3SM"
 python e3sm_rad_aavg.py
 
-fi
-
 echo "Areal annual avgd values of radiation vars: RACMO"
 python racmo_rad_aavg.py
 
@@ -203,12 +144,5 @@ python plot_files/E3SM_RACMO23_smb.py
 python plot_files/TimeSeries_smblike.py 
 fi
 
-if [ "$version" = "cesm1" ]; then
-echo "make the glc plots"
-mkdir -p "$outpath/glc/plots"
-
-python plot_files/TimeSeries_acab.py
-fi # cesm version
-
 fi # to make plots or not
 

postproc/e3sm1p0/plot_files/E3SM_cldhgh_ANN.py

@@ -6,8 +6,8 @@ import Nio, Ngl, os
 
 # ---------------- Data source in TITAN ------------------------
 input_file0 = '/lustre/atlas1/cli115/world-shared/4ue/Greenland_5km_v1.1_SacksRev_c110629.nc'
-input_file8 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/lnd/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_JJA_climo.nc'
-input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/atm/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_JJA_climo.nc'
+input_file8 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
+input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/atm/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
 input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/ISCCP_ANN_climo.nc'
 input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/CLOUDSAT_ANN_climo.nc'
 input_file22 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/remapped_racmo/racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.JJA.nc'

postproc/e3sm1p0/plot_files/E3SM_cldlow_ANN.py

@@ -6,8 +6,8 @@ import Nio, Ngl, os
 
 # ---------------- Data source in TITAN ------------------------
 input_file0 = '/lustre/atlas1/cli115/world-shared/4ue/Greenland_5km_v1.1_SacksRev_c110629.nc'
-input_file8 = '/lustre/atlas1/cli115/world-shared/4ue//postproc/lnd/climos/_JJA_climo.nc'
-input_file1 = '/lustre/atlas1/cli115/world-shared/4ue//postproc/atm/climos/_JJA_climo.nc'
+input_file8 = '/lustre/atlas1/cli115/world-shared/4ue//postproc/lnd/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
+input_file1 = '/lustre/atlas1/cli115/world-shared/4ue//postproc/atm/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
 input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/ISCCP_ANN_climo.nc'
 input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/CLOUDSAT_ANN_climo.nc'
 input_file22 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/remapped_racmo/racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.JJA.nc'

postproc/e3sm1p0/plot_files/E3SM_cldtot_ANN.py

@@ -6,8 +6,8 @@ import Nio, Ngl, os
 
 # ---------------- Data source in TITAN ------------------------
 input_file0 = '/lustre/atlas1/cli115/world-shared/4ue/Greenland_5km_v1.1_SacksRev_c110629.nc'
-input_file8 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/lnd/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_JJA_climo.nc'
-input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/atm/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_JJA_climo.nc'
+input_file8 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
+input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/atm/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc'
 input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/ISCCP_ANN_climo.nc'
 input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/CLOUDSAT_ANN_climo.nc'
 input_file22 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/remapped_racmo/racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.JJA.nc'

postproc/e3sm1p0/plot_files/TimeSeries_smblike.py

@@ -9,10 +9,10 @@ plt.switch_backend('agg')
 from netCDF4 import Dataset
 
 # ---------------- Data source in rhea, the two files save the same data ------------------------
-input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/tseries/B_1PAL.clm2.aavg.h0.yrly.RAIN.nc'
-input_file12 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/tseries/B_1PAL.clm2.aavg.h0.yrly.SNOW.nc'
-input_file13 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/tseries/B_1PAL.clm2.aavg.h0.yrly.QRUNOFF.nc'
-input_file14 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/tseries/B_1PAL.clm2.aavg.h0.yrly.QSOIL.nc'
+input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.RAIN.nc'
+input_file12 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.SNOW.nc'
+input_file13 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QRUNOFF.nc'
+input_file14 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QSOIL.nc'
 input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc'
 # --------------------------------------------------------------
 
@@ -93,7 +93,7 @@ ax.set_xlabel('(years)',fontsize=14)
 ##ax.text(-12,5,'T-2m, JJA',fontsize=18)
 ##ax.legend([s1,s2],['>20%','>99%'],loc='lower right')
 
-plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/plots/TimeSeries_smblike.ps')
+plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/TimeSeries_smblike.ps')
 #plt.savefig('TimeSeries_SMB.ps')
 plt.close()
 
@@ -124,7 +124,7 @@ ax.set_ylabel('RACMO SMB(km/m2/yr)',fontsize=14)
 ax.set_ylim([0,500])
 ax.set_xticklabels('')
 
-plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/plots/BoxPlot_smblike.ps')
+plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/BoxPlot_smblike.ps')
 #plt.savefig('BoxPlot_SMB.ps')
 plt.close()