## This program for a data have 6 hour time resolution


# Set a series of shell variables:
# File containing Streamfunction at 850 (PSI)
psifile="./psi.nc"

# File containing mean sea level pressure (msl)
pslfile="./msl.nc"

# File containing relative humidity at 850 hPa (RH)
rhfile="./RH850.nc"

# File containing surface wind field  (u10, v10)
uvfile="./uvsurf.nc"

# File containing surface geopotential (zs)
zsfile="./geopotential_surface.nc"

# File containing land-sea fraction (lsm)
lsmfile="./land_sea.nc"

# Candidate file
candidatefile="./LPS_candidates.txt"

# Output file
outfile="./LPS_output.txt"

# Track file
trackfile="./LPS_track.txt"


# Detect LPS:
#   1) Search for candidates as minima in the streamfunction field
#   2) Merge candidates within 5 degrees
#   3) Output streamfunction minima at the candidate point,
#      sea level pressure minima within 3 degrees of the candidate point,
#      wind speed maxima within 3 degrees of candidate, 
#      relative humidity averaged over the area within 3 degrees of candidate,
#      surface geopotential maximum within 2 degrees of candidate, 
#      land-sea ratio at candidate point, and
#      minimum sea level pressure location

$path-to-TEMPEST-bin/DetectNodes  
    --in_data "$psifile;$pslfile;$rhfile;$uvfile;$zsfile;$lsmfile"
    --out $candidatefile
    --searchbymin "PSI" --mergedist 5.0 
    --closedcontourcmd "PSI,12.5e5,10,0"
    --outputcmd "PSI,min,0;msl,min,3;RH,avg,3;zs,max,2;lsm,max,0;
     _VECMAG(u10,v10),max,3;msl,minix,3.0" 

#Stitch candidates to make tracks:
#   1) With a maximum distance of 3.0 degrees between in-sequence candidates
#   2) With a minimum length of 5 time points
#       (ensures track lasts longer than 1 day in six-hourly data)
#   3) With a maximum gap size of 2 time points (12 hours)
#   4) With the surface geopotential maximum less than 8000 m^2 s^{-2}
#       for at least 4 time steps along the track
#   5) With the relative humidity higher than 85% for at least 4 time steps 
#       along the track.

$path-to-TEMPEST-bin/StitchNodes --in $candidatefile --out $outfile 
    --format "i,j,lon,lat,strf850,slp,rh,zs,lsm,sp,minmslix" 
    --range 3.0 --minlength 5 --maxgap 2 
    --threshold "zhi,<=,8000.,4;rh,>=,85,4"

############################################################################
## for 1 hour time interval data 	                                   #
## $path-to-TEMPEST-bin/StitchNodes --in $candidatefile --out $outfile     #
##    --format "i,j,lon,lat,strf850,slp,rh,zhi,lsm,sp,minmslix"            #
##    --range 3.0 --minlength 30 --maxgap 12                               #
##    --threshold "zhi,<=,8000.,24;rh,>=,85,24"                            #
#                                                                          #
## for 3 hour time interval data         		                   #
## $path-to-TEMPEST-bin/StitchNodes --in $candidatefile --out $outfile     #
##    --format "i,j,lon,lat,strf850,slp,rh,zhi,lsm,sp,minmslix"            #
##    --range 3.0 --minlength 10 --maxgap 4                                #
##    --threshold "zhi,<=,8000.,8;rh,>=,85,8"                              #
############################################################################

# Calculate some derived quantities using track data:
# count the closed contours of MSLP, calculate ACEPSL, ACE, PDI, IKE
# ACEPSL = 1e-4 * 3.92 * (1.94384 * (1016 - PSL)^0.644)^2 if PSL < 1016.  Here PSL is the  
  minimum surface pressure at the point of detection in hPa.  This formula 
  leverages the TC wind-pressure relationship to provide a proxy for ACE.
# ACE = 1e-4 * (1.94384 kt/(m/s) * Umax)^2, where Umax is the maximum 
  instantaneous wind speed magnitude within the three degree radius (in m/s).
# PDI = Umax^3, where Umax is the maximum instantaneous windspeed magnitude 
  within the three degree radius (in m/s)
# IKE = sum_i 0.5 * U_i^2 * area_i, where the sum is taken over all grid boxes  
  within the given radius, U_i is the wind magnitude in that grid cell,  
  and area_i is the area of the grid cell (in m^2).
 
$path-to-TEMPEST-bin/NodeFileEditor --in_file $outfile --out_file $trackfile 
    --in_data "$pslfile;$uvfile" 
    --in_fmt "lon,lat,strf850,slp,rh,zs,lsm,sp,minmslix" 
    --calculate "deltaslp=max_closed_contour_delta(msl,10,minmslix);
        acepsl=eval_acepsl(msl,10.0);
        ace=eval_ace(u10,v10,3.0);
        pdi=eval_pdi(u10,v10,3.0);
        ike=eval_ike(u10,v10,3.0)" 
    --out_fmt "lon,lat,strf850,slp,deltaslp,sp,rh,zs,lsm,acepsl,ace,pdi,ike"

