;
; plot profiles averaged over all of level 2 sonde data
; 	theta_e, theta_es, u and v
;
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "~/Library/NCL/thermo.ncl"
;**************************************************
 begin

   FlightDays = False
   ShipDays = True
   if (FlightDays) then
     wks  = gsn_open_wks ("pdf", "profiles_fdays")
   else
     if (ShipDays) then
       wks  = gsn_open_wks ("pdf", "profiles_sdays")
     else
       wks  = gsn_open_wks ("pdf", "profiles_adays")
     end if
   end if

   spnt = addfile("sondes_spnt_3.0.nc","r")
   c130 = addfile("sondes_rvsj_3.0.nc","r")

   rico_day = spnt->rico_day
   nsondes = dimsizes(rico_day)
;
; construct a pressure field for plvl data
;
   prsi = spnt->tair
   prsi@long_name = "pressure"
   prsi@units = "Pa"
   do k=0,nsondes-1
     prsi(:,k) = spnt->SurfacePressure(k) - spnt->DeltaP
   end do
   prsi = (/prsi/100./)

   tair = spnt->tair
   rhum = spnt->rhum
   uvel = -spnt->wspd * sin(spnt->wdir * 3.14159/180.)
   vvel = -spnt->wspd * cos(spnt->wdir * 3.14159/180.)
;
; only average data from flight days
;
   if (FlightDays) then
     do n = 0, nsondes-1
       if (ismissing(c130->SurfacePressure(n))) then
         tair(:,n) = tair@_FillValue
         rhum(:,n) = rhum@_FillValue
         uvel(:,n) = uvel@_FillValue
         vvel(:,n) = vvel@_FillValue
       end if
     end do
   end if
;
; only average data from ship days
;
   if (ShipDays) then
     do n = 0, nsondes-1
       if (ismissing(c130->SurfacePressure(n))) then
         tair(:,n) = tair@_FillValue
         rhum(:,n) = rhum@_FillValue
         uvel(:,n) = uvel@_FillValue
         vvel(:,n) = vvel@_FillValue
       end if
     end do
   end if
;
; get moisture and derive theta_e
;
   theta_e = spnt->tair
   theta_e@long_name = "Equivalent Potential Temperature"
   rv      = spnt->tair
   rv@long_name = "Water Vapor Mixing Ratio"
   rv@units = "g/kg"

   rv      = (/ mixhum_ptrh(prsi,tair,rhum,1) * 1000./)
   rl      = rv
   rs      = rv
   rl      = 0.
   rs      = (/ get_rs(tair,prsi*100.)*1000. /)
   theta_e  = (/ get_theta_e(tair,prsi,rv,rl)/)
   theta_es = (/ get_theta_e(tair,prsi,rs,rl)/)
;
; construct y axis
;
   yi   = (/ -avg(spnt->SurfacePressure) + spnt->DeltaP/)/100.
   yi@long_name = "Pressure"
   yi@units = "Pa"
;
; create x-variables from theta_e and theta_es
;
   xi = yi
   xi!0 = "pressure" 
   xi&pressure = yi

   x1  = xi
   x2  = xi
   x1  = (/dim_avg(theta_e)/)
   x2  = (/dim_avg(theta_es)/)

   delete(rv)
   delete(rl)
   delete(rs)
   delete(tair)
   delete(rhum)
   delete(prsi)
;
; repeat analysis for dropsondes
;
;
; construct a pressure field for plvl data
;
   prsi = c130->tair
   do k=0,dimsizes(c130->rico_day)-1
     prsi(:,k) = c130->SurfacePressure(k) - c130->DeltaP
   end do
   prsi = (/prsi/100./)

   tair = c130->tair
   rhum = c130->rhum
   uveld = -c130->wspd * sin(c130->wdir * 3.14159/180.)
   vveld = -c130->wspd * cos(c130->wdir * 3.14159/180.)
;
; get rv
;
   z1 = c130->tair
   rv = c130->tair

   rv      = (/ mixhum_ptrh(prsi,tair,rhum,1) * 1000./)
   rl      = rv
   rs      = rv
   rl      = 0.
   rs      = (/ get_rs(tair,prsi*100.)*1000. /)
   z1 = (/ get_theta_e(tair,prsi,rv,rl)/)
   z2 = (/ get_theta_e(tair,prsi,rs,rl)/)

   z3   = (/ -avg(c130->SurfacePressure) + c130->DeltaP/)/100.
   z4  = (/dim_avg(z1)/)
   z5  = (/dim_avg(z2)/)

;
; set up plot, and plot resources
;
   xyres= True
   xyres@gsnDraw = False
   xyres@gsnFrame = False
   xyres@tiYAxisString = "[hPa]"
   xyres@tiXAxisString = "[K]"
   xyres@xyLineColors = (/"dodgerblue","red"/)
   xyres@tmBorderThicknessF = 2.
   xyres@xyLineThicknessF = 2.

   if (FlightDays) then
     xyres@tiMainString = "Flight Days, :F8:Q:BF0:e:NF8:, Q:BF0:es:N"
   else
     if (ShipDays) then
       xyres@tiMainString = "Ship Days, :F8:Q:BF0:e:NF8:, Q:BF0:es:N"
     else
       xyres@tiMainString = "All Days, :F8:Q:BF0:e:NF8:, Q:BF0:es:N"
     end if
   end if

   xyres@tmYRBorderOn = False
   xyres@tmYROn = False

   xyres@tmXTBorderOn = False
   xyres@tmXTOn = False
 
   x3 = x2
   x3 = (/x1(1) - x3/) 
   x3 = mask(yi,x3.gt.0.,True)
   x3@FillValue = -999.

   xx  = mask(x2, yi.gt.-800. .and. yi .lt.-400., True)
   xbar  = round(avg(xx)*10.,1)/10.
   xmin  = round(min(x1)*10.,1)/10.
   xmax  = round(max(x2(0:10))*10.,1)/10.
   ybase = round(min(yi)*10.,1)/10
   ythsm = round(yi(maxind(x2(100:300))+100)*10.,1)/10.
   ythem = round(yi(minind(x1))*10.,1)/10.
   ylnb  = round(max(x3)*10.,1)/10.

   xyres@tmXBMode = "Explicit"
   xyres@tmXBMajorOutwardLengthF = 0.0
   xyres@tmXBMinorOutwardLengthF = 0.0
   xyres@tmXBLabelFontHeightF = .015
   xyres@tmXBLabelFont = "times-roman"
   xyres@trXMaxF = 355.
   xyres@trXMinF = 320.
   xyres@tmXBValues = (/xmin,xbar,xmax/)
   xyres@tmXBLabels = xyres@tmXBValues
   xyres@tmXBMinorValues = fspan(320.,355.,26)

   xyres@tmYLMode = "Explicit"
   xyres@tmYLMajorOutwardLengthF = 0.0
   xyres@tmYLMinorOutwardLengthF = 0.0
   xyres@tmYLLabelFontHeightF = .015
   xyres@tmYLLabelFont = "times-roman"
   xyres@trYMaxF = -200.
   xyres@trYMinF = -1020.
   xyres@tmYLValues = (/ybase,ythsm,ythem,ylnb/)
   xyres@tmYLLabels = - xyres@tmYLValues

   xyres@tiXAxisFontHeightF = 0.015
   plots = new ( (/10/), "graphic")
   xyres@gsnXRefLine = avg(x1(0:5))
   xyres@gsnXRefLineDashPattern = 0
   xyres@gsnXRefLineThicknessF = 1.0
   plots(0) = gsn_csm_xy(wks, (/x1,x2/), yi, xyres)

   xyres@xyLineThicknessF = 1.
;   xyres@xyDashPatterns = (/2,2/)
   z4  = mask(z4,z3.lt.-50., True)
   z5  = mask(z5,z3.lt.-50., True)
   ovrly = gsn_csm_xy(wks, (/z4,z5/), z3, xyres)
   overlay(plots(0),ovrly)

   if (FlightDays) then
     xyres@tiMainString = "Flight Days, :F10:u,v"
   else
     xyres@tiMainString = "All Days, :F10:u,v"
   end if

   xyres@tiXAxisString = "[ms:S:-1:N:]"
   xyres@xyLineColors = (/"black","black"/)

   x1  = (/dim_avg(uvel)/)
   x2  = (/dim_avg(vvel)/)
   z4  = (/dim_avg(uveld)/)
   z5  = (/dim_avg(vveld)/)

   z4  = mask(z4,z3.lt.-50., True)
   z5  = mask(z5,z3.lt.-50., True)

   xyres@trXMaxF =  30.
   xyres@trXMinF = -10.
   xmin = round(min(x1)*10.,1)/10.
   xmax = round(max(x1)*10.,1)/10.
   xyres@tmXBValues = (/xmin,0.,xmax/)
   xyres@tmXBLabels = xyres@tmXBValues
   delete(xyres@tmXBMinorValues)
   xyres@tmXBMinorValues = fspan(-10.,30.,41)
   xyres@xyLineThicknessF = 2.
   xyres@xyDashPatterns = (/0,1/)
   xyres@gsnXRefLine = 0.
   xyres@gsnXRefLineDashPattern = 0
   xyres@gsnXRefLineThicknessF = 1.0
   plots(1) = gsn_csm_xy(wks, (/x1,x2/), yi, xyres)

   xyres@xyLineThicknessF = 1.
;   xyres@xyDashPatterns = (/2,2/)
   ovrly = gsn_csm_xy(wks, (/z4,z5/), z3, xyres)
   overlay(plots(1),ovrly)

   pnl_res = True
   pnl_res@gsnMaximize = True
   pnl_res@gsnFrame = False
   gsn_panel(wks,plots(0:1),(/1,2/),pnl_res)
   frame(wks)

 end