;*****************************************************************************************
; Script that reads .dat files from the axisymmetric version of Brian Tang's ASPECH model
; and plots vertical velocity.
;
; Created by: Rosimar Rios Berrios (15 July 2014)
; Edited by: Josh Alland (22 July 2014)
;*****************************************************************************************

load "/nfs/kc11/jalland/Downdrafts_Project/load.ncl"

path2files = "/kc11/btang/aspechout/xdef100/"

;******OBTAIN RADIUS COORDINATES FROM STAGGERED C-GRID************************************
  rcoord = asciiread(path2files+"r.dat",-1,"float")
  dims   = dimsizes(rcoord)
  r      = rcoord(0:dims-1:2)
  rs     = rcoord(1:dims-1:2)    ;---Vertical velocity calculated on rs coordinates
  rs_km  = rs/1000.

;******OBTAIN HEIGHT COORDINATES FROM STAGGERED C-GRID************************************
  zcoord = asciiread(path2files+"z.dat",-1,"float")
  dims   = dimsizes(zcoord)
  z      = zcoord(0:dims-1:2)     ;---Vertical velocity calculated on z coordinates
  z_km   = z/1000.
  zs     = zcoord(1:dims-1:2) 
  zs_km  = zs/1000. 

;*****OBTAIN VERTICAL VELOCITY DATA FOR EACH TIME*****************************************
  d1   = dimsizes(rs)
  d2   = dimsizes(z)
;  var  = "w360.dat"
   var = "w3000.dat"
  w    = fbinread(path2files+var, (/d1,d2/),"float")

  w!0      = "radius"
  w!1      = "height"
  w&radius = rs_km
  w&height = z_km
  w2plt    = w(height|:,radius|:)

;*****RESOURCES***************************************************************************
  res                          = True
  res@cnFillOn                 = True
  res@gsnDraw                  = False
  res@gsnFrame                 = False
  res@gsnYAxisIrregular2Linear = True
  res@gsnXAxisIrregular2Linear = True
;  res@cnFillOn                = True
;  res@cnLinesOn               = False
;  res@lbLabelBarOn            = True
;  res@trXMinF                 = 0.    ;-1.0
;  res@trXMaxF                 = 500.     ;1.0
;  res@cnInfoLabelOn           = False
;  res@cnLineLabelsOn          = False
;  res@lbBoxLinesOn            = False
;;  res@cnLevelSelectionMode    = "ManualLevels"     ; set manual contour levels
;;  res@cnMinLevelValF	       = -0.5                ; set min contour level
;;  res@cnMaxLevelValF	       = 0.5  ;100                ; set max contour level
;;  res@cnLevelSpacingF	       = 0.1    ;5           ; set contour spacing
;  res@tiXAxisString           = "Tangential Velocity (m/s)"
;  res@tiYAxisString           = "Height (km)"

;*****PLOT VERTICAL VELOCITY**************************************************************
  wks = gsn_open_wks("png","test")
  gsn_define_colormap(wks,"rainbow.rgb")
  plot = gsn_csm_contour(wks,w2plt,res)

  draw(plot)
  frame(wks)