#!/bin/bash -e
## ensmean and multiyear average 
#;DIAG_NORCPM; CASENAME: 
#;DIAG_NORCPM; CASEPREFIX: 
#;DIAG_NORCPM; LEADYEAR: 0 
#;DIAG_NORCPM; INITMONTH: 10
#;DIAG_NORCPM; INITDAY: 15
#;DIAG_NORCPM; OUTPUTFILE: 
#;DIAG_NORCPM; PLEVHPA:
#;NeedBeReplaced: YEARS, TIMETAG, VARNAME, COMPONENT, ARCHIVEPATH, LEADYEARS, LYTAG
source /cluster/projects/nn9039k/people/pgchiu/conda/miniconda3/etc/profile.d/conda.sh && conda activate py3env_v2   

## https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-statistics/what-is-a-z-score-what-is-a-p-value.htm
## 
archivepath=/cluster/shared/NS9039K/norcpm_ana_hindcast/archive/
#archivepath=/cluster/shared/NS9039K/norcpm_ana_hindcast/archive/
## caseprefix means case=$caseprefix$year$month$day
## keep it empty to use fixed case, like hist or analysis run
caseprefix=norcpm_ana_f09_tn14_decadal_hindcast_
#caseprefix=norcpm_ana_f09_tn14_decadal_hindcast_
## case can be empty, will retrive from $archivepath
case=norcpm_ana_f09_tn14_decadal_hindcast
#case=norcpm_ana_f09_tn14_decadal_hindcast
test -z "$case" && case=$(basename $archivepath) && archivepath=$(dirname $archivepath)
## lead year, hist and analysis run should be 0
## can be multple years, seprate with ','
leadyrs=1
lytag=1

workdir="work_$$"

## time range
years='1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2014'
years=$(echo $years | sed 's/,/ /g')
## timealias will show in file name, cannot be 4 character
timealias=1999-2013
month=10
day=15

var=TREFHT
plev_hPa=''

dfn='atm/hist/'

## check run or not
outputfile='ensmean_hind_TREFHT_ly1_1999-2013.nc'
#outputfile=''
test -z "$outputfile" && outputfile="ensmean_${case}_${var}_${levtag}leadyr${lytag}_${timealias}.nc"
test -f "$outputfile" && echo "$outputfile present, skip." && exit 0

## create output vertical levels (need imply plev_hPa)
if [ -z "$plev_hPa" ] ; then
    vrtarg=''
    levtag=''
else
    vrt=vrt_prs_${plev_hPa}hPa.nc
    test -f $vrt || \
        ncap2 -O -v \
            -s 'defdim("plev",1);plev[$plev]={'$plev_hPa'00};' \
            $vrt
    vrtarg="--vrt_fl=$vrt"
    levtag="plev${plev_hPa}_"
fi

## lead year 1: year+1, jan-dec
for year in $years; do
    ifilelist=''
    if [ ! -z "$caseprefix" ] ;then
        case1=$caseprefix$year$month$day
    fi
    path="$archivepath/$case1"
    echo -n "cal ensmean of ${case1} init ${year} "
    for leadyr in $(echo $leadyrs | tr ',' ' ');do
        echo -n ", ${yr} for $leadyr"
        yr=$(($year+$leadyr))
        t0=$(date +%s)
        for mm in $(seq -w 01 12) ; do 
            datafns="atm/hist/*.h0.${yr}-${mm}.nc"
            ifilelist="$ifilelist $(ls $path/${case1}_mem*/$datafns)"
            ## for nco, most of operations are need be from single file (and output to single file)
        done
        ## yearly mean
    done
    ncea -O -h -v ${var} ${ifilelist} -o ensmean_${case}_${var}_leadyr${leadyr}_i${year}.nc
    echo "Done in "$(($(date +%s)-${t0})) sec.
done

ncrcat -O -h -v ${var} ensmean_${case}_${var}_leadyr${leadyr}_i????.nc -o $outputfile 
rm ensmean_${case}_${var}_leadyr${leadyr}_i????.nc