Commit e64eb4ba authored by Evans, Katherine's avatar Evans, Katherine
Browse files

Merge remote-tracking branch 'origin' into develop

parents 78cd3784 51dcf26a
## Postprocessing of files of coupled or ice sheet only runs for use in LIVVkit
This directory holds subdirectories with scripts that will postprocess model data for use in creating the plots and data used by LIVVkit. There are also scripts to processed observation data useful for any model analysis.
These script environments take raw output from either a stand alone ice sheet model (e.g. cism-albany)
or a coupled climate model (e.g. cesm1p2, e3sm (in beta testing))
To enable or update the scripts to work on a given machine, alter the main scripts in the target model
subdirectory to load the necessary software.
To create scripts for a model not yet included here, create a new subidirectory with README files and all
the necessary files to complete the postprocessing
#### Some notes about nco.
The ncra and ncwa commands are mostly used for averaging, but in some cases they compute a sum.
Here are some alternative commands with ncra and ncwa using the -y op_typ flag that can be done within LIVVkit:
* avg returns time average
* min time min
* max time max
* ttl temporal sum
* sdn temporal std deviation
(see 3.35 in nco documentation for complete descrtion of op_types
These features could also be used with the nces (ens evg feature)
## Postprocessing of files of coupled CESM version CESM version 1.2 runs for use in LIVVkit
(note: a test version for E3SM data 'main_e3sm.sh' is also here, but its not production yet)
This directory holds subdirectories with scripts that will postprocess model data for use in creating the plots and data used by LIVVkit.
### What is needed to run postproc scripts:
*Create a working directory (with write permissions) to put all target files
*Refer to check list below to ensure all the files are accessable
*Run on OLCFs rhea system, to ensure enough memory for postprocessing CISM-Albany data
### Files needed to run coupled model postprocessing (right now a copy exists in the $WORLDWORK/cli115//4ue space)
1. b.e10.BG20TRCN.f09_g16.002
* /lnd/hist
* /atm/hist
* /glc/hist
2. surfdata_0.9x1.25_simyr1850_c110725.nc
3. /obs_data
* /CLOUDSAT_* for all months, DJF, JJA, and ANN
* /ISCCP_* for all months, DJF, JJA, and ANN
4. /racmo23_GRN_monthly (all this data is already separated by variable )
5. Greenland_5km_v1.1_SacksRev_c110629.nc
### Main script for analysis
#### Usage:
*Copy over main.sh to separate file (e.g. main_deal.sh)
*Set parameters at the top of main_deal.sh to match your environment
*Below that, opt out, if desired, to postprocessing, remapping, and plots/tables
```sh
bash main_deal.sh
```
#### Postprocessing
* loads computing environment for rhea, with hooks to add other machines
* makes necessary mask files if they do not already exist
* if 'postprocess_obs' is selected
* calls atm_obs.sh
* remaps ISCCP and CLOUDSAT climo data to model grid for selected variables for plots
* makes ice sheet area averaged climos from original climo data
* calls lnd_obs.sh
* makes climos from original racmo climo data
* makes yearly time series from original racmo monthly files
* makes ice sheet area average (but not ice sheet %) of yearly time series of racmo data
* calls atm.sh
* makes climos and ice sheet area averages
* calls monthly_atm_ts.sh
* make monthly and yearly averaged time series
* make ice sheet area average of time series
* calls lnd.sh
* makes climos and ice area average of climos
* makes daily avergages over time record
* submits rhea_tseries_lnd.pbs
* makes time series of daily avg data
* makes ice sheet average of daily time series
* makes annualized time series of daily avg data
* makes ice sheet average of annualized daily time series
* calls tseries_lnd_h0.sh
* makes time series from monthly averaged data
* makes yearly averages of monthly averaged data
* makes ice sheet average of monthly averaged data
* calls glc.sh
* select flag to compute velnorm (turn off if velocities do not exist)
* makes annual climos
* makes ice sheet average of climos
* makes time series of yearly averaged data
* area average time series (it is not masked right now)
#### Make plots and tables
* contour plots of 2m T MODEL, RACMO and difference for DJF and JJA
* plots of ANN low, med, and high clouds for MODEL, ISCCP, and CLOUDSAT
* plots of the annual cycle of low, med, and high clouds for MODEL, ISCCP, and CLOUDSAT
* plot of time series of 2m T MODEL and RACMO
* print ice sheet are and time averages of suite of variables for model and RACMO
* time series of downscaled CISM SMB
#!/bin/bash
# atmosphere
if [ "$postprocess" = "true" ]; then
# regridding step to regrid to lat lon, currently regrids all years even if not requested
if [ "$version" == "e3sm" ]; then
##ncwa -a -time $outpath/masks/mask_ne30np4.nc tmp.nc
##ncks -x -v time tmp.nc $outpath/masks/mask_ne30np4_notime.nc
#for fl in $path/atm/hist/*h0*; do
#echo "first add mask in original grid to monthly atm output in $fl"
##ncks -A -C -v gris_mask $outpath/masks/mask_ne30np4_notime.nc $fl
#ncks -A -C -v gris_mask $mask_file $fl
#done
h0atm_path=$outpath/atm/hist_regrid
#mkdir -p "$h0atm_path"
#ls $path/atm/hist/*h0* | ncremap -m $map_file -O $h0atm_path
else
h0atm_path=$path/atm/hist
fi
#run a check, is the requested data directory there?
if test -d $h0atm_path; then
echo "Directory $h0atm_path exists"
else
echo "Target directory $h0atm_path does not exist, you prob need to activate remapping"
exit
fi
echo "make atm climos"
mkdir -p $outpath/atm/climos
ncclimo -a sdd -c ${case} -s ${year_start} -e ${year_end} -m ${atm} -i ${h0atm_path} -o ${outpath}/atm/climos
if test -a $outpath/atm/climos/${case}_ANN_climo.nc; then
echo "File $outpath/atm/climos/${case}_ANN_climo.nc exists"
else
echo "File $outpath/atm/climos/${case}_ANN_climo.nc does not exist"
exit
fi
if [ "$version" == "cesm1" ]; then
echo "add mask and area variables to monthly atm climos"
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_ANN_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_DJF_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_JJA_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_01_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_02_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_03_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_04_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_05_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_06_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_07_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_08_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_09_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_10_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_11_climo.nc
ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/atm/climos/${case}_12_climo.nc
fi
echo "average over ice sheet region of atm climos "
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_ANN_climo.nc $outpath/atm/climos/${case}_ANN_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_JJA_climo.nc $outpath/atm/climos/${case}_JJA_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_DJF_climo.nc $outpath/atm/climos/${case}_DJF_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_01_climo.nc $outpath/atm/climos/${case}_01_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_02_climo.nc $outpath/atm/climos/${case}_02_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_03_climo.nc $outpath/atm/climos/${case}_03_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_04_climo.nc $outpath/atm/climos/${case}_04_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_05_climo.nc $outpath/atm/climos/${case}_05_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_06_climo.nc $outpath/atm/climos/${case}_06_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_07_climo.nc $outpath/atm/climos/${case}_07_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_08_climo.nc $outpath/atm/climos/${case}_08_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_09_climo.nc $outpath/atm/climos/${case}_09_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_10_climo.nc $outpath/atm/climos/${case}_10_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_11_climo.nc $outpath/atm/climos/${case}_11_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/atm/climos/${case}_12_climo.nc $outpath/atm/climos/${case}_12_aavg_climo.nc
echo "make yearly time series for atm h0 data"
mkdir -p $outpath/atm/tseries
for ivar in "${atm_vars[@]}"
do
. tseries_atm_h0.sh
done
echo "ATM POSTPROCESSING COMPLETE "
echo " "
fi # to postprocess or not
##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 ------------------------
#CESM1
#input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc'
#input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc'
#input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc'
#CESM2
input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/climos/B_1PAL_JJA_climo.nc'
#input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/atm/climos/B_1PAL_ANN_climo.nc'
input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/B_1PAL/postproc/lnd/climos/B_1PAL_ANN_climo.nc'
#for CESM1 and CESM2
#input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/surfdata_0.9x1.25_simyr1850_c110725.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,:,:]
#CESM1
GrisMask = ncid1.variables['gris_mask'][0,:,:]
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
##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 ------------------------
#CESM1
#input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc'
#input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc'
#input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc'
#CESM2
input_file1 = '/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_file11 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/atm/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_ANN_climo.nc'
input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/lnd/climos/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison_ANN_climo.nc'
#for CESM1 and CESM2
#input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/surfdata_0.9x1.25_simyr1850_c110725.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
#!/bin/bash
# ice sheet (CESM1 and CESM2 have the same naming conventions)
if [ "$postprocess" = "true" ]; then
calcvn=1 # create velnorm if uvel and vvel exist
#run a check, is the requested data directory there?
if test -d $path/glc/hist; then
echo "Directory $path/glc/hist exists"
else
echo "Target directory $path/glc/hist does not exist"
exit
fi
if [ $calcvn = 1 ]; then
mkdir -p $outpath/glc/hist
for (( year=${year_start}; year<=${year_end}; year++ ))
do
#echo "make velnorm from just uvel and vvel (memory constraints require separating uvel and vvel first)"
#ncrcat -O -v time,level,uvel,vvel $path/glc/hist/${case}.$glc.h.${year}-01-01-00000.nc $outpath/glc/hist/${case}.$glc.h.${year}-01-01-00000.vels.nc
#echo "created" $outpath/${case}.$glc.h.${year}-01-01-00000.vels.nc
echo "creating $outpath/${case}.$glc.h.velnorm.${year}-01-01-00000.nc file with velnorm"
ncap -O -s "velnorm=sqrt(uvel^2+vvel^2)" $path/glc/hist/${case}.$glc.h.${year}-01-01-00000.nc $outpath/glc/hist/${case}.$glc.h.velnorm.${year}-01-01-00000.nc
done
fi
echo "make ice sheet climos"
mkdir -p $outpath/glc/climos
echo /ccs/home/zender/bin_rhea/ncclimo -C ann -c $case -s $year_start -e $year_end -m $glc -h h -i $path/glc/hist -o $outpath/glc/climos
/ccs/home/zender/bin_rhea/ncclimo -C ann -c $case -s $year_start -e $year_end -m $glc -h h -i $path/glc/hist -o $outpath/glc/climos
if test -a $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc; then
echo "File $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc exists"
else
echo "File $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc does not exist"
exit
fi
echo "add area variable to file for averaging/summing"
ncks -A -v AREA $area_file $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc
echo "add masking variable to file for averaging/summing"
ncks -A -v landcover $base/Greenland_5km_v1.1_SacksRev_c110629.nc $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc
echo "average annual ice sheet climos over ice sheet"
ncwa -O -a x1,y1,lev,time -w AREA -B "landcover >= 1" $outpath/glc/climos/${case}_ANN_${year_start}01_${year_end}12_climo.nc $outpath/glc/climos/${case}_ANN_aavg.nc
#
if test -a $outpath/glc/climos/${case}_ANN_aavg.nc; then
echo "Global annual averages of $glc variables made"
else
echo "Global annual averages of $glc variables failed"
exit
fi
echo "make time_series of annual averages"
mkdir -p $outpath/glc/tseries
for i in "${glc_vars[@]}"
do
echo "make time series of yearly averages"
ncrcat -O -v $i $path/glc/hist/$case.$glc.h.*.nc $outpath/glc/tseries/$case.$glc.h.$i.nc
#ncrcat -O -v time,$i $path/glc/hist/$case.$glc.h.*.nc $outpath/glc/tseries/$case.$glc.h.$i.nc
#ncea -F -d time,first,last in.out out.nc
echo "add area variable to file for averaging/summing"
ncks -v AREA -A $area_file $outpath/glc/tseries/${case}.$glc.h.$i.nc
# this line screws up time dimension
#ncks -v landcover -A /lustre/atlas1/cli106/proj-shared/4ue/CESM-CISM/Greenland_5km_v1.1_SacksRev_c110629.nc $outpath/glc/tseries/${case}.$glc.h.$i.nc
echo "average variable over whole ice sheet domain"
ncwa -O -a x1,y1 -w AREA $outpath/glc/tseries/${case}.$glc.h.$i.nc $outpath/glc/tseries/${case}.$glc.h.${i}_aavg.nc
#ncwa -O -a x1,y1 -w AREA -B "landcover >= 1" $outpath/glc/tseries/${case}.$glc.h.$i.nc $outpath/glc/tseries/${case}.$glc.h.${i}_aavg.nc
if test -a $outpath/glc/tseries/$case.$glc.h.${i}_aavg.nc; then
echo "Timeseries of global annual averages of $glc for $i made"
else
echo "Timeseries of global annual averages of $glc for $i failed"
exit
fi
done
fi # to postprocess or not
#if [ "$make_plots" = "true" ]; then
#echo "make the plots"
#mkdir -p "$outpath/glc/plots"
#pwd
#python plot_files/TimeSeries_acab.py
#fi # to make plots or not
#!/bin/bash
# land
if [ "$postprocess" = "true" ]; then
if [ "$version" == "e3sm" ]; then
#ncwa -a -time $outpath/masks/mask_ne30np4_lnd.nc $outpath/masks/mask_ne30np4_lnd_notime.nc
#for fl in $path/lnd/hist/*h0*; do
#echo "first add mask in original grid to monthly lnd output in $fl"
#ncks -A -C -v gris_mask $outpath/masks/mask_ne30np4_lnd_notime.nc $fl
#done
# regridding step to regrid to lat lon, currently regrids all years even if not requested
#h0lnd_path=$outpath/lnd/hist_regrid
#mkdir -p "$h0lnd_path"
#ls $path/lnd/hist/*h0* | ncremap -m $map_file -O $h0lnd_path
h0lnd_path=$path/archive/lnd/hist
else
h0lnd_path=$path/lnd/hist
fi
echo "make land h0 climos"
mkdir -p "$outpath/lnd/climos"
ncclimo -c $case -s ${year_start} -e ${year_end} -m ${lnd} -i ${h0lnd_path} -o ${outpath}/lnd/climos
if test -a $outpath/lnd/climos/${case}_ANN_climo.nc; then
echo "File $outpath/lnd/climos/${case}_ANN_climo.nc exists"
else
echo "File $outpath/lnd/climos/${case}_ANN_climo.nc does not exist"
exit
fi
if [ "$version" == "cesm1" ]; then
h0lnd_path=$outpath/lnd/climos
echo "average land monthly climos over ice sheet region"
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/lnd/climos/${case}_ANN_climo.nc $outpath/lnd/climos/${case}_ANN_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/lnd/climos/${case}_JJA_climo.nc $outpath/lnd/climos/${case}_JJA_aavg_climo.nc
ncwa -O -a lat,lon -w area -B "gris_mask == 1" $outpath/lnd/climos/${case}_DJF_climo.nc $outpath/lnd/climos/${case}_DJF_aavg_climo.nc
mkdir -p "$outpath/lnd/davg"
#echo "make annualized daily averages for land for all variables"
#ls $path/lnd/hist/$case.$lnd.h1.*.nc | /ccs/home/zender/bin_rhea/ncclimo --job_nbr=12 -c $case -C dly -s ${year_start} -e ${year_end} -o $outpath/lnd/davg -m $lnd
#/ccs/home/zender/bin_rhea/ncclimo --job_nbr=12 -c $case -C dly -s ${year_start} -e ${year_end} -o $outpath/lnd/davg -m $lnd -i $path/lnd/hist
#echo "loop over each variable to create and process time series for h1 vars"
#qsub rhea_tseries_lnd.sh
#for ivar in "${lnd_h1_vars[@]}"
#do
#echo $ivar
#echo "concatenate all daily dates to be averaged into annualized land time series for $ivar"
#ncrcat -O -F -v $ivar $outpath/lnd/davg/*$year_start*$year_end*_climo.nc $outpath/lnd/davg/$case.$lnd.h1.$ivar.davg.${year_start}-${year_end}.nc
#
#echo "add mask and area variables to annualized daily land averages for $ivar"
#ncks -A -v gris_mask,area $outpath/masks/$case.$lnd.gris_mask.nc $outpath/lnd/davg/$case.$lnd.h1.$ivar.davg.${year_start}-${year_end}.nc