! Silent mode:
!
\cancel mode verify
!
! Ferret GO tool
!
! Part of POMfer Plotting Package for Plotting the Results of OzPOM
!
! John Hunter,
! Antarctic Climate & Ecosystems Cooperative Research Centre,
! Hobart, Tasmania, Australia
!
! john.hunter@utas.edu.au
!
! contour_path_section.jnl
!
! Contours:
!
!    1. Section of "xyz" scalar model variable ,
!    2. Section of "xyzt" scalar model variable,
!
! along a path defined by a series of waypoints.
!
! ($1) ..... window number
! ($2) ..... variable name
! ($3) ..... filename for file of longitude and latitude of waypoints
! ($4) ..... time (days) for scalar or "n" if time-independent
!
! Note that the layer index is reversed by Ferret from the way POM defines it (1 at surface, kb at bottom).
!
! Initial version, JRH 25/11/2004
!
! Cancel all:
!
cancel data/all
cancel data_set/all
cancel memory/all
cancel symbol/all
cancel variable/all
set window/aspect=.75:axis `($1)`
!
! No Ferret logo:
!
cancel mode logo
!
! Input special data for these model results:
!
go special_data
!
! Check smoothing option:
!
if `smop ne ""` then
  if `smop ne "sm"` then
    message/continue "***** smop must be '' or 'sm' *****"
    exit
  endif
endif
!
! Change variable name if "t"
!
if `"($2)" ne "t"` then
  if `"($2)" eq "psi"` then
    define symbol var = ($2)
  elif `"($2)" eq "psi_o"` then
    define symbol var = ($2)
  else
    define symbol var = ($2)`smop`
    define symbol filename_var = ($2)`smop`
  endif
else
  if `smop ne "sm"` then
    define symbol var = 'T'
    define symbol filename_var = t
  else
    define symbol var = TSM
    define symbol filename_var = tsm
  endif
endif
!
set data_set/format=cdf "`nc_file`"
!
! Find size of arrays (use salinity as typical 3-D array):
!
let im = `s,return=isize`
let jm = `s,return=jsize`
let kb = `s,return=ksize`
!
! Get prefix for metafile:
!
let slash = strrindex("`nc_file`","/")
let dot = strrindex("`nc_file`",".")
let len = dot - slash - 1
!
let prefix = substring("`nc_file`",`slash+1`,`len`)
!
set mode metafile path_section_`prefix`_($filename_var)_($3)_($4).plt
!
! Make plot title and define bias:
!
let units = "`($var),return=unit`"
!
if `"($2)" eq "t"` then
  if `smop ne "sm"` then
    let title = "Temperature (`units`)"
  else
    let title = "Smoothed temperature (`units`)"
  endif
  let bias = tbias
  let levels = tlevels
  let keymod = tkeymod
elif `"($2)" eq "s"` then
  if `smop ne "sm"` then
    let title = "Salinity (`units`)"
  else
    let title = "Smoothed salinity (`units`)"
  endif
  let bias = sbias
  let levels = slevels
  let keymod = skeymod
else
  let title = "($var), `($var),return=title` (`units`)"
  let bias = obias
  let levels = olevels
  let keymod = okeymod
endif
!
! Read waypoints (longitude, latitude; degrees):
!
define axis/z=1:20480:1 way_axis
define grid/z=way_axis way_grid
set data/ez/grid=way_grid/var="long_way,lat_way" "($3)"
!
let east_way = (long_way[d=2]-long0)*d2r*re*coslatscale
let north_way = (lat_way[d=2]-lat0)*d2r*re
!
! Number of waypoints:
!
let n_waypoints = `east_way,return=ksize`
!
! Generate files of indices of waypoints:
!
list/quiet/nohead/clobber/file=i_way/format=(a22) {"i-indices of waypoints"}
list/quiet/nohead/clobber/file=j_way/format=(a22) {"j-indices of waypoints"}
!
! Find the grid point with a minimum distance from each waypoint:
!
repeat/k=1:`n_waypoints` (let fun = int(((east_e[d=1]-east_way[d=2])^2+(north_e[d=1]-north_way[d=2])^2)^0.5+0.5);\
  let min_val = fun[i=@min,j=@min];\
  let idum = fun[i=@loc:`min_val`];let imin = idum[j=@max];\
  let jdum = fun[j=@loc:`min_val`];let jmin = jdum[i=@max];\
  list/quiet/nohead/append/file=i_way/format=(f5.0) `imin`;\
  list/quiet/nohead/append/file=j_way/format=(f5.0) `jmin`)
!
! Read back indices of waypoints:
!
! set data/ez/grid=way_grid/var="i_way"/skip=1 "i_way"
! set data/ez/grid=way_grid/var="j_way"/skip=1 "j_way"
!
set data/ez/var="i_way"/skip=1 "i_way"
set data/ez/var="j_way"/skip=1 "j_way"
!
! Interpolate to 10 times as many points:
!
define axis/x=1:`n_waypoints`:0.1 fine_axis
!
let i_fine = i_way[d=3,gx=fine_axis@lin]
let j_fine = j_way[d=4,gx=fine_axis@lin]
!
! Calculate horizontal coordinate:
!
let east_sample = samplexy(east_e[d=1],i_fine,j_fine)
let north_sample = samplexy(north_e[d=1],i_fine,j_fine)
!
! Distances between adjacent points:
!
let dum1 = ((east_sample[x=@shf:-1]-east_sample)^2+(north_sample[x=@shf:-1]-north_sample)^2)^0.5
let incr_dist = if missing(dum1,-9999) eq -9999 then 0 else dum1
!
! Integrate distance:
!
let dist_sample = (zz[d=1,k=2:`kb`]/zz[d=1,k=2:`kb`]) * incr_dist[x=@iin]
!
! Interpolate depth:
!
let elev_mask = fsm[d=1]
set variable/bad = 0 elev_mask
!
let dum2 = h[d=1]
set variable/bad = 1 dum2
if `"($4)" ne "n"` then
  if `smop ne "sm"` then
    let elev = elb[d=1,t=($4)]*elev_mask
  else
    let elev = elbsm[d=1,t=($4)]*elev_mask
  endif
else
!
! If no time-dependence, assume free surface at -e_atmos[l=1]:
!
  let elev = (0-e_atmos[d=1,l=1])*elev_mask
!
endif
! 
let dep = zz[d=1,k=2:`kb`] * (dum2+elev)+elev
let dep_1 = z[d=1,g=zz] * (dum2+elev)+elev
!
let dep_sample = samplexy(dep,i_fine,j_fine)
let dep_1_sample = samplexy(dep_1,i_fine,j_fine)
!
! Approximate vertical axis length (note "k=2:2:1" instead of just "k=2" to fix an apparent bug):
!
let axis_vert = `dep_sample[i=@min,k=2:2:1]`
!
! Interpolate scalar variable along track:
!
if `"($4)" ne "n"` then
  let dum3 = samplexy(($var)[d=1,k=2:`kb`,t=($4)],i_fine,j_fine)
else
  let dum3 = samplexy(($var)[d=1,k=2:`kb`],i_fine,j_fine)
endif
!
! Mask out any scalar which does not have appropriate coordinates:
!
let var_sample = if missing(dep_sample,9999) ne 9999 then dum3 else (-1.e+34)
!
!.............................................................................
!
! Plot title:
!
if `"($4)" ne "n"` then
!
! Time-dependent:
!
  let plot_title = "`title`, after `($4)` days, from track file ($3)"
!
else
!
! Not time-dependent:
!
  let plot_title = "`title`, from track file ($3)"
!
endif
!
!.............................................................................
!
! Plot:
!
contour/d=1/fill/title="`plot_title`"/levels=`levels`/vlimits=`axis_vert`:0/set_up var_sample+bias,dist_sample,dep_sample
!
! Key labelling modifications:
!
ppl shakey,`keymod`
ppl fill
!
!.............................................................................
!
! Plot boundaries:
!
! Upper and lower:
!
plot/overlay/nolabels/vs/line/color=black/thick=2 dist_sample[k=2],dep_1_sample[k=1]
plot/overlay/nolabels/vs/line/color=black/thick=2 dist_sample[k=2],dep_1_sample[k=`kb`]
!
!.............................................................................
!
! Verify mode:
!
set mode/last verify