#! /bin/csh

module load cray-python/3.6.1.1

# Directory where finished files will be stored
set SMBdatahome = /data/matlar/hirham/smb_offline/Data/CARRA_DataALLVAR
# Directory where intermediate calculations are done. Is emptied along the way
set dir_temp    = /data/matlar/hirham/smb_offline/Data/Python-CDS/TEMPDIRS/TEMP`date +%Y%m%d-%H%M%S`
#set dir_temp    = /data/pla/hirham/smb_offline/Data/Python-CDS/TEMPDIRS/TEMP20220112-145612
mkdir -p $SMBdatahome
mkdir -p $dir_temp

# SHift date tool:
set shift_date = /home/pla/tools/shift_date

#Which CDO to use:
set cdo     = "/home/pla/tools/cdo-1.6.6-static -s"
#set cdo     = "/netapp/dmisys/swrepo/cdo/1.9.5-0/bin/cdo -s"
#Which NCO to use:
set ncodist = /netapp/dmisys/swrepo/nco/4.6.0-1/GENERIC/bin
setenv LD_LIBRARY_PATH /netapp/dmisys/tl/temp_xc30_libs/opt/cray/netcdf/4.3.3.1/gnu/4.9/lib:/netapp/dmisys/tl/temp_xc30_libs/opt/cray/hdf5/1.8.14/gnu/4.9/lib


# Python script to download variables from CDS
set cds_download_py = /data/matlar/hirham/smb_offline/Data/Python-CDS/Download-CDS_v2.py

#set start = 20010427 # Included
#set end   = 20030101 # Not included

set start = 20181001 # Included
set end   = 20190101 # Not included

#set start = 20030101 # Included
#set end   = 20050101 # Not included

#set start = 20050101 # Included
#set end   = 20070101 # Not included

#set start = 20070101 # Included
#set end   = 20090101 # Not included

#set start = 20090101 # Included
#set end   = 20110101 # Not included

#set start = 20110101 # Included
#set end   = 20130101 # Not included

#set start = 20130101 # Included
#set end   = 20150101 # Not included

#set start = 20150101 # Included
#set end   = 20170101 # Not included

#set start = 20170101 # Included
#set end   = 20190101 # Not included

#set start = 20190101 # Included
#set end   = 20210101 # Not included

cd $dir_temp

set now = $start
while ($now < $end)
    set Y = `date --date=$now +%Y`
    set M = `date --date=$now +%m`
    set D = `date --date=$now +%d`
    echo "  Date "$Y,$M,$D
    # Find yesterday
    set yesterday = `$shift_date $now -1`
    set Yy = `date --date=$yesterday +%Y`
    set My = `date --date=$yesterday +%m`
    set Dy = `date --date=$yesterday +%d`
    # FOr cleanup along the way, find day before yesterday
    set dbyesterday = `$shift_date $now -2`
    set Yyy = `date --date=$dbyesterday +%Y`
    set Myy = `date --date=$dbyesterday +%m`
    set Dyy = `date --date=$dbyesterday +%d`

    set outdir = $SMBdatahome/$Y/$M/$D
    mkdir -p $outdir

   # File with time stamp hh shall include accumulated
   # fluxes from hh-3hr to hh. So:
   #         left        -  right
   # 00-03: 12yesterday:15  - 12yesterday:12
   # 03-06: 12yesterday:18  - 12yesterday:15
   # 06-09: 00today:09      - 00today:06
   # 09-12: 00today:12      - 00today:09
   # 12-15: 00today:15      - 00today:12
   # 15-18: 00today:18      - 00today:15
   # 18-21: 12today:09      - 12today:06
   # 21-00: 12today:12      - 12today:09

   # Times to make files for each day
   set times = ( 00-03   03-06   06-09   09-12   12-15   15-18   18-21   21-00 )
   set tstamp= ( 00      03      06      09      12      15      18      21     )
   # Dates of forecast start
   set Yf    = ( $Yy     $Yy     $Y      $Y      $Y      $Y      $Y      $Y)
   set Mf    = ( $My     $My     $M      $M      $M      $M      $M      $M)
   set Df    = ( $Dy     $Dy     $D      $D      $D      $D      $D      $D)
   # Times of forecast start
   set Tf    = ('12:00' '12:00' '00:00' '00:00' '00:00' '00:00' '12:00' '12:00')
   # Lead-times left and right files
   set Ll    = ( 15      18      09      12      15      18      09      12)
   set Lr    = ( 12      15      06      09      12      15      06      09 )
  
   # If it doesn't exist yet, download yesterday's 12:00 file, and split it:
   if ( ! -e  CARRA_{$Yy}{$My}{$Dy}_12:00.nc ) then
      python $cds_download_py $Yy $My $Dy '12:00' 
      $cdo splithour CARRA_{$Yy}{$My}{$Dy}_12:00.nc CARRA_{$Yy}{$My}{$Dy}_12:00+
   endif
   # Download today's  two files, and split them
   if ( ! -e  CARRA_{$Y}{$M}{$D}_00:00.nc ) then
      python $cds_download_py $Y $M $D '00:00'  
      $cdo splithour CARRA_{$Y}{$M}{$D}_00:00.nc CARRA_{$Y}{$M}{$D}_00:00+
   endif
   if ( ! -e  CARRA_{$Y}{$M}{$D}_12:00.nc ) then
      python $cds_download_py $Y $M $D '12:00'  
      $cdo splithour CARRA_{$Y}{$M}{$D}_12:00.nc CARRA_{$Y}{$M}{$D}_12:00+
   endif
   

   set tn = 0 # Counter for times
   foreach time ($times)
     @ tn += 1
     # Subtract the two
     set diff_time = {$Yf[$tn]}{$Mf[$tn]}{$Df[$tn]}_{$times[$tn]}
     $cdo sub \
	CARRA_{$Yf[$tn]}{$Mf[$tn]}{$Df[$tn]}_{$Tf[$tn]}+{$Ll[$tn]}.nc \
	CARRA_{$Yf[$tn]}{$Mf[$tn]}{$Df[$tn]}_{$Tf[$tn]}+{$Lr[$tn]}.nc \
      diff_{$diff_time}.nc
     
     # Calculate mean 2m temp as average over the two slices:
     $cdo mulc,0.5 -selname,t2m -add \
	CARRA_{$Yf[$tn]}{$Mf[$tn]}{$Df[$tn]}_{$Tf[$tn]}+{$Ll[$tn]}.nc \
	CARRA_{$Yf[$tn]}{$Mf[$tn]}{$Df[$tn]}_{$Tf[$tn]}+{$Lr[$tn]}.nc \
      CARRA_2m_temperature_{$diff_time}.nc
     
     # Extract energy fluxes and divide by 3 hr (= 10800 s) to go from J/m2/3hr to W/m2:
     $cdo chname,slhf,ahfl   -divc,10800 -selname,slhf diff_{$diff_time}.nc CARRA_ahfl_{$diff_time}.nc
     $cdo chname,sshf,ahfs   -divc,10800 -selname,sshf diff_{$diff_time}.nc CARRA_ahfs_{$diff_time}.nc
     $cdo chname,ssrd,dswrad -divc,10800 -selname,ssrd diff_{$diff_time}.nc CARRA_dswrad_{$diff_time}.nc
     $cdo chname,strd,dlwrad -divc,10800 -selname,strd diff_{$diff_time}.nc CARRA_dlwrad_{$diff_time}.nc

     # Extract mass fluxes
     # Evaporation and sublimation: Multiply by -1 to get positive downwards (as in HH and SMB-code):
     #$cdo chname,eva,evspsbl -mulc,-1    -selname,eva  diff_{$diff_time}.nc CARRA_evspsbl_{$diff_time}.nc
     #Temporary fix to get evap+subl:
     # AS an intermediate fix, calculate total evspsbl (evap+subl) from
    # the total latent heat flux (ahfl) assuming that it is all due to
    # sublimation, ie., evspsbl =  ahfl /Ls * 10800 s, where Ls = 2834000 J/kgCARRA_evspsbl_{$diff_time}.nc
    # Sign convention:  evspsbl should be positive downwards (as in HH and SMB-code)
    # so we keep the sign of ahfl which is positive when mass flux is down (deposition).
    # Multiply by 10800 s and divide by Ls(J/kg)  to go from W/m2 to W/m2 * s /(J/kg) = kg/m2
    # That is, the total vapor mass flux over 3 hours, given a ahfl flux measured in W/m2.
     cdo -expr,'evspsbl = ahfl/2834000* 10800' CARRA_ahfl_{$diff_time}.nc CARRA_evspsbl_{$diff_time}.nc
     $ncodist/ncatted -O -a units,evspsbl,c,c,"kg m**-2" -a long_name,evspsbl,c,c,"evap+subl" CARRA_evspsbl_{$diff_time}.nc 


     # Get total precip
     $cdo selname,tp diff_{$diff_time}.nc CARRA_tp_{$diff_time}.nc

     # Tentative calculation of snow fraction from temperature field (to be replaced by actual snow/rain later)
     # Snow fraction is 1 below -2, 0 above +4, and ramped linearly between -2 and 4:
     # e.g. Feiccabrino and Lundberg 2009
     $cdo -expr,'tc=t2m-273.15'                 CARRA_2m_temperature_{$diff_time}.nc CARRA_tc_{$diff_time}.nc
     $cdo -expr,'S0= 1-(tc-(-2))/(4-(-2))'     CARRA_tc_{$diff_time}.nc CARRA_S0_{$diff_time}.nc
     # In the following we need a file with 1's and one with 0's:
     $cdo -expr,'S0= tc-tc +1 '     CARRA_tc_{$diff_time}.nc CARRA_ONES_{$diff_time}.nc
     $cdo -expr,'S0= tc-tc    '     CARRA_tc_{$diff_time}.nc CARRA_ZEROS_{$diff_time}.nc
     # Find places where S0<1
     $cdo ltc,1 CARRA_S0_{$diff_time}.nc CARRA_S0_lt1_{$diff_time}.nc
     # Make S1: Where S0<1, set to S0, elsewhere set to 1
     $cdo ifthenelse CARRA_S0_lt1_{$diff_time}.nc CARRA_S0_{$diff_time}.nc CARRA_ONES_{$diff_time}.nc CARRA_S1_{$diff_time}.nc
     # Find where S1>0
     $cdo gtc,0 CARRA_S1_{$diff_time}.nc CARRA_S1_gt0_{$diff_time}.nc
     # Where S1>0, set to S1, elsewhere set to 0
     $cdo ifthenelse CARRA_S1_gt0_{$diff_time}.nc CARRA_S1_{$diff_time}.nc CARRA_ZEROS_{$diff_time}.nc CARRA_snowfracS_{$diff_time}.nc
     # Rename result to "snowfrac"
     $cdo chname,S0,snowfrac CARRA_snowfracS_{$diff_time}.nc CARRA_snowfrac_{$diff_time}.nc
     # Let rainfrac = 1 - snowfrac
     $cdo selname,rainfrac \
          -expr,'rainfrac = 1 - snowfrac' CARRA_snowfrac_{$diff_time}.nc CARRA_rainfrac_{$diff_time}.nc
     # Multiply total precip by snowfrac and rainfrac to get snow and rain
     $cdo chname,tp,snfall   -mul CARRA_tp_{$diff_time}.nc CARRA_snowfrac_{$diff_time}.nc CARRA_snfall_{$diff_time}.nc
     $cdo chname,tp,rainfall -mul CARRA_tp_{$diff_time}.nc CARRA_rainfrac_{$diff_time}.nc CARRA_rainfall_{$diff_time}.nc
 
     # Merge all varaibles into one file in the correct final destination:
     $cdo merge CARRA_{ahfl,ahfs,dswrad,dlwrad,evspsbl,snfall,rainfall,2m_temperature}_{$diff_time}.nc ALLVARS_{$diff_time}.nc
     # Set time stamp and put in final destination
     $cdo setdate,${Y}-${M}-${D} -settime,{$tstamp[$tn]}:00:00 ALLVARS_{$diff_time}.nc $outdir/ALLVARS_{$times[$tn]}.nc

   end
   rm -f *{$Yyy}{$Myy}{$Dyy}*.nc 

   set now = `$shift_date $now 1`
end #while over days

