! 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
!
! plot_k_section.jnl
!
! Plots:
!
!    1. "xy" scalar model variable,
!    2. Specified time for "xyt" scalar model variable,
!    3. k-slice of "xyz" scalar model variable or
!    4. Specified time for k-slice of "xyzt" scalar model variable,
!
! with option of overlayed 2-D or 3-D vector field.
!
! ($1) ..... window number
! ($2) ..... plot type ("shade" or "contour")
! ($3) ..... variable name
! ($4) ..... "v" for velocity vectors, otherwise "n"
! ($5) ..... k-index for scalar or "n" if independent of z
! ($6) ..... k-index for vector or  "n" if independent of z or no vector
! ($7) ..... time (days) for scalar or "n" if time-independent
! ($8) ..... time (days) for vector or "n" if time-independent or no vector
!
! Note that the k index should be input as POM defines it (1 at surface, kb at bottom).
! This is reversed to the Ferret convention by this script prior to plotting.
!
! Initial version, JRH 25/11/2004
!
! Cancel all:
!
cancel data/all
cancel data_set/all
cancel memory/all
cancel symbol/all
cancel variable/all
!
! No Ferret logo:
!
cancel mode logo
!
! Check inputs (complicated forms are to solve parsing bug involving "and"):
!
if `"($2)" ne "shade"` then
  if `"($2)" ne "contour"` then
    message/continue "***** Second argument should be shade or contour *****"
    exit
  endif
endif
!
if `"($4)" ne "v"` then
  if `"($4)" ne "n"` then
    message/continue "***** Fourth argument should be v or n *****"
    exit
  endif
endif
!
if `"($7)" ne "n"` then
  if `"($8)" ne "n"` then
    if `"($7)" ne "($8)"` then
      message/continue "***** Times for scalar and vector data should be the same *****"
      exit
    endif
  endif
endif
!
if `"($5)" ne "n"` then
  if `"($6)" ne "n"` then
    if `"($5)" ne "($6)"` then
      message/continue "***** k-indices for scalar and vector data should be the same *****"
      exit
    endif
  endif
endif
!
if `"($4)" eq "v"` then
  if `"($8)" eq "n"` then
    message/continue "***** Vector must be time-dependent *****"
    exit
  endif
endif
!
! 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 `"($3)" ne "t"` then
  if `"($3)" eq "psi"` then
    define symbol var = ($3)
    define symbol filename_var = ($3)
  elif `"($3)" eq "psi_o"` then
    define symbol var = ($3)
    define symbol filename_var = ($3)
  else
    define symbol var = ($3)`smop`
    define symbol filename_var = ($3)`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`
!
! Find aspect ratio:
!
let ydx = `(north_v[x=@max,y=@max]-north_v[x=@min,y=@min])/(east_u[x=@max,y=@max]-east_u[x=@min,y=@min])`
!
! 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 k_section_`prefix`_($2)_($filename_var)_($4)_($5)_($6)_($7)_($8).plt
!
! Make plot title and define bias:
!
let units = "`($var),return=unit`"
!
if `"($3)" 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 `"($3)" 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
!
! Mask out land:
!
let var_mask = if fsm eq 1 then ($var) else (-1.e+34)
!
! Scalar type:
!
if `"($7)" ne "n"` then
!
! Time-dependent scalar:
!
  if `"($5)" ne "n"` then
!
! Vertically-varying:
!
    if `"($2)" eq "shade"` then
      let scalar_spec = "[i=1:`im-1`,j=1:`jm-1`,k=`kb-($5)+1`,t=($7)]"
    else
      let scalar_spec = "[k=`kb-($5)+1`,t=($7)]"
    endif
!
  else
!
! Not vertically-varying:
!
    if `"($2)" eq "shade"` then
      let scalar_spec = "[i=1:`im-1`,j=1:`jm-1`,t=($7)]"
    else
      let scalar_spec = "[t=($7)]"
    endif
!
  endif
!
else
!
! Not time-dependent scalar:
!
  if `"($5)" ne "n"` then
!
! Vertically-varying:
!
    if `"($2)" eq "shade"` then
      let scalar_spec = "[i=1:`im-1`,j=1:`jm-1`,k=`kb-($5)+1`]"
    else
      let scalar_spec = "[k=`kb-($5)+1`]"
    endif
!
  else
!
! Not vertically-varying:
!
    if `"($2)" eq "shade"` then
      let scalar_spec = "[i=1:`im-1`,j=1:`jm-1`]"
    else
      let scalar_spec = ""
    endif
!
  endif
!
endif
!
! Depth-variation of vector:
!
if `"($4)" eq "v"` then
!
  if `"($6)" ne "n"` then
!
! Vertically-varying:
!
    if `smop ne "sm"` then
      define symbol udata = u
      define symbol vdata = v
    else
      define symbol udata = usm
      define symbol vdata = vsm
    endif
    let vector_spec = "[i=1:`im`:`iskip`,j=1:`jm`:`jskip`,k=`kb-($6)+1`,t=($8)]"
    list kb
!
  else
!
! Not vertically-varying:
!
    if `smop ne "sm"` then
      define symbol udata = uab
      define symbol vdata = vab
    else
      define symbol udata = uabsm
      define symbol vdata = vabsm
    endif
    let vector_spec = "[i=1:`im`:`iskip`,j=1:`jm`:`jskip`,t=($8)]"
!
  endif
!
endif
!
! Plot title:
!
if `"($5)" eq "n"` then
  let p_layer = "($6)"
else
  let p_layer = "($5)"
endif
!
if `"($7)" eq "n"` then
  let p_time = "($8)"
else
  let p_time = "($7)"
endif
!
if `p_time ne "n"` then
!
! Time-dependent:
!
  if `p_layer ne "n"` then
!
! Vertically-varying:
!
    let plot_title = "`title`, level `p_layer`, after `p_time` days"
!
  else
!
! Not vertically-varying:
!
    let plot_title = "`title`, after `p_time` days"
!
  endif
!
else
!
! Not time-dependent:
!
  if `p_layer ne "n"` then
!
! Vertically-varying:
!
    let plot_title = "`title`, level `p_layer`"
!
  else
!
! Not vertically-varying:
!
    let plot_title = "`title`"
!
  endif
!
endif
!
! Plot:
!
! Scalar plot:
!
if `"($2)" eq "shade"` then
!
  set window/aspect=`ydx`:axis `($1)`
  shade/d=1/title="`plot_title`"/levels=`levels`/set_up var_mask`scalar_spec`+bias,east_c,north_c
!
! Key labelling modifications:
!
  ppl shakey,`keymod`
  ppl shade
!
else
!
  set window/aspect=`ydx`:axis `($1)`
  contour/d=1/fill/title="`plot_title`"/levels=`levels`/set_up var_mask`scalar_spec`+bias,east_e,north_e
!
! Key labelling modifications:
!
  ppl shakey,`keymod`
  ppl fill
!
endif
!
! Vector plot:
!
if `"($4)" eq "v"` then
!
! Regrid:
!
  define axis/x=1.5:`im+.5`:1 xe
  define axis/y=1.5:`jm+.5`:1 ye
!
! Rotate to east and north:
!
  let ($udata)e = ($udata)[gx=xe@lin]*cos(rot)-($vdata)[gy=ye@lin]*sin(rot)
  let ($vdata)e = ($udata)[gx=xe@lin]*sin(rot)+($vdata)[gy=ye@lin]*cos(rot)
!
  go my_poly_vectors east_e[i=1:`im`:`iskip`,j=1:`jm`:`jskip`] north_e[i=1:`im`:`iskip`,j=1:`jm`:`jskip`] ($udata)e`vector_spec` ($vdata)e`vector_spec` `vlength` "arrow"
!
! Vector key (note that it apparently has to be before the polygon call to work properly):
!
  go poly_arrow_key "","`vlength` m/s",black,line,`arrow_key_x`,`arrow_key_y`
!
  let fill = plt_x_arrow[i=1]/plt_x_arrow[i=1]
  polygon/over/noax/nolab/palette=black plt_x_arrow,plt_y_arrow,fill
!
endif
!
! Mask land etc:
!
file/var="longpol,latpol" "land_pol"
!
let eastpol = (longpol-long0)*d2r*re*coslatscale
let northpol = (latpol-lat0)*d2r*re
!
polygon/overlay/nolabels/fill/levels=(0,16,1) eastpol,northpol,5.5
plot/overlay/nolabels/vs/line=7 eastpol,northpol
!
! Ice front:
!
file/var="longfront,latfront" "AIS_front.dat"
!
let eastfront = (longfront-long0)*d2r*re*coslatscale
let northfront = (latfront-lat0)*d2r*re
!
plot/overlay/nolabels/vs/line/color=green/thick=2 eastfront,northfront
!
! Trial polar grid:
!
! file/var="eastgrid,northgrid" "polar_grid_2"
! plot/overlay/nolabels/vs/line=1 eastgrid,northgrid
!
! Track for section:
!
file/var="longtrack,lattrack" "waypoints"
!
let easttrack = (longtrack-long0)*d2r*re*coslatscale
let northtrack = (lattrack-lat0)*d2r*re
!
plot/overlay/nolabels/vs/line/color=black/thick=2 easttrack,northtrack
!
! Verify mode:
!
set mode/last verify