;
; plot_dns.ncl: ncl routines for regridding and plotting level 2 data
;
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;**************************************************
 begin

   source = "rvsj"
   fi = addfile("sondes_"+source+"_2.0.nc","r")

   i0 = 0
   i1 = min((/229,dimsizes(fi->DeltaP)-1/))
;
   rico_day = fi->rico_day
;
; define coordinates of uniform grid and a template variable for gridded data
;
   rday = fspan( -24.5,24.5,150)
   rday@long_name = "Days after 2005-01-01"
   rday@units = "days"

   plvl = new( (/dimsizes(fi->DeltaP)/), "float")
   plvl = (/avg(fi->SurfacePressure) -fi->DeltaP/)
   plvl@long_name = "Pressure"
   plvl@units = "Pa"

   template = new ( (/dimsizes(plvl),dimsizes(rday)/), "float")
   template@_FillValue = -999.
   template!0 = "plvl"
   template&plvl = plvl
   template!1 = "rday"
   template&rday = rday
;
; define coordinates of input data 
;
   xi   = rico_day
   yi   = avg(fi->SurfacePressure) - fi->DeltaP
   yi@long_name = "Pressure"
   xi@long_name = rday@long_name
   yi@units = "Pa"

   fldi = (/fi->tair/)
   fldi@_FillValue = -999.
   fldi!0 = "yi"
   fldi&yi = yi
   fldi!1 = "xi"
   fldi&xi = xi
;
; construct a pressure field for plvl data
;
   prsi = (/fi->tair/)
   do k=0,dimsizes(xi)-1
     prsi(:,k) = fi->SurfacePressure(k) - fi->DeltaP
   end do
;
; calculate potential temperature on original data filling missing values along
; the time axis
;
   theta = fldi
   gamma = 287.05/1004.
   theta = (/fldi*(100000./prsi)^gamma/)
;
;
; assemble all the data from a given period, and deduce the average value.
;
  tht = template
  q = template
  u = template
  v = template

  qi = fldi
  ui = fldi
  vi = fldi

  qi = (/mixhum_ptrh(prsi/100.,fi->tair,fi->rhum,-2)/)
  ui = (/-fi->wspd*sin(fi->wdir*3.14159/180.)/)
  vi = (/-fi->wspd*cos(fi->wdir*3.14159/180.)/)

  if (source .ne. "c130") then
    dt = 0.5 ;(rday(1)-rday(0))/2.
  else
    dt = (rday(1)-rday(0))/2.
  end if
  do k=0,dimsizes(rday)-1
    xt =mask(rico_day,rico_day.gt.rday(k)-dt.and.rico_day.lt.rday(k)+dt,True)
    ia =minind(xt)
    ib =maxind(xt)
    if (.not.ismissing(ia)) then
      tht(:,k) = dim_avg(theta(yi|:,xi|ia:ib))
      q(:,k)   = dim_avg(qi(yi|:,xi|ia:ib))
      u(:,k)   = dim_avg(ui(yi|:,xi|ia:ib))
      v(:,k)   = dim_avg(vi(yi|:,xi|ia:ib))
    end if
  end do
;
; create plots
;
  plots = new( (/12/), "graphic")
  wks  = gsn_open_wks ("pdf", "time_height_"+source)
  gsn_define_colormap (wks,"BlAqGrYeOrRe")

  cnres = True
  cnres@gsnFrame = False
  cnres@gsnDraw  = False
  cnres@gsnSpreadColors     = True  
  cnres@gsnSpreadColorStart = 2
  cnres@gsnSpreadColorEnd   = -2

  cnres@vpHeightF      = 0.15
  cnres@cnRasterModeOn = True
  cnres@cnFillOn       =  True  
  cnres@cnLinesOn      =  False
  cnres@cnLineLabelsOn =  False
  cnres@cnLevelSelectionMode = "ManualLevels"
  cnres@lbLabelStride  = 5

  cnres@tiYAxisFontHeightF = 0.01
  cnres@tiXAxisFontHeightF = 0.01
  fldo = tht
  fldo_avg = dim_avg(fldo)
  fldop = fldo
  do k=0,dimsizes(rday)-1
    fldop(:,k) = (/fldo(:,k) - fldo_avg/)
  end do

  cnres@tiXAxisOn = False
  cnres@cnMinLevelValF = 297.
  cnres@cnMaxLevelValF = 317.
  if (source.ne."c130") then
    cnres@cnMaxLevelValF = 367.
  end if
  cnres@cnLevelSpacingF = (cnres@cnMaxLevelValF - cnres@cnMinLevelValF)/10.
  fldo@long_name = "Potential Temperature"
  fldo@units = "K"
  plots(0) = gsn_csm_pres_hgt(wks,fldo, cnres)
  cnres@tiXAxisOn = True
  cnres@cnMaxLevelValF =  2.
  cnres@cnMinLevelValF = -cnres@cnMaxLevelValF
  cnres@cnLevelSpacingF = cnres@cnMaxLevelValF/10.
  fldop@units = "K"
  fldop@long_name = "Anomaly"
  plots(1) = gsn_csm_pres_hgt(wks,fldop(i0:i1,:), cnres)

  pnl_res = True
  pnl_res@gsnMaximize = True
  pnl_res@gsnFrame = False
  gsn_panel(wks,plots(0:1),(/2,1/),pnl_res)
  frame(wks)
;
; repeat analysis for specific humidity
;
  fldo  = q
  fldo_avg = dim_avg(fldo)
  fldop = fldo
  do k=0,dimsizes(rday)-1
    fldop(:,k) = (/fldo(:,k) - fldo_avg/)
  end do

  cnres@tiXAxisOn = False
  fldo@long_name = "Specific Humidity"
  fldo@units = "g/kg"
  cnres@cnLevelSpacingF = 1.
  cnres@cnMinLevelValF =  2.
  cnres@cnMaxLevelValF =  15.
  plots(0) = gsn_csm_pres_hgt(wks,fldo, cnres)
  cnres@tiXAxisOn = True
  cnres@cnMaxLevelValF =  3.
  cnres@cnMinLevelValF = -cnres@cnMaxLevelValF
  cnres@cnLevelSpacingF = cnres@cnMaxLevelValF/10.
  fldop@long_name = "Anomaly"
  fldop@units = "g/kg"
  plots(1) = gsn_csm_pres_hgt(wks,fldop(i0:i1,:), cnres)

  gsn_panel(wks,plots(0:1),(/2,1/),pnl_res)
  frame(wks)
;
; zonal wind 
;
  fldo  =u
  fldo_avg = dim_avg(fldo)
  fldop = fldo
  do k=0,dimsizes(rday)-1
    fldop(:,k) = (/fldo(:,k) - fldo_avg/)
  end do

  cnres@tiXAxisOn = False
  fldo@long_name = "Zonal Wind"
  fldo@units = "m/s"
  cnres@cnLevelSpacingF = 3.
  cnres@cnMinLevelValF =  -15.
  if (source.ne."c130") then
    cnres@cnMaxLevelValF = 45.
  end if
  plots(0) = gsn_csm_pres_hgt(wks,fldo, cnres)
  cnres@tiXAxisOn = True
  cnres@cnMaxLevelValF =  10.
  cnres@cnMinLevelValF = -cnres@cnMaxLevelValF
  cnres@cnLevelSpacingF = cnres@cnMaxLevelValF/10.
  fldop@long_name = "Anomaly"
  fldop@units = "m/s"
  plots(1) = gsn_csm_pres_hgt(wks,fldop(i0:i1,:), cnres)

  gsn_panel(wks,plots(0:1),(/2,1/),pnl_res)
  frame(wks)
;
; meridional wind 
;
  fldo  =v
  fldo_avg = dim_avg(fldo)
  fldop = fldo
  do k=0,dimsizes(rday)-1
    fldop(:,k) = (/fldo(:,k) - fldo_avg/)
  end do

  cnres@tiXAxisOn = False
  fldo@long_name = "Meridional Wind"
  fldo@units = "m/s"
  cnres@cnMaxLevelValF =  15.
  if (source.ne."c130") then
    cnres@cnMaxLevelValF = 30.
  end if
  cnres@cnMinLevelValF = -cnres@cnMaxLevelValF
  cnres@cnLevelSpacingF = cnres@cnMaxLevelValF/10.
  plots(0) = gsn_csm_pres_hgt(wks,fldo, cnres)
  cnres@tiXAxisOn = True
  cnres@cnMaxLevelValF =  7.5
  cnres@cnMinLevelValF = -cnres@cnMaxLevelValF
  cnres@cnLevelSpacingF = cnres@cnMaxLevelValF/10.
  fldop@long_name = "Anomaly"
  fldop@units = "m/s"
  plots(1) = gsn_csm_pres_hgt(wks,fldop(i0:i1,:), cnres)

  gsn_panel(wks,plots(0:1),(/2,1/),pnl_res)
  frame(wks)


 end