;----------------------------------------------------------------------
; This script creates plots from "wind.ascii" and "temp.ascii",
; downloaded from:
;
; http://www.sparc-climate.org/data-center/data-access/reference-climatologies/randels-climatologies/temperature-wind-climatology/
;
; The data and README file can be found at:
;    ftp://sparc-ftp1.ceda.ac.uk/sparc/ref_clim/randel/temp_wind/
;----------------------------------------------------------------------
; This script creates a very customized set of plots. See
; plot_sparc_data_simple.ncl for a very basic set of plots.
;----------------------------------------------------------------------
begin

;---Set up dimension information based on README mentioned above
  ntim  = 12   ; JAN, FEB, ..., DEC
  nlat  = 41   
  ntlev = 33
  nwlev = 46

;----------------------------------------------------------------------
; Create time, lat, level coordinate arrays for 3D arrays
;----------------------------------------------------------------------
  time      = ispan(1,ntim,1)
  lat       = ispan(-80,80,4)                 ; nlat
  temp_pres = 10^(3.-ispan(0,ntlev-1,1)/6.)
  wind_pres = 10^(3.-ispan(0,nwlev-1,1)/6.)

  time!0          = "time"
  time&time       = time
  time@long_name  = "months"

  lat!0         = "lat"
  lat&lat       = lat
  lat@long_name = "latitude"
  lat@units     = "degrees_north"

  temp_pres!0         = "lev"
  temp_pres&lev       = temp_pres
  temp_pres@long_name = "pressure"
  temp_pres@units     = "hPa"

  wind_pres!0         = "lev"
  wind_pres&lev       = wind_pres
  wind_pres@long_name = "pressure"
  wind_pres@units     = "hPa"

;----------------------------------------------------------------------
; Read ASCII files containing temperature and wind data
;----------------------------------------------------------------------
  temp_filename = "temp.ascii"
  wind_filename = "wind.ascii"

  temp_dims     = (/ntlev,nlat,ntim/)     ; This is in Fortran order
  wind_dims     = (/nwlev,nlat,ntim/)
  
  temp = asciiread(temp_filename,temp_dims,"float")
  wind = asciiread(wind_filename,wind_dims,"float")

;----------------------------------------------------------------------
; This section assigns metadata as described in the README file
;----------------------------------------------------------------------
  temp@long_name  = "zonal mean temperature"
  temp@units      = "degK"
  wind@long_name  = "zonal mean zonal wind"
  wind@_FillValue = 1e36    ; missing value
  wind@units      = "m/s"

  temp!0    = "lev"
  temp!1    = "lat"
  temp!2    = "time"
  temp&lev  = temp_pres
  temp&lat  = lat
  temp&time = time

  wind!0    = "lev"
  wind!1    = "lat"
  wind!2    = "time"
  wind&lev  = wind_pres
  wind&lat  = lat
  wind&time = time

;----------------------------------------------------------------------
; Debug prints to make sure everything looks ok.
;----------------------------------------------------------------------
  printVarSummary(temp)
  printMinMax(temp,0)

  printVarSummary(wind)
  printMinMax(wind,0)

;----------------------------------------------------------------------
; Start the graphics
;----------------------------------------------------------------------
  wtype = "png"           ; x11, ps, pdf, svg
  if(wtype.eq."png")
    wtype@wkWidth  = 1500
    wtype@wkHeight = 1500
  end if
  wks = gsn_open_wks(wtype,"sparc_custom")

;---Set some graphical resources (options) common for both wind and temp plots
  res                      = True

  res@gsnDraw              = False    ; don't draw individual plots
  res@gsnFrame             = False    ; don't advance frame

  res@trYReverse           = True     ; reverse Y axis

;---Customize contours
  res@cnFillOn             = True     ; turn on contour fill
  res@cnLineLabelsOn       = False    ; turn off line labels
  res@cnInfoLabelOn        = False    ; turn off info label
  res@cnFillPalette        = "MPL_jet" ; change the color map
  res@lbLabelBarOn         = False    ; will add labelbar in panel
  res@cnLevelSelectionMode = "ManualLevels"

;---Various titles
  res@tiYAxisString        = temp&lev@long_name + " (" + temp&lev@units + ")"
  res@tiXAxisString        = lat@long_name
  res@gsnRightString       = ""
  res@gsnLeftString        = ""

;---Customize tickmarks
  res@tmXTOn               = False       ; turn off top tickmarks
  res@tmYROn               = False       ; turn off right tickmarks
  res@tmYLMode             = "Explicit"  ; will explicitly label left axis (later)

;---Now create individual resource lists for each variable
  tres = res
  wres = res

;---Set explicit labels on left Y axis
  tres@tmYLValues          = (/0.005,0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,\
                               20,50,100,200,500,1000/)
  tres@tmYLLabels          = "" + tres@tmYLValues

  wres@tmYLValues          = (/5e-5,1e-4,0.0002,0.0005,0.001,0.002,0.005,\
                               0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,\
                               100,200,500,1000/)
  wres@tmYLLabels          = "" + wres@tmYLValues

;---Fix the contour levels for both plots
  temp_mnmxint = nice_mnmxintvl( min(temp), max(temp), 18, False)
  wind_mnmxint = nice_mnmxintvl( min(wind), max(wind), 18, False)

  tres@cnMinLevelValF      = temp_mnmxint(0)
  tres@cnMaxLevelValF      = temp_mnmxint(1)
  tres@cnLevelSpacingF     = temp_mnmxint(2)

  wres@cnMinLevelValF       = wind_mnmxint(0)
  wres@cnMaxLevelValF       = wind_mnmxint(1)
  wres@cnLevelSpacingF      = wind_mnmxint(2)

;---Loop through all time steps and create a temp/wind plot for each
  plot_temp = new(ntim,graphic)    ; array to hold each plot
  plot_wind = new(ntim,graphic)
  do nt=0,ntim-1
    plot_temp(nt) = gsn_csm_contour(wks,temp(:,:,nt),tres)
    plot_wind(nt) = gsn_csm_contour(wks,wind(:,:,nt),wres)
  end do

;---Panel both sets of plots
  pres                       = True
  pres@gsnPanelLabelBar      = True
  pres@gsnMaximize           = True

  pres@gsnPanelFigureStrings            = str_upper(month_name(0))  ; "JAN", "FEB", ...
  pres@gsnPanelFigureStringsFontHeightF = 0.01

  pres@txString         = temp_filename + " - " + temp@long_name + " (" + temp@units + ")"
  gsn_panel(wks,plot_temp,(/3,4/),pres)

  pres@txString         = wind_filename + " - " + wind@long_name + " (" + wind@units + ")"
  gsn_panel(wks,plot_wind,(/3,4/),pres)
end