preprocessor.ESGF

#!/bin/bash
# Sample preprocessor in shell
#--------------------------------------------------------------------------
#+ Settings
#--------------------------------------------------------------------------
set -eo pipefail     # abort execution on error (e), abort execution on undefined variable (u), abort execution on error in pipe command (o pipefail)
BCTOOL_DEBUG=${BCTOOL_DEBUG:-0}
test "${BCTOOL_DEBUG}" -ne 0 && set -x

# Read arguments
sformateddate=$1     # Initial date to process as YYYY-MM-DD_HH:MM:SS (Ex. 1983-08-27_00:00:00)
eformateddate=$2     # end date to process
BCTABLE=$3           # BCtable file

# Set a full path (if relative)
test "${BCTABLE}" = ${BCTABLE#/} && BCTABLE=`pwd`/${BCTABLE}

# Set default filters to be applied, before and after filter given in the BCTABLE
wrfpostfilter="convert2grb|set_grb_code|set_grb_ltype|set_extension(grb)|unleap"
prefilter="time_slice"
postfilter="maskregion|${wrfpostfilter}"

#--------------------------------------------------------------------------
#Functions
#--------------------------------------------------------------------------
function timemsg(){
  echo "[$(date '+%Y-%m-%d %H:%M:%S')] $*"
}

function errormsg(){
  echo "Error: $*"
  exit 1
}

function debugmsg(){
  echo "<debugmsg $(date '+%Y-%m-%d %H:%M:%S')> $*" >&2
}

function bctable_header(){
  val=$1
  awk '/\<'${val}'\>/ {print $2}' ${BCTABLE}
}

function bctable_variables(){
  awk '{if (a) print; if (substr($1,1,3)=="---") a=1}' ${BCTABLE}
}

function findfile(){
  local filespec=$1
  local year=$2
  find ${drsdir} -type f -name "${filespec}" \
    | sort \
    | while read f; do
      idate=${f//*_/}; idate=${idate//-*/}
      lendate=${#idate}
      iyear="${year}0101000000"; iyear="${iyear:0:${lendate}}"
      fdate=${f//*_/}; fdate=${fdate//*-/}; fdate=${fdate//.*/}
      fyear="${year}1231180000"; fyear="${fyear:0:${lendate}}"
      if test ${iyear} -ge ${idate} -a ${fyear} -le ${fdate}; then
        echo ${f}
        break
      fi
    done
}

function copy(){
  orig=$1
  dest=$2
  cp ${orig} ${dest}
  chmod ug+w ${dest}
}

function apply_filter(){
  filters=$1
  filterlist=${filters##*BEGIN}
  filterlist=${filterlist%%END*}
  for filt in ${filterlist//|/ }
  do
    timemsg "    Applying filter: ${filt}"
    read filter_name filter_arguments <<< $(echo ${filt} | tr ')(,' '   ')
    ${filter_name} ${tmpDir}/${variable}_${year} ${filter_arguments}
    if test $? -ne 0; then
      timemsg "    FAILED"
      exit
    fi
  done
}


function time2Interval(){
  ifile=$1
  ofile=$2
  nrec=$(cdo -s ntime ${ifile})
  cdo -s seltimestep,1 ${ifile} s1
  cdo -s settime,00:00:00 s1 s0
  cdo -s seltimestep,${nrec} ${ifile} s1
  cdo -s setdate,$(date --utc +%Y-%m-%d --d "`cdo -s showtimestamp s1` 1day") s1 s2
  cdo -s settime,00:00:00 s2 s1
  cdo -s mergetime s0 ${ifile} s1 s3
  cdo -s inttime,$(cdo -s showdate s0 | tr -d ' '),00:00:00,${interval} s3  ${ofile}
  rm s0 s1 s2 s3
}
#To test: time2Interval tsl_2033.nc pp.nc


function timeRange2Interval(){
  ifile=$1
  ofile=$2
  ststamp=$3
  etstamp=$4
  nrec=$(cdo -s ntime ${ifile})
  test ${nrec} -eq 0 && nrec=1
  cdo -s seltimestep,1 ${ifile} s1
  cdo -s setdate,${ststamp} s1 s2
  cdo -s settime,$(echo ${ststamp} | awk -FT '{print $2}') s2 s0
  cdo -s seltimestep,${nrec} ${ifile} s1
  cdo -s setdate,${etstamp} s1 s2
  cdo -s settime,$(echo ${etstamp} | awk -FT '{print $2}') s2 s1
  cdo -s mergetime s0 s1 s3
  cdo -s inttime,$(cdo -s showdate s0 | tr -d ' '),00:00:00,${interval} s3  ${ofile}
  rm s0 s1 s2 s3
}
#To test: timeRange2Interval tsl_2033.nc pp.nc 2033-12-24T00:00:00 2034-01-04T00:00:00

function recordOfClosestDate(){
  ifile=$1
  targetdate=$2
  tdate=$(date +%s -d ${targetdate})
  cdo -s showtimestamp ${ifile} | tr ' ' '\n' | grep -v '^$' | while read datestr; do date +%s -d $datestr; done \
    | awk '$1-'${tdate}'<0 {a+=1} END {print a}'
}

function dump_coord_bnds(){
  coord=$1
  thefile=$2
  ncdump -v ${coord} ${thefile} \
    | awk -F'[,;]' '/\}/{a=0} {if(a) printf("%.5f %.5f\n", $1,$2)} /^ '${coord}' =/{a=1}' 
}

function write_static_files(){
cat << EOF > myzaxis2m.dat
zaxistype = height
size      = 1 
levels    = 2
EOF
cat << EOF > myzaxis10m.dat
zaxistype = height
size      = 1 
levels    = 10
EOF
cat << EOF > region.dat      
-40.0 10.0
-80.0 90.0
100.0 90.0
60.0 10.0
EOF
}

#--------------------------------------------------------------------------
#+ Filter functions
#--------------------------------------------------------------------------
# PIPE-START FUNCTIONS
# "Input" file in these functions is the output file
#
function only_closest_to_sdate(){
  ofile=$1
  ifile=$(findfile ${filespec} ${year})
  cdo seltimestep,$(recordOfClosestDate ${ifile} ${sformateddate/_/T}) ${ifile} s59
  cdo settime,${sformateddate//*_/} -setdate,${sformateddate//_*/} s59 ${ofile}
  rm s59
}

function time_slice(){
  ofile=$1
  if test "${frequency}" == "fx"; then
    ifile=$(find ${drsdir} -type f -name "${filespec}" | head -n 1)
    copy ${ifile} ${ofile}
  else
    ifile=$(findfile ${filespec} ${year})
    cdo -s seldate,${sformateddate/_/T},${eformateddate/_/T} ${ifile} s84
    cdo -s selyear,${year} s84 ${ofile} ; rm s84
    calendar=$(cdo -s -sinfo ${ofile} | grep Calendar | awk '{print $10}')
  fi
}

function use_era_interim(){
  ofile=$1
  ERAINPATH=$(bctable_header era_interim_path)
  if test -z "${ERAINPATH}"; then
    timemsg "Error: 'era_interim_path' not set in ${BCTABLE}, that calls use_era_interim filter"
    exit 2
  fi
  arbitrary_year=2014 # It doesn't really matter as long as it is available in ERAINPATH
  ifile="${ERAINPATH}/${arbitrary_year}/${imm}/INTERIM075_${arbitrary_year}${imm}_SFC_${gribcode}.128_an.grb"
  calendar=$(cdo -s -sinfo ${ifile} | grep Calendar | awk '{print $10}')
  copy ${ifile} ${ofile}
}


# PIPE FUNCTIONS
# Input file is overwritten by the output
#
function maskregion(){
  # region.dat contains the outer mask polygon
  ifile=$1
  cdo -s maskregion,region.dat ${ifile} s47
  mv s47 ${ifile}
}

function ringregion(){
  # region.dat contains the outer mask polygon
  # Here, the interior is also set to missings (except for the first record)
  # /!\ Not usable for climate simulations (requires sst_update = .F. and non-lat/lon projection).
  # /!\ Fail at real.exe step (see details in dyn_em/module_initialize_real.F)
  ifile=$1
  cdo seltimestep,1 ${ifile} s41
  nrec=$(cdo -s ntime ${ifile})
  cdo seltimestep,$(seq 2 ${nrec}| tr '\n' ',') ${ifile} s42
  cdo -s maskregion,region.dat s41 s43
  cdo -s -setctomiss,-9999 -setclonlatbox,-9999,-5,15,40,50 \
      -maskregion,region.dat s42 s44
  rm ${ifile} # To avoid appending
  cdo cat s43 s44 ${ifile}
  rm -f s4[1234]
}

function is_land_mask(){
  ifile=$1
  land_mask="${ifile}"
}

function only_ic(){
  ifile=$1 # unused
  is_only_ic=1
}

function percent2one(){
  ifile=$1
  cdo -s mulc,0.01 ${ifile} s73
  mv s73 ${ifile}
}

function celsius2K(){
  ifile=$1
  cdo -s addc,273.15 ${ifile} s74
  mv s74 ${ifile}
}

function set_start_time(){
  ifile=$1
  cdo -s setdate,${sformateddate//_*/} ${ifile} s58
  cdo -s settime,${sformateddate//*_/} s58 ${ifile}
  rm s58
}


function fixed2interval(){
  ifile=$1
  timeRange2Interval ${ifile} s53 ${sformateddate/_/T} ${eformateddate/_/T}
  mv s53 ${ifile}
}

function tinterp2interval(){
  ifile=$1
  time2Interval ${ifile} s54
  mv s54 ${ifile}
}

function shift_time(){
  ifile=$1
  tshift=$2
  cdo -s shifttime,${tshift} ${ifile} s83
  mv s83 ${ifile}
}

function unleap(){
  ifile=$1.grb
  if [[ $noleap == "true" ]] && [[ $frequency != "fx" ]] ; then 
  	nonstandard="360_day" # not functional for the "360_day" calendar
        check=`echo $(cdo -s -info ${ifile} | grep "02-29 " && echo true || echo false) | awk '{print $NF}'` # checking if date in the file
        if [[  $check == "true" ]] && [[  $calendar != $nonstandard ]] ; then
          timemsg "    		   removing Feb 29"
          cdo -s del29feb ${ifile} s86; mv s86 ${ifile}
        else
          timemsg "   		   Feb 29 not in the file or nonstandard calendar - not done"  
        fi
  fi
}

function set_hybrid_levels(){
  ifile=$1
  cdo -s invertlev ${ifile} s56
  dump_coord_bnds ap_bnds s56 | awk '{print $1}' > ap_bnds.txt 
  dump_coord_bnds ap_bnds s56 | awk '{print $2}' | tail -1 >> ap_bnds.txt
  dump_coord_bnds b_bnds  s56 | awk '{print $1}' > b_bnds.txt
  dump_coord_bnds b_bnds  s56 | awk '{print $2}' | tail -1 >> b_bnds.txt
  # Create hybrid coordinate coefficients for WRF
  paste -d ' ' ap_bnds.txt b_bnds.txt | nl -v0 > ecmwf_coeffs
  # Create a Z-axis description file (myzaxisinvert.dat) to convert nc files to grib ones.                            
  nlevels=$(ncdump -h s56 | grep lev | awk '{print $3}' | head -1)  # number of levels read from the netcdf file
  cat << EOF > myzaxisinvert.dat
zaxistype = hybrid 
size      = ${nlevels}
levels    = $(seq -s ' ' 1 ${nlevels})
vctsize   = $((2*${nlevels}+2))
vct       = $(cat ap_bnds.txt b_bnds.txt)
EOF
  cdo -s setzaxis,myzaxisinvert.dat s56 s57
  ncks -O -C -x -v ps s57 ${ifile}
  rm s56 s57 *_bnds.txt myzaxisinvert.dat 
}

function set_hybrid_levels_python(){
  ifile=$1
  cdo -s invertlev ${ifile} s56
  python3 - << EOF
from netCDF4 import Dataset
import numpy as np
with Dataset("${ifile}","r") as ncid:
  ap_bnds = ncid["ap_bnds"]
  ap_bnds = np.concatenate((ap_bnds[:,1],ap_bnds[-1:,0]))
  b_bnds = ncid["b_bnds"]
  b_bnds = np.concatenate((b_bnds[:,1],b_bnds[-1:,0]))
  hlev = np.vstack((ap_bnds,b_bnds))
  with open("ecmwf_coeffs") as fid:
    for i,a in enumerate(hlev.T):
      fid.write("{:d} {:10.3f} {:10.3f}".format(i,a[0],a[1]))
  nlevels = ncid['lev'].size
  with open("myzaxisinvert.dat") as fid:
    fid.write("zaxistype = hybrid")
    fid.write("size      = {}".format(nlevels))
    fid.write("levels    = {}".format(" ".join(map(str,range(1,nlevels+1))))
    fid.write("vctsize   = {}".format(2*nlevels+2))
    fid.write("vct       = {} {}".format(" ".join(map(str,ap_bnds)), " ".join(map(str,b_bnds)))
EOF
  cdo -s setzaxis,myzaxisinvert.dat s56 s57
  cdo -s delname,ps,ap_bnds,b_bnds s57 ${ifile} 
  rm s56 s57 *_bnds.txt myzaxisinvert.dat 
}

function remapnn(){
  ifile=$1
  # Fill missing value areas since this tool is usually applied before a masking operation
  export REMAP_EXTRAPOLATE=on
  cdo -s remapnn,${land_mask} ${ifile} s60
  mv s60 ${ifile} 
}

function sea_masked(){
  ifile=$1
  # Assumes .grb extension in the mask file
  cdo -s ifnotthen ${land_mask} ${ifile} s58
  mv s58 ${ifile} 
}

function rename(){
  ifile=$1
  bname=$2
  newname="$(dirname ${ifile})/${ifile//*_/${bname}_}" # preserve _year.ext
  current_bname=$(basename ${ifile})
  if test "${current_bname//_*/}" != "${bname}"; then
    # It could be that the name was already updated in a previous year processing
    # This dangerous stuff needs to be re-implemented in the main loop.
    mv ${ifile} ${newname}
    test "${ifile}" = "${land_mask}" && land_mask=${newname}
    variable=${bname} # Poor programming practice: renames loop variable within a function call
  fi
}

function split_soil_temp_grb(){
  ifile=$1
  ndepth=$(ncdump -h ${ifile} | grep depth | awk '{print $3}' | head -1)
  grib_params=( 139 170 183 236 ) 	# grib codes from ECMWF parameter database (https://apps.ecmwf.int/codes/grib/param-db); for up to 4 layers 
  depths=$(cdo -s showlevel ${ifile})
  cdo -s splitlevel ${ifile} soilt
  # Calculate boundaries of the layers
  lower_bnd=0
  for nd in `seq $ndepth`; do
    depth=`echo $depths | awk -v var=$nd '{print $var}'` 
    upper_bnd=`awk "BEGIN {x=$lower_bnd; y=$depth; z=(2*y-x); print z}"`
    lower_cm=`awk "BEGIN {x=$lower_bnd; z=(100*x); print z}"`
    upper_cm=`awk "BEGIN {x=$upper_bnd; z=(100*x); print z}"`
    delta_cm=`awk "BEGIN {x=$upper_cm; y=$lower_cm; z=1/(10*(x-y)); print z}"`
cat << EOF > myzaxis${lower_cm}_${upper_cm}.dat      
zaxistype = depth_below_land 
size      = 1
name      = depth
longname  = depth_below_land
units     = cm
lbounds    = $lower_cm
ubounds    = $upper_cm
EOF
    lower_bnd=$upper_bnd
    grib_param=${grib_params[`expr $nd-1`]}
    # Convert to grib format
    cdo -s -f grb setzaxis,myzaxis${lower_cm}_${upper_cm}.dat -setltype,${ltype} -chparam,-1,$grib_param soilt*$depth.nc soilt0${nd}_${year}.grb
  done  
  rm ${ifile} myzaxis${lower_cm}_${upper_cm}.dat soilt*.nc
}

function split_soil_mois_grb(){
  # NOTE: Depth variable inicating the layer depth  (middle) should exist in the file
  ifile=$1
  ndepth=$(ncdump -h ${ifile} | grep depth | awk '{print $3}' | head -1)
  depths=$(cdo -s showlevel ${ifile})
  cdo -s splitlevel ${ifile} soilm 
  # Calculate boundaries of the layers
  lower_bnd=0
  for nd in `seq $ndepth`; do
    depth=`echo $depths | awk -v var=$nd '{print $var}'` 
    upper_bnd=`awk "BEGIN {x=$lower_bnd; y=$depth; z=(2*y-x); print z}"`
    lower_cm=`awk "BEGIN {x=$lower_bnd; z=(100*x); print z}"`
    upper_cm=`awk "BEGIN {x=$upper_bnd; z=(100*x); print z}"`
    delta_cm=`awk "BEGIN {x=$upper_cm; y=$lower_cm; z=1/(10*(x-y)); print z}"`
cat << EOF > myzaxis${lower_cm}_${upper_cm}.dat      
zaxistype = depth_below_land 
size      = 1
name      = depth
longname  = depth_below_land
units     = cm
lbounds    = $lower_cm
ubounds    = $upper_cm
EOF
    lower_bnd=$upper_bnd
    grib_param=`expr 38 + $nd`
    # Convert to grib format
    cdo -s -f grb setzaxis,myzaxis${lower_cm}_${upper_cm}.dat -setltype,${ltype} -chparam,-1,$grib_param soilm*$depth.nc soilm0${nd}_${year}.grb
    # Convert to m3/m3
    cdo -s mulc,${delta_cm} soilm0${nd}_${year}.grb tmp.grb
    mv tmp.grb soilm0${nd}_${year}.grb
  done   
  rm ${ifile} myzaxis${lower_cm}_${upper_cm}.dat soilm*.nc
}

# set grid codes with ecCodes package developed by ECMWF
function ecc_set_grb_code(){
  ifile=$1
  grib_set -s indicatorOfParameter=${gribcode} ${ifile} s51
  mv s51 ${ifile}
}

# set grid level type with ecCodes package developed by ECMWF
function ecc_set_grb_ltype(){
  ifile=$1
  grib_set -s indicatorOfTypeOfLevel=${ltype} ${ifile} s52
  mv s52 ${ifile}
}

function convert2grb(){
  ifile=$1
  cdo -s -a -f grb copy ${ifile} s56
  mv s56 ${ifile}
}

function set_grb_code(){
  ifile=$1
  cdo -s chparam,1,${gribcode} ${ifile} s51
  mv s51 ${ifile}
}

function set_grb_ltype(){
  ifile=$1
  cdo -s setltype,${ltype} ${ifile} s52
  mv s52 ${ifile}
}

function set_extension(){
  ifile=$1
  extension=$2
  test "${ifile}" = "${land_mask}" && land_mask=${ifile}.${extension}
  mv ${ifile} ${ifile}.${extension}
}


#--------------------------------------------------------------------------
#+ MAIN
#--------------------------------------------------------------------------
timemsg BEGIN

read iyy imm trash <<< `echo $sformateddate | tr '_T:-' '    '`
read fyy fmm trash <<< `echo $eformateddate | tr '_T:-' '    '`
tmpDir=`pwd`/BCdata
mkdir -p ${tmpDir}
cd ${tmpDir}

write_static_files

project=$(bctable_header project)
model=$(bctable_header model)
experiment=$(bctable_header experiment)
extdata_path=$(bctable_header data_path)
interval=$(bctable_header interval)
noleap=$(bctable_header remove_leap)

if test -z "${extdata_path}"; then
  timemsg "No data_path in BCtable. Trying to find out where are the data ..."
  institute=$(bctable_header institute)
  activity=$(bctable_header activity)
  product=$(bctable_header product)
  case "$(echo ${project} | tr '[:upper:]' '[:lower:]')" in
    cmip5) extdata_path="${ESGF_DRS_BASEDIR}/${project}/${activity}/${product}/${institute}/${model}";;
    cmip6) extdata_path="${ESGF_DRS_BASEDIR}/${project}/${activity}/${institute}/${model}";;
  esac
fi
test -d "${extdata_path}" || errormsg "Non-existing data path: ${extdata_path}"
timemsg "Retrieving global data from: ${extdata_path}"

bctable_variables ${BCTABLE} | while read \
  variable gribcode ltype version run frequency realm table filter
do
  test "${variable:0:1}" = "#" && continue
  timemsg "Processing ${variable}"
  case "$(echo ${project} | tr '[:upper:]' '[:lower:]')" in
    cmip5) drsdir="${extdata_path}/${experiment}/${frequency}/${realm}/${table}/${run}/${version}";;
    cmip6) drsdir="${extdata_path}/${experiment}/${run}/${table}/*/gn/${version}";;
  esac
  filespec="${variable}_${table}_${model}_${experiment}_${run}*.nc"
  is_only_ic=0
  for year in $(seq $iyy $fyy)
  do
    timemsg "  YEAR: ${year}"
    apply_filter "${prefilter}|${filter}|${postfilter}"
    test "${is_only_ic}" -eq 1 && break
  done # year
done # while read

# Deleting dat files
rm ${tmpDir}/*.dat
timemsg END