;**************************************************
; HAZWX SOUNDINGS!!!! 
;
; Concepts illustrated:
;   - Reading GRIB data
;   - Using getind_latlon2d to determine grid locations
;   - Drawing Skew-T plots at nearest grid locations  
;                    to user specified locations
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/home/Soundings/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "/home/Soundings/names_grib2.txt"
load "/home/Soundings/plot_hodo.ncl"
external CONV "/home/Soundings/Extra/test.so"
;***********************************************
; The GRIB file's contents can be examined via:
;     ncl_filedump -itime hrr.grb2  | less
;***********************************************
;  MODEL NAME MUST BE ONE OF THE FOLLOWING:
;  "Rapid_Refresh"
;  "NAM_12KM"
;  "NAM_4KM"
;  "HRRR"
;  "Hi_Res_GFS"
;  "GFS_13KM"
;************************************************
;Arguements coming in!

  a0 = getenv("NCL_ARG_2")  ; nc file

  a1 = getenv("NCL_ARG_3")  ; model name

  a2 = getenv("NCL_ARG_4")  ; valid time

  a5 = getenv("NCL_ARG_5")  ; forecast hour

  a3 = getenv("NCL_ARG_6")  ; lat

  a4 = getenv("NCL_ARG_7")  ; lon

; optional argument: parcel_type
; parcel_type could define the type of parcel to be used
; examples: surface-based, mixed-layer, most-unstable,
; custom
; a5 = getenv("NCL_ARG_8")

;Load this for extra cool sounding stuff

;GAME FACE!
print("Making sounding for " + a1 + " model...")
; --- Read GRIB file------------;
  diri  = ""
  fili  = a0
 
  ; force a 'time' dimension
  setfileoption("nc","Format","LargeFile")

  f     = addfile(diri+fili,"r") 
                                     ; Model grid point locations
  if (isStrSubset(a1,"GFS")) then
   lat1d = f->lat_0
   lon1d = f->lon_0
   lat2d = new((/dimsizes(lat1d),dimsizes(lon1d)/),float)
   lon2d = new((/dimsizes(lat1d),dimsizes(lon1d)/),float)
   do a = 0,dimsizes(lat1d)-1
    lon2d(a,:) = lon1d
   end do
   do b = 0,dimsizes(lon1d)-1
    lat2d(:,b) = lat1d
   end do
  else
   lat2d = f->gridlat_0               ; ( ygrid_0, xgrid_0 ) 
   lon2d = f->gridlon_0
  end if

; initial diagnostics
;  print("lat2d: min="+min(lat2d)+"  ;  max="+max(lat2d))
;  print("lon2d: min="+min(lon2d)+"  ;  max="+max(lon2d))

; --- Specify one or more locations
 
  lat   = (/ stringtofloat(a3) /)
  lon   = (/ stringtofloat(a4) /)

;*************************
; create plot(s)
;*************************

  skewtOpts = True
  skewtOpts@DrawIsotherm = True
  skewtOpts@DrawIsobar = True
  skewtOpts@DrawMixRatio = True
  skewtOpts@DrawDryAdiabat = True
  skewtOpts@DrawMoistAdiabat = True
  skewtOpts@DrawWind = True
  skewtOpts@DrawStandardAtm = False
  skewtOpts@DrawColLine = True 
  skewtOpts@DrawColAreaFill = False
  skewtOpts@DrawFahrenheit = True
  skewtOpts@DrawHeightScale = False
  skewtOpts@DrawHeightScaleFt = True 
 
  localOpts = True
  localOpts@Cape      = True 
  localOpts@ThermoInfo= True
  ;localOpts@Parcel    = 0   
  gsOpts = True
 
  dataOpts = True
  dataOpts@colTemperature="red"
  dataOpts@colDewPt = "green"
  dataOpts@colCape = "orange"
  if (isStrSubset(a1,"Hi_Res_GFS").or.isStrSubset(a1,"GFS_13KM")) then
   dataOpts@Wthin     = 1        ; plot every n-th wind barb
  else
   dataOpts@Wthin     = 2        ; plot every n-th wind barb
  end if
    
  ; find grid point nearest the user specified location
  nm  = getind_latlon2d (lat2d,lon2d, lat(0), lon(0))
  nn  = nm(0,0)
  mm  = nm(0,1)
  print("location=("+lat(0)+","+lon(0)+")  grid=("+lat2d(nn,mm)+","+lon2d(nn,mm)+")")

 ; Load in our variables
 ; ALL models pull in the following  
  p_pa     = f->lv_ISBL0         ; [Pa]
  p_mb = 0.01*p_pa ; convert to [hPa]
  p_mb@units = "hPa"
 ; Now to account for the differences
  if (isStrSubset(a1,"Hi_Res_GFS").or.isStrSubset(a1,"GFS_13KM")) then               
   z = f->$GFSGHNAME$(:,nn,mm)
   u = f->$GFSUNAME$(:,nn,mm)
   v = f->$GFSVNAME$(:,nn,mm)
   u10 = f->$GFSU10MNAME$(0,nn,mm)
   v10 = f->$GFSV10MNAME$(0,nn,mm)
  ; t2m = f->$GFST2MNAME$(nn,mm)
  ; td2m = f->GFSDPT2MNAME$(nn,mm)
   hsfc = f->$GFSGHSFCNAME$(nn,mm)
   psfc = f->$GFSPSFCNAME$(nn,mm)
   cape = f->$GFSCAPENAME$(nn,mm)  
   cin = f->$GFSCINNAME$(nn,mm)        
   li = f->$GFSLINAME$(nn,mm)        
   ustm = f->$GFSUSTMNAME$(nn,mm)
   vstm = f->$GFSVSTMNAME$(nn,mm)    
   srh03 = f->$GFSHLCYAGLNAME$(nn,mm)
   pw = f->$GFSPWNAME$(nn,mm)        
  else ; not-GFS names
   z = f->$GHNAME$(:,nn,mm) 
   u = f->$UNAME$(:,nn,mm)  
   v = f->$VNAME$(:,nn,mm)   
   if (isStrSubset(a1,"HRRR")) then    
    u10 = f->$U10MNAME$(nn,mm)
    v10 = f->$V10MNAME$(nn,mm)
   else
    u10 = f->$U10MNAME$(0,nn,mm)
    v10 = f->$V10MNAME$(0,nn,mm)
   end if
 ;  t2m = f->$T2MNAME$(nn,mm)
  ; td2m = f->$DPT2MNAME$(nn,mm)
   hsfc = f->$GHSFCNAME$(nn,mm)
   psfc = f->$PSFCNAME$(nn,mm)
   cape = f->$CAPENAME$(nn,mm)
   cin = f->$CINNAME$(nn,mm)
   li = f->$LINAME$(nn,mm)
;   if (isStrSubset(a1,"NAM_4KM") .or. isStrSubset(a1,"HRRR")) then
;    lclz = f->$LCLZNAME$(nn,mm)
;   end if
   pw = f->$PWNAME$(nn,mm)
   if (isStrSubset(a1,"HRRR")) then
    ustm = f->$USTMNAME$(nn,mm)
    vstm = f->$VSTMNAME$(nn,mm)
   end if
   if (isStrSubset(a1,"HRRR")) then
;    ushr01 = f->$USHRNAME$(0,nn,mm)
;    vshr01 = f->$VSHRNAME$(0,nn,mm)
;    ushr06 = f->$USHRNAME$(1,nn,mm)
;    vshr06 = f->$VSHRNAME$(1,nn,mm)
   else if (isStrSubset(a1,"RAP")) then
    ushr06 = f->$USHRNAME$(nn,mm)
    vshr06 = f->$VSHRNAME$(nn,mm)
   end if
   end if
   if (isStrSubset(a1,"NAM_12KM")) then
    ; GFS already taken care of, so no need to add it to this IF statement
    srh03 = f->$HLCYAGLNAME$(nn,mm) ; NCL switch
   end if
  end if

  ; We need dewpoint in C
  ; Also need to cut off a level in GFS data
  if (isStrSubset(a1,"NAM_12KM").or.isStrSubset(a1,"NAM_4KM") .or. isStrSubset(a1,"Rapid_Refresh")) then
   tk1 = f->$TNAME$(:,nn,mm)
   rh = f->$RHNAME$(:,nn,mm)
   tdk1 = dewtemp_trh(tk1,rh)
  else if (isStrSubset(a1,"HRRR")) then
   tk1 = f->$TNAME$(:,nn,mm)       
   rh = f->$RHNAME$(:,nn,mm)     
   tdk1 = f->$DPTNAME$(:,nn,mm)
  else if (isStrSubset(a1,"Hi_Res_GFS").or.isStrSubset(a1,"GFS_13KM")) then
   ; In the GFS, temperature exists at lv_ISBL0 levels, whereas 
   ; RH only exists at lv_ISBL6 levels, which is the same as lv_ISBL0
   ; except not 20 mb. So...have to remove the element that corresponds to 
   ; the 20 mb level in ALL arrays
   tk1 = f->$GFSTNAME$(:,nn,mm)
   rh = f->$GFSRHNAME$(:,nn,mm)
   nz = dimsizes(tk1)
   temp_tk1 = new(dimsizes(tk1)-1,float)
   temp_p_pa = new(dimsizes(tk1)-1,float)
   temp_p_mb = new(dimsizes(tk1)-1,float)
   temp_u = new(dimsizes(tk1)-1,float)
   temp_v = new(dimsizes(tk1)-1,float)
   temp_z = new(dimsizes(tk1)-1,float)
   count = 0
   do k = 0,nz-1
    if (p_pa(k) .ne. 2000) then
     temp_tk1(count) = tk1(k)
     temp_p_pa(count) = p_pa(k)
     temp_p_mb(count) = p_mb(k)
     temp_v(count) = v(k)
     temp_u(count) = u(k)
     temp_z(count) = z(k)
     count = count + 1
    end if
   end do
   delete([/tk1,p_pa,p_mb,u,v,z/])
   tk1 = temp_tk1
   p_pa = temp_p_pa
   p_mb = temp_p_mb
   u = temp_u
   v = temp_v
   z = temp_z
   delete([/temp_tk1,temp_p_pa,temp_p_mb,temp_u,temp_v,temp_z/])
   tdk1 = dewtemp_trh(tk1,rh)
  end if
  end if
  end if
  qv = mixhum_ptd(p_pa,tdk1,1) ; obtain water vapor mixing ratio

;  parcel_type = stringtofloat(a5)
;  if (parcel_type == 2) then ; mixed-layer
;   ; assume mixing over lowest 1000 m
;   p_sum = 0.
;   tk_sum = 0.
;   qv_sum = 0.
;   nz = dimsizes(z)
;   do k = nz(0)-1,0,-1
;    index = nz(0)-1-k
;    if (z(k) - hsfc < 1000.0) then
;     p_sum = p_sum + p_mb(k)
;     tk_sum = tk_sum + p_mb(k)*tk1(k)
;     qv_sum = qv_sum + p_mb(k)*qv(k)
;    else
;     break
;    end if
;   end do
;   ml_tk = tk_sum/psum
;   ml_qv = qv_sum/psum
;  end if



; Convert variables to appropriate units (celcius for temps, mb for pressure, kts for wind)
  r2d = 45.0/atan(1.0)
  tc = tk1 - 273.15
  tdc = tdk1 - 273.15
  if (isStrSubset(a1,"NAM_12KM") .or. isStrSubset(a1,"NAM_4KM")) then
   li = li - 273.15
  end if
  li = toint(li)
  plcl = lclvl(p_mb,tk1,tdk1)
  plcl = plcl(::-1)
  plcls = toint(plcl)
  cape = toint(cape)
  cin = tointeger(cin)
  Uspd1 = 1.94384449*sqrt(u*u + v*v)
  Udir1 = atan2(u, v) * r2d + 180.0
  pw = tofloat(pw*0.0393701)
  pws = sprintf("%4.2f",pw)

  ;invert column data so it is going from surface to TOA
  p_mb = p_mb(::-1)
  u = u(::-1)      
  v = v(::-1) 
  z = z(::-1)
  tc = tc(::-1)
  tdc = tdc(::-1)
  Uspd1 = Uspd1(::-1)
  Udir1 = Udir1(::-1)

 ;If shear isn't in the file, need to compute it
  nz = dimsizes(z)
  if (.not. isStrSubset(a1,"RAP")) then
   ;RAP already has 6 km shear
   ;find wind at 6 km AGL
   do k = 0,nz(0)-1
    ; going up from surface
    if (z(k) - hsfc .gt. 6000.0) then
     ; at the first level up from the surface that we're > 6 km AGL, stop
     ; linearly interpolate between this and the previous level (down)
     ; to get 6 km wind components
     u6000 = (u(k)-u(k-1))/(z(k)-z(k-1))*(6000.0+hsfc-z(k)) + u(k)
     v6000 = (v(k)-v(k-1))/(z(k)-z(k-1))*(6000.0+hsfc-z(k)) + v(k)
     break
    end if
   end do
   end if ; RAP has 6 km shear, but not 1 km shear
   ; find wind at 1 km AGL
   do k = 0,nz(0)-1 
    ; the order of the looping means we're going up from the bottom of the model output
    if (z(k) - hsfc .gt. 1000.0) then
     ; linearly interpolate between this and the previous level (below)
     ; to get 1 km wind components 
     u1000 = (u(k)-u(k-1))/(z(k)-z(k-1))*(1000.0+hsfc-z(k)) + u(k)
     v1000 = (v(k)-v(k-1))/(z(k)-z(k-1))*(1000.0+hsfc-z(k)) + v(k)
     break
    end if
   end do 
   ushr06 = u6000-u10
   vshr06 = v6000-v10
   ushr01 = u1000-u10              
   vshr01 = v1000-v10
   ; shear magnitude is calculated and converted later

  ;Determine storm motion vector using ID method from Bunkers et al. (2000), WAF
  if (.not. isStrSubset(a1,"Hi_Res_GFS")) then
   ; Compute 0-6 km non-pressure weighted mean wind
   nz = dimsizes(z)
   count05 = 0
   count55 = 0
   count6 = 0
   u05_sum = 0.0
   v05_sum = 0.0
   u55_sum = 0.0
   v55_sum = 0.0
   u6_sum = 0.0
   v6_sum = 0.0
   do k = 0,nz(0)-1
    if (z(k) - hsfc .le. 6000.0 .and. z(k) .ge. hsfc) then
     u6_sum = u6_sum + u(k)
     v6_sum = v6_sum + v(k)
     count6 = count6 + 1
     if (z(k) - hsfc .le. 500.0 .and. z(k) .ge. hsfc) then
      u05_sum = u05_sum + u(k)
      v05_sum = v05_sum + v(k)
      count05 = count05 + 1
     end if
     if ((z(k) - hsfc .ge. 5500.0) .and. (z(k) - hsfc .le. 6000.0)) then
      u55_sum = u55_sum + u(k)
      v55_sum = v55_sum + v(k)
      count55 = count55 + 1
     end if
    else if (z(k) - hsfc .gt. 6000.0)
     break
    end if
    end if
   end do
   u6_mean = u6_sum/count6
   v6_mean = v6_sum/count6
   ; make sure we don't get INFs for the narrow mean winds that have no values in a 500 m thick layer
   if (count05 .eq. 0) then
    u05_mean = u10
    v05_mean = v10
   else
    u05_mean = u05_sum/count05
    v05_mean = v05_sum/count05
   end if
   if (count55 .eq. 0) then
    ; I hope this never happens because to account for it, i'll have to allow u6000 and v6000 to be computed
    ; for every model type, and that consumes time
    u55_mean = u6000
    v55_mean = v6000
   else
    u55_mean = u55_sum/count55
    v55_mean = v55_sum/count55
   end if
   ; next, find the shear between the mean 0-0.5 km wind and the 5.5-6 km wnid
   ushear05 = u55_mean - u05_mean
   vshear05 = v55_mean - v05_mean
   ; Now determine storm motion by moving 7.5 m/s away from the shear vector at the 
   ; location of the mean wind
   CONST = sqrt(ushear05*ushear05+vshear05*vshear05)
   uadd = 7.5*vshear05 / CONST
   vadd = -7.5*ushear05 / CONST
   ustm = u6_mean + uadd
   vstm = v6_mean + vadd
  end if

  ; Compute storm relative helicity for models that don't already have it
  nz = dimsizes(u)
  curve = new((/50/),float,0.0)
  DZ = new((/50/),float,0.0)
  i = 0
  do k = 0,nz(0)-1
   if (z(k) .ge. hsfc .and. z(k) .lt. 3000.0+hsfc) then
     DZ(k) = z(k+1)-z(k)
     dudz = (u(k+1)-u(k)) / DZ(k)
     dvdz= (v(k+1)-v(k)) / DZ(k)
     u_sr = u(k) - ustm
     v_sr = v(k) - vstm
     curve(i) = -u_sr*dvdz + v_sr*dudz
     i = i + 1
   else if (z(k) .gt. 3000.0+hsfc) then
     DZ(k) = z(k+1)-z(k)
     dudz = (u(k+1)-u(k)) / DZ(k)
     dvdz= (v(k+1)-v(k)) / DZ(k)
     u_sr = u(k) - ustm            
     v_sr = v(k) - vstm       
     curve(i) = -u_sr*dvdz + v_sr*dudz
    break
   end if
   end if
  end do
  SRH_sum = 0.0
  i = 0
  do k = 0,nz(0)-1
   if (z(k) .ge. hsfc .and. z(k) .lt. 1000.0+hsfc) then
    SRH_sum = SRH_sum + 0.5*DZ(k)*(curve(i)+curve(i+1))
    i = i + 1
   else if (z(k) .gt. 1000.0+hsfc) then
    srh01 = SRH_sum
    break
   end if
   end if
  end do
  delete(SRH_sum)
  SRH_sum = 0.0
  i = 0
  do k = 0,nz(0)-1
   if (z(k) .ge. hsfc .and. z(k) .lt. 3000.0+hsfc) then
    SRH_sum = SRH_sum + 0.5*DZ(k)*(curve(i)+curve(i+1))
    i = i + 1
   else if (z(k) .gt. 3000.0+hsfc) then
    srh03 = SRH_sum
    break
   end if
   end if
  end do
  srh01 = toint(srh01)
  srh03 = toint(srh03)

  ;Truncate data to levels above the surface
  nz = dimsizes(z)
  count = 0
  do k = 0,nz-1
   if (z(k) .ge. hsfc) then
    count = count + 1
   end if
  end do
  new_tc = new(count,float)
  new_tdc = new(count,float)
  new_z = new(count,float)
  new_spd = new(count,float)
  new_dir = new(count,float)
  new_p_mb = new(count,float)
  do k = 0,count
   if (z(k) .ge. hsfc) then
    new_tc(k) = tc(k)
    new_tdc(k) = tdc(k)
    new_z(k) = z(k)-hsfc
    new_spd(k) = Uspd1(k)
    new_dir(k) = Udir1(k)
    new_p_mb(k) = p_mb(k)
   end if
  end do


ctemp = 0.0
CONV:: sndg_anal(nlevs,t,p,r,u,v,ctemp)


  ;Formatting - AFTER shear and storm motion are calculated
  SMspd = tointeger(1.94384449*sqrt(ustm*ustm + vstm*vstm))          
  SMdir = tointeger(atan2(ustm, vstm) * r2d + 180.0)
  SMdir_string = sprinti("%0.3i",SMdir)
  shear01 = tointeger(1.94384449*sqrt(ushr01*ushr01 + vshr01*vshr01))
  shear06 = tointeger(1.94384449*sqrt(ushr06*ushr06 + vshr06*vshr06))

  itime = a2
  uscore = "_"
  space = " "

  skewtOpts@tiMainString  = str_sub_str(a1,uscore,space)+" | Valid "+str_sub_str(a2,uscore,space)+" | ("+lat(0)+","+lon(0)+")"

  ;Lets setup all these wind barb settings first
  wmsetp("wbt", 0.40) ; 0.33 is the default 
  wmsetp("blw", 2.0) ; 2.0 is our default 

  ;set up our image output
  wks_type = "png"
  wks_type@wkWidth = 685
  wks_type@wkHeight = 650

  ;each location will have a different file name  
    
  wks  = gsn_open_wks (wks_type, "skewt_"+a1+"_"+a2+"_f"+a5+"_"+lat+"_"+lon)
  skewt_bkgd = skewT_BackGround (wks, skewtOpts)
  skewt_data = skewT_PlotData(wks, skewt_bkgd, new_p_mb, new_tc, new_tdc, new_z, new_spd, new_dir, dataOpts)
  draw(skewt_data)
  draw(skewt_bkgd)

; Draw the inset info box background
;
bpRes                  = True
bpRes@tmXBOn           = False
bpRes@tmXBBorderOn     = False
bpRes@tmXTOn           = False
bpRes@tmXTBorderOn     = False
bpRes@tmYROn           = False
bpRes@tmYRBorderOn     = False
bpRes@tmYLOn           = False
bpRes@tmYLBorderOn     = False
bpRes@gsnFrame         = False
infobox = gsn_blank_plot(wks,bpRes)
setvalues infobox
  "vpXF" : 0.567
  "vpYF" : 0.9245
  "vpWidthF" : 0.195
  "vpHeightF" : 0.284
end setvalues
gsRes                  = True
gsRes@gsFillColor      = "Gray"
gsn_polygon(wks,infobox,(/0.0,1.0,1.0,0.0/),(/1.0,1.0,0.35,0.35/),gsRes)
;
; Draw the inset info box text
;
txidcape = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "Sounding Information"
  "txFontHeightF" : 0.017
  "txFontThicknessF" : 2.0
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidcape)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    -0.06
  "amOrthogonalPosF":  -0.08
end setvalues
txidcape = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "CAPE: "+cape+" J/kg"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidcape)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.01
end setvalues
;
txidcin = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "CIN: "+cin+" J/kg"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
ianno = NhlAddAnnotation(infobox,txidcin)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.07
end setvalues
;
;no mu values for now
;txidmucape = create "TextItems" textItemClass wks
  ;"txFuncCode" : "~"
  ;"txString" : "MUCAPE: "+mucape+" J/kg"
  ;"txFontHeightF" : 0.012
  ;"txFontColor" : "black"
;end create
;ianno = NhlAddAnnotation(infobox,txidmucape)
;setvalues ianno
  ;"amTrackData":      False
  ;"amZone":           1
  ;"amSide":           "top"
  ;"amJust":           "topleft"
  ;"amParallelPosF":    0.05
  ;"amOrthogonalPosF":  0.13
;end setvalues
;
;txidmucin = create "TextItems" textItemClass wks
  ;"txFuncCode" : "~"
  ;"txString" : "MUCIN: "+mucin+" J/kg"
  ;"txFontHeightF" : 0.012
  ;"txFontColor" : "black"
;end create
;ianno = NhlAddAnnotation(infobox,txidmucin)
;setvalues ianno
  ;"amTrackData":      False
  ;"amZone":           1
  ;"amSide":           "top"
  ;"amJust":           "topleft"
  ;"amParallelPosF":    0.05
  ;"amOrthogonalPosF":  0.19
;end setvalues
;
txidli = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "LI: "+li+" C"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidli)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.13
end setvalues
txidlcl = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "LCL: "+plcls+" mb"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidlcl)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.19
end setvalues
;
if (False) then
txidlpl = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "LPL: "+lpl+" mb"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidlpl)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.19
end setvalues
end if
;
txidconvt = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "Conv T: "+ctemp+" F"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidconvt)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.25
end setvalues
;
txidsrh03 = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "0-3km SRH: "+srh03+" m~S~2~N~/s~S~2~E~ "
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidsrh03)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.31
end setvalues
;
txidsrh01 = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "0-1km SRH: "+srh01+" m~S~2~N~/s~S~2~E~ "
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidsrh01)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.37
end setvalues
;
txidshear06 = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "0-6km Shear: "+shear06+" kt"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidshear06)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.43
end setvalues
;
txidshear01 = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "0-1km Shear: "+shear01+" kt"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidshear01)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.49
end setvalues
;
txidcell = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "CELL: "+SMdir_string+"@"+SMspd+" kt"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidcell)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.54
end setvalues
;
txidpw = create "TextItems" textItemClass wks
  "txFuncCode" : "~"
  "txString" : "PW: "+pws+" in"
  "txFontHeightF" : 0.012
  "txFontColor" : "black"
end create
;
ianno = NhlAddAnnotation(infobox,txidpw)
setvalues ianno
  "amTrackData":      False
  "amZone":           1
  "amSide":           "top"
  "amJust":           "topleft"
  "amParallelPosF":    0.05
  "amOrthogonalPosF":  0.60
end setvalues
;
draw(infobox)
;
; Draw the hodograph on the sounding
;
hodores = True
hodores@gsnDraw  = False
hodores@gsnFrame = False
hodo1 = plot_hodo(wks, u, v, z, ustm,vstm, hodores)
;
setvalues hodo1
  "vpXF" : 0.07
  "vpYF" : 0.925
  "vpWidthF" : 0.27
  "vpHeightF" : 0.27
end setvalues
;
draw(hodo1)
;
frame(wks)