*
* Demonstrates the calculation of streamfunction and velocity
* potential using the following OpenGrADS extensions:
*
*     - libbjt (advection functions)
*     - fish (poison solver)
*
* Restrictions:
*
*    - global and uniform grid (use re() if grid is not uniform)
*    - first and latidudinal gridpoints cannot be at the poles
*
* This script is in the public domain.

function utfish(lev)

  if ( lev='' )
       lev = 200
  endif

* GFS rotating forecasts
* ----------------------
 'sdfopen http://nomad2.ncep.noaa.gov:9090/dods/gfs/rotating/gfs_00z'

* Auxiliaty fields
* ----------------
  'one    = 1 + 0 * lat '
  'cosphi = cos(0.0174532925199433 * lat)'

* B.-J.'s advection functions: intrinsic functions hcurl/hdivg
* do not handle boundary conditions properly
* ------------------------------------------------------------
  'set lev ' lev
  'vor = - ( madvu(vgrd,one) - madvv(ugrd,cosphi) / cosphi )'
  'div = - ( madvu(ugrd,one) + madvv(vgrd,cosphi) / cosphi )'

*  Solve poisson equation
*  ----------------------
  'psi = fish(vor) '
  'chi = fish(div) '

*  Standardize arbritary constant by forcing psi/chi
*  to have zero global mean
*  -------------------------------------------------
  'psi = psi - aave(psi,global)'
  'chi = chi - aave(chi,global)'

* Rotational and divergent winds
* ------------------------------
  'upsi =   mvadv(one,psi)'
  'vpsi = - muadv(one,psi)'
  'uchi = - muadv(one,chi)'
  'vchi = - mvadv(one,chi)'

* Plotting
* --------
  'set gxout shaded'
  'display psi/1e7'
  'set gxout stream'
  'display upsi;vpsi'
  'draw title Streamfunction/Rotational Wind - 'lev' hPa'
  'cbarn'
  'printim psi.png'
  'clear'

  'display chi/1e6'
  'set gxout vector'
  'set cmin  2'
  'set cmax 20'
  'display skip(uchi,6,6);vchi'
  'draw title Velocity Potential/Divergent Wind - 'lev' hPa'
  'cbarn'
  'printim chi.png'

return