;Written by Xiaoming Hu,CAPS, OU, Oct. 3, 2011 
; contact yuanfangcan@gmail.com
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"


begin

     S54511_Lat = 39.8 
     S54511_Lon = 116.4667

  latitude_3S54511  = (/S54511_Lat, S54511_Lat, S54511_Lat/)
  longitude_3S54511 = (/S54511_Lon, S54511_Lon, S54511_Lon/)
  names_3S54511     = (/"S54511",  "S54511",   "S54511"/)

 lat = S54511_Lat 
 lon = S54511_Lon 

 timewithsounding = (/5, 11, 23, 29/)
; soundingfiles = systemfunc("ls Z_UPAR_I_54342_201?????05*.SECOND")
 soundingfiles = systemfunc("ls Z_UPAR_I_54342_201?????11*.SECOND") ; also turn off adding the PBL line

  resp                  = True                      ; polyline mods desired
  resp@gsLineColor      = "red"                     ; color of lines
  resp@gsLineThicknessF = 4.0                       ; thickness of lines
  resp@gsLineDashPattern= 2                       ; thickness of lines
  tres       =  True
  tres@txFontHeightF = 0.01322


do ihours =  0,  dimsizes(soundingfiles)-1

  figurename =  "TS_profile_Theta_RH_11z2km_"+ihours
;  figurename =  "TS_profile_Theta_RH_5z_"+ihours
;  figurename =  "TS_profile_Theta_withPBL_q_"+ihours


    a4_height = 29.7 ; in centimeters, if my 
    a4_width = 23.0 ; reference is correct 
    cm_per_inch = 2.54 
 
  res1                       = True            ; plot mods desired
  res1@gsnFrame               = False                     ; don't draw yet

;  res1@mpShapeMode  = "FreeAspect"
  res1@vpHeightF             =0.42
  res1@vpWidthF              = 0.42 
;  res1@vpYF              =0.65

;  res1@gsnPaperHeight        = a4_height/cm_per_inch 
;  res1@gsnPaperWidth        = a4_width/cm_per_inch  
  res1@gsnPaperOrientation       = "portrait"
;  res1@gsnMaximize       = True
  res1@gsnPaperMargin    = 0.01
;  res1@tiXAxisString         = "~F8~q~F21~~B~v~N~, K";"O~B~3~N~, ppb";"Potential Temperature, K";
;  res1@xyExplicitLegendLabels = leg            ; explicit labels created above
  res1@xyDashPatterns    = (/0,0,16,14/)       ; 0 is default (solid)
  res1@xyLineColors           = (/"red","red","black","red","Magenta2","red","brown"/)
;  res1@xyMarkLineModes   = (/"MarkLines","Lines","Lines","Lines","Lines"   /)
  res1@xyMarkLineModes   = (/"Markers","Lines","Lines","Lines","Lines"   /)
 res1@xyMarkers         =  (/16,16/)                      ; choose type of marker  
; res1@xyMarkerColor     = "Turquoise3"       
 res1@xyMarkerColor     = "Cyan3"       
 res1@xyMarkerColor     = "red"       
 res1@xyMarkerSizeF     = 0.00855  
  res1@xyLineThicknessF       = 8 
  res1@pmLegendSide           = "Top"               ; Change location of
  res1@pmLegendParallelPosF   = .73                  ; move units right
  res1@pmLegendOrthogonalPosF = -1.1
  res1@pmLegendWidthF         = 0.12                ; Change width and
  res1@pmLegendHeightF        = 0.15                ; height of legend.
  res1@tiXAxisFuncCode = "~"
;  res1@tmXLMode = "Manual"
;  res1@trXMaxF = Tpotential_ea_sitesAve(1, 20) 
  res1@trYMaxF = 3. 
  if (ihours.eq.29) then 
   res1@trXMaxF = 310 
;   res1@trXMaxF = 330 
   res1@trXMinF = 301 
  end if 
  res1@tmYLValues     = fspan(0,5,11)
  res1@tmYLMinorValues= fspan(0,5,21)
  res1@tmYRMinorValues= fspan(0,5,21)
  res1@tmYLLabels     = fspan(0,5,11) 
  res1@tmYLMode = "Explicit"
  res1@trYMinF = 0
  res1@tmYLLabelStride = 1
  res1@tmYRLabelStride = 2
;  res1@tmXBLabelStride = 2
  res1@tmXTLabelStride = 2
  res1@tmYLFormat    = "f" 
  res1@tmYRFormat    = "f" 
  res1@tmXBFormat    = "f" 
  res1@tmYLTickSpacingF = 0.5
  res1@tmXBTickSpacingF = 5.0
  res1@tmYRMajorOutwardLengthF = 0
  res1@tmYLMajorOutwardLengthF = 0
  res1@tmXBMajorOutwardLengthF = 0
  res1@tmXBMinorOutwardLengthF = 0
  res1@tmYLMinorOutwardLengthF = 0
  res1@tmXTMinorOutwardLengthF = 0
  res1@tmYRMinorOutwardLengthF = 0
  res1@tmYLMinorOn = True 
  res1@tmYRMinorOn = True 
  res1@tmXBMinorOn = True 
  res1@tmXTMinorOn = True 

  res1@lgPerimOn=False
  res1@lgLabelFontHeightF = 0.0205
  res1@tiXAxisFontHeightF= 0.028
  res1@tiYAxisFontHeightF= 0.028
  res1@tmXBLabelFontHeightF = 0.028
  res1@tmYLLabelFontHeightF = 0.028

  res1@lgLabelFontHeightF = 0.0205
  res1@tiXAxisFontHeightF= 0.0128
  res1@tiYAxisFontHeightF= 0.0128
  res1@tmXBLabelFontHeightF = 0.0128
  res1@tmYLLabelFontHeightF = 0.0128
  res1@tmXMajorGrid                = True          ; implement x grid
  res1@tmXMajorGridThicknessF      = 1.0           ; 2.0 is default
  res1@tmXMajorGridLineDashPattern = 2             ; select short dash lines


  res1@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
;  res1@tmXTOn = False            ; turn off top   labels
;  res1@tmYROn = False            ; turn off right labels
;  res1@trYReverse = True
  res1@pmLegendDisplayMode    = "Never"       ; create a legend
;  res1@xyExplicitLegendLabels =  (/"05 UTC","09 UTC"/)

;  res1@lgItemOrder  = (/3,2,1,0/)
  res2=res1
  res2@xyMarkLineModes   = (/"Markers","Lines","Lines","Lines","Lines"   /)
 res2@xyMarkerColor     = "black"       
  res2@xyLineColors           = (/"green","red","black","red","Magenta2","red","brown"/)
  res2@xyDashPatterns    = (/16,0,16,14/)       ; 0 is default (solid)
;  res2@tmXTMode = "Explicit"
;  res2@tmXTValues     = ispan(0,20,2)
;  res2@tmXTLabels     = ispan(0,20,2) 
  res2@tmYRLabelStride = 2
  res2@tmYRTickSpacingF = 0.5
  res2@tiXAxisString = ""
  res2@tiYAxisString = ""
  res2@tmYRMode = "Manual"

   wks_type = "png"
   wks_type@wkWidth = 811
   wks_type@wkHeight = 811 
;   wks_type@wkWidth = 1011
;   wks_type@wkHeight = 1011 
   wks = gsn_open_wks(wks_type,figurename)
;  wks = gsn_open_wks("eps" ,figurename)          ; eps,pdf,x11,ncgm,eeps


  res2@tiDeltaF = 0.5
;  res1@tiMainOffsetYF=-0.01
  res2@tmYRLabelsOn = False
  res2@tmXTFormat    = "f" 
;  print(Tpotential_ea_sitesAve(:,:3))


   print("gonna reading "+soundingfiles(ihours))
   system("dos2unix "+soundingfiles(ihours))
   res2@xyExplicitLegendLabels =  (/"Radiosonde"/)
   data_second = readAsciiTable(soundingfiles(ihours), 12, "float", 7)
   print("finish reading table "+soundingfiles(ihours))

  file_char = stringtocharacter(soundingfiles(ihours))
;   print(file_char)
;  res1@xyExplicitLegendLabels =  chartostring(file_char(9:13)) 
   data_second@_FillValue = -999.
;   Height_second = sin(data_second(:,4)*3.1415926535/180.)*data_second(:,6)
  data_all = asciiread(soundingfiles(ihours), -1, "float")
 height_site = data_all(4)
  print("height_site="+height_site)
 delete(data_all)
   Height_second = (data_second(:,11)-height_site)/1000.

   index_bad = ind(Height_second.eq.0)
;     print(Height_second(index_bad))
   if(.not.all(ismissing(index_bad)))
      Height_second(index_bad) = Height_second@_FillValue
   end if 
   delete(index_bad)
   Press = data_second(:,2)
   Press = where(Press.ne.0, Press, Press@_FillValue)
   Theta_second = (data_second(:,1)+273.15)*(1000./Press)^0.286
;   print(Theta_second(:10))
;   print(data_second(:10,1))
;;;;;;;;;;;;;;;;;;;;diagnose PBLH from LBand Radiosound;;;;;;;;;;;;;;;;;;;;;;;;;;
     minPT_Obs = min(Theta_second)
      ilevel = minind(Theta_second)+1
          ismissing_xhu = all(ismissing(Theta_second(ilevel)))
           do while (ismissing_xhu)
 ;            print("ilevel "+ilevel +" "+Theta_second(ilevel))
             ilevel = ilevel + 1
             ismissing_xhu = all(ismissing(Theta_second(ilevel)))
           end do ;
 ;       if (icases.eq.15) then
 ;            print("ilevel "+ilevel +" "+Theta_second(ilevel))
 ;            print("(Theta_second(ilevel)-minPT_Obs)=" + (Theta_second(ilevel)-minPT_Obs))
 ;       end if
       do while ((Theta_second(ilevel)-minPT_Obs).lt.1.5)
        ilevel = ilevel + 1
          ismissing_xhu = all(ismissing(Theta_second(ilevel)))
           do while (ismissing_xhu)
 ;            print("ilevel "+ilevel +" "+Theta_second(ilevel))
             ilevel = ilevel + 1
             ismissing_xhu = all(ismissing(Theta_second(ilevel)))
           end do ;
       end do
       PBLH_1_5_Obs = Height_second(ilevel-1) ;+ (1.5-(Theta_second(ilevel-1)-minPT_Obs))/(Theta_second(ilevel)-Theta_second(ilevel-1))* (Height_second(ilevel)-Height_second(ilevel-1))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   q_sounding = mixhum_ptrh (Press, data_second(:,1)+273.15, data_second(:,3),-1)
;   print(data_second(:10,3))
;  res1@trXMaxF = 330. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) 
  res1@trXMinF = 290.; max(Tpotential_ea_sitesAve(1, :29)) 
  res1@trXMaxF = 320. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) 
  res1@trXMinF = 258.; max(Tpotential_ea_sitesAve(1, :29)) 
  res1@trXMaxF = 292. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) 
  res1@trYMaxF = 2. 
  res2@trYMaxF = 2.
;  res2@trXMaxF = max(Tpotential_ea_sitesAve(1, :29))
;  res2@trXMinF = min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/))
   res1@vpXF             = 0.1
  res1@tiXAxisString         = "Potential temperature, K";"Potential Temperature, K";
  res1@tiYAxisString         = "Height, km"; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3"      ; label bottom axis with units
  res1@pmLegendDisplayMode    = "Never"       ; create a legend
  res2@pmLegendDisplayMode    = "Never"       ; create a legend

   res1@tiMainString = chartostring(file_char(15:24)) +" at " + chartostring(file_char(9:13)) + " Shenyang_Liaoning"
  if (file_char(9:13).eq."55299") then 
    res1@trXMinF = 310.; max(Tpotential_ea_sitesAve(1, :29)) 
    res1@trXMaxF = 345. ; min((/Tpotential_ea_sitesAve(1, 0)-6,TH2_Mesonet_mean-0.5/)) 
  end if 
   plot                      = gsn_csm_xy(wks,Theta_second,Height_second,res1)
;    gsn_polyline (wks, plot, (/0,1000/), (/PBLH_1_5_Obs,PBLH_1_5_Obs/), resp)
    resp@gsLineColor      = "red"                     ; color of lines
;    deltaPBLH = min((/abs(PBLH-PBLH_1_5_Obs), abs(PBLH_1_5-PBLH_1_5_Obs)/))
   res1@vpXF             = 0.53
   res1@trXMaxF = 100. 
   res2@trXMaxF = 100. 
   res1@trXMinF = 10.02
   res2@trXMinF = 10.02
   res1@tmYLLabelsOn = False
   res2@tmXTLabelsOn = True 
;   res1@tiXAxisString         = "Mixing ratio, g kg~S~-1"
   res1@tiXAxisString         = "RH, %"
  res1@tiYAxisString         = ""; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3"      ; label bottom axis with units
  res1@pmLegendDisplayMode    = "Never"       ; create a legend
  res2@pmLegendDisplayMode    = "Never"       ; create a legend
;   print(q_sounding(:))
;   print(Height_second(:))
;   print(Height_second)
;   print(q_sounding)
   res1@xyMarkerColor     = "blue"       
   res1@tiMainString ="" ; chartostring(file_char(15:24)) +" at " + chartostring(file_char(9:13))
;   plot                      = gsn_csm_xy(wks,q_sounding,Height_second,res1)
   plot                      = gsn_csm_xy(wks,data_second(:,3),Height_second,res1)
;    gsn_polyline (wks, plot, (/0,1000/), (/PBLH_1_5_Obs,PBLH_1_5_Obs/), resp)
   res1@tmYLLabelsOn = True 
   delete(data_second)
   delete(Height_second)
   delete(Press)
   delete(Theta_second)
   delete(q_sounding)
      gsres = True
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerIndex = 13 
      gsres@gsMarkerSizeF = 0.015
      gsres@gsMarkerColor      = "pink" 
;    print(TH2_Mesonet_mean)
;      gsn_polymarker(wks,plot,TH2_Mesonet_mean,0.002,gsres)


  frame(wks)


;  system("convert -trim -resize 100% "+figurename+".eps /nsftor/xhu/public_html/BeijingO3_MaZhiQiang/Simulations/3.6/wrfYMiao_YSU.20100912/"+figurename+".png")
;  system("rm "+figurename+".eps")
;  print("finish "+figurename+".eps")
end do  ; ihour

 end