#!/bin/bash -e
#;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
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/archive/norcpm_ana_f09_tn14
## caseprefix means case=$caseprefix$year$month$day
## keep it empty to use fixed case, like hist or analysis run
caseprefix=
## case can be empty, will retrive from $archivepath
case=
test -z "$case" && case=$(basename $archivepath) && archivepath=$(dirname $archivepath)
## lead year, hist and analysis run should be 0
leadyr=0 

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_analysis_TREFHT_1999-2013.nc'
test -z "$outputfile" && outputfile="ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${timealias}.nc"
test -f "$outputfile" && echo "$outputfile present, skip." && exit 0

## goto workdir
mkdir -p $workdir
cd $workdir

## 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
    if [ ! -z "$caseprefix" ] ;then
        case=$caseprefix$year$month$day
    fi
    path="$archivepath/$case"
    yr=$(($year+$leadyr))
    echo -n "cal ensmean of ${case} ${yr} for lead year $leadyr of ${year}..."
    t0=$(date +%s)
    for mm in $(seq -w 01 12) ; do 
        datafns="atm/hist/*.h0.${yr}-${mm}.nc"
        for mem in $(seq -w 01 10) ; do 
            fn=$(ls $path/${case}_mem${mem}/$datafns)
            basefn=$(basename ${fn})
            ## interpolate to plev
            if [ -z "$plev_hPa" ] ;then
                ncks    -v $var ${fn} -o ${levtag}${var}_${basefn}
            else
                ncremap $vrtarg  -v $var -i ${fn} -o ${levtag}${var}_${basefn} > ncremap.log
            fi
        done
        ## for nco, most of operations are need be from single file (and output to single file)
        ncecat -u member -v ${var} ${levtag}${var}_*.nc mems_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc
        rm -f ${levtag}${var}_*.nc
        ## ensemble mean and stddev
        ncwa -O -a member  mems_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc \
                        ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc 
        rm -f     mems_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc  

        #ncks -O --mk_rec_dmn time ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc \
        #                       -o ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-${mm}.nc 
    done
    ## yearly mean
    nces -O -v $var \
               ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}-??.nc \
            -o ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}.nc


    ## set time axis as unlimit, for combine years to 1 file
    ncks -O --mk_rec_dmn time ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}.nc  \
                           -o ensmean_${case}_${var}_${levtag}leadyr${leadyr}_${yr}.nc
    echo "Done in "$(($(date +%s)-${t0})) sec.
done

## combine into 1 file with t,y,x to correlate it with obs t,y,x data
##                                              remap needed
ncrcat -O ensmean_*_${var}_${levtag}leadyr${leadyr}_????.nc $outputfile
rm -f ensmean_*_${var}_${levtag}leadyr${leadyr}_????.nc 
rm -f ensmean_*_${var}_${levtag}leadyr${leadyr}_????-??.nc 


## end, back
mv *.nc ..
cd ..
rmdir $workdir