;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Plot data on a cross section
;   This script will plot data at a set angle through a specified point
;   Add some label info to the Y-axis
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "/home/xhu2/tools/nclscripts/wrf/WRFUserARW.ncl"

begin

  elevation_array = readAsciiTable("ACTAMERICA-Elevation_B200_20160816_R0_L1.ict", 4, "float", 36)
  elevation_array@_FillValue = -9999.
;  Altitude_AGL_m  = elevation_array(:,3)
  Altitude_GPS    = elevation_array(:,1)
  Altitude_time    = elevation_array(:,0)

  hours_since_midnight_UTC_simple = round(elevation_array(:,0)/3600.*10, 3)/10.

  
  file_data = "ACTAMERICA-PDS_B200_20160816_R1_L1.ICT"
 if (file_data.eq."ACTAMERICA-PDS_B200_20160711_R0.ict".or.file_data.eq."ACTAMERICA-PDS_B200_20160715_R1.ict".or.file_data.eq."ACTAMERICA-Hskping_c130_20160721_R0.ict") then
  data_array = readAsciiTable(file_data, 33, "float", 65)
 else
  data_array = readAsciiTable(file_data, 34, "float", 67)
 end if
  data_array@_FillValue = -9999.

  lat = data_array(:,2)
  lon = data_array(:,3)
    delete(data_array)
    data_array = readAsciiTable("ACTAMERICA-PICARRO_B200_20160816_R0_L1.ict", 6, "float", 38)
    data_array@_FillValue = -9999.
    CO2_obs_raw_temp = data_array(:,3)
   index_good = ind(.not.ismissing(CO2_obs_raw_temp))
    CO2_obs_raw = CO2_obs_raw_temp(index_good) 
    CO2_time_raw = data_array(index_good,2)
;  print(CO2_obs_raw)
   delete(index_good)

    CO_obs_raw_temp = data_array(:,5)
   index_good = ind(.not.ismissing(CO_obs_raw_temp))
    CO_obs_raw = CO_obs_raw_temp(index_good)*1000. ; convert ppm to ppb
    CO_time_raw = data_array(index_good,2)
   delete(index_good)
;  print(CO_obs_raw)
    delete(data_array)
    data_array = readAsciiTable("ACTAMERICA-Ozone_B200_20160816_R0_L1.ict", 4, "float", 36)
    data_array@_FillValue = -9999.
    O3_obs_raw_temp = data_array(:,3)
   index_good = ind(.not.ismissing(O3_obs_raw_temp))
    O3_obs_raw = O3_obs_raw_temp(index_good)
    O3_time_raw = data_array(index_good,2)
   delete(index_good)


      CO2_obs = linint1 (CO_time_raw, CO_obs_raw, False, Altitude_time, 0)
     print(CO2_time_raw(:10))
     print(Altitude_time(:10))

    if (dimsizes(Altitude_GPS).ne.dimsizes(CO2_obs)) then 
      print("records in two files doesn't match !!!!!!!!!!!!!!!!" + "data line "+dimsizes(CO2_obs)+" elevation line "+dimsizes(Altitude_GPS))
    end if 

      gsres               = True
      gsres@gsMarkerIndex = 16          ; circle at first
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.01
      gsres@gsMarkerColor      = "black"
      gsres@gsMarkerColor      = "red"
; 32.453037,-93.834152 Shreveport Regional Airport
     Meteo_Lat = 32.453037 
     Meteo_Lon = -93.834152 
idomain =  1
;   files = systemfunc("ls wrfout_d0"+idomain+"_2011-08-??_*:00:00")
   files = systemfunc("ls wrfout_d01_2016-08-16_1[5-9]* wrfout_d01_2016-08-16_20*:00") 

do ifiles =  0 ,    dimsizes(files)-1
  a = addfile(files(ifiles)+".nc","r")
  opt = True
  ij = wrf_user_ll_to_ij(a,Meteo_Lon, Meteo_Lat,opt)
;  print(ij)

; We generate plots, but what kind do we prefer?
 type = "png"
 type@wkWidth = 998
 type@wkHeight = 998

  figurename = "wrfout_d0"+idomain+"_CO2total_crossSNpartial_overlayObs_B200_co_"+ifiles
  wks = gsn_open_wks(type,figurename)
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map


; Set some basic resources
  res = True
;  res@MainTitle = "REAL-TIME WRF"
  res@Footer = False
  


  ter_res = True
  opts_ter = ter_res
  opts_ter@gsnYRefLine = 0.0
  opts_ter@gsnDraw = False
  opts_ter@gsnFrame = False

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  ntimes = dimsizes(times)          ; number of times in the file
  FirstTime = True

  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
  nd = dimsizes(mdims)

;---------------------------------------------------------------

  do it = 0,ntimes-1,2                  ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)           ; Set Valid time to use on plots
  ter  = wrf_user_getvar(a,"ter",0)  
  xlon = wrf_user_getvar(a, "XLONG",0)
  xlat = wrf_user_getvar(a, "XLAT",0)
  PBLH  = a->PBLH(0,:,:)
    tc  = wrf_user_getvar(a,"tc",it)     ; T in C
    theta = wrf_user_getvar(a,"theta",it)      ; relative humidity
    z   = wrf_user_getvar(a, "z",it)     ; grid point height
  va  = wrf_user_getvar(a,"va",it)   ; u on mass points
  temp  = a->CO2_BIO(it,:,:,:)   ; u on mass points
  temp  = (/a->CO2_BIO(0,:,:,:) + a->CO2_ANT(0,:,:,:) - a->CO2_BCK(0,:,:,:)/) 
  CO2total = temp 
  CO2total@units = "ppmv";(/temp*1000/) 
  wa  = wrf_user_getvar(a,"wa",it)   ; v on mass points

    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = max(z)/1000.
      nz   = floattoint(zmax/2 + 1)
      FirstTime = False
    end if
;---------------------------------------------------------------
; Plot a cross session that run south-north through the middle of the plot
; For this we need a pivot point and a angle
;                   |
;       angle=0  ; is  |

;        angle = 90  ; E-W
        angle = -14 
        plane = new(2,float)
        plane = (/ ij(0), ij(1) /)    ; pivot point is center of domain (x,y)
        opts = False                                  ; start and end points not specified

        theta_plane = wrf_user_intrp3d(theta,z,"v",plane,angle,opts)
        tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
        CO2total_plane = wrf_user_intrp3d(CO2total,z,"v",plane,angle,opts)
        va_plane = wrf_user_intrp3d(va,z,"v",plane,angle,opts)
        wa_plane = wrf_user_intrp3d(wa,z,"v",plane,angle,opts)
        ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
        PBLH_plane = wrf_user_intrp2d(PBLH,plane,angle,opts)
        X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
        lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
   distance = (lon_plane( :)+96.8)^2+(X_plane(:)-32.7833333)^2
   id_d = minind(distance)

        dim = dimsizes(theta_plane)                      ; Find the data span - for use in labels
       print(dim)
        zspan = dim(0)
        print(zspan)
         zcordinate = fspan(zmin,zmax,zspan)

      ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@cnMissingValPerimOn     = False ; True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz*2*2)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz*2*2))  ; Create labels
        opts_xy@tmXBMode                = "Explicit"
        X_plane_round = round(X_plane*10,3)/10.
        opts_xy@tmXBValues              = ispan(0,260,20)
        opts_xy@tmXBLabels              = X_plane_round(ispan(0,260,20)) 

        opts_xy@tiXAxisFontHeightF      = 0.020
        opts_xy@tiYAxisFontHeightF      = 0.020
        opts_xy@tmXBMajorLengthF        = 0.02
        opts_xy@tmYLMajorLengthF        = 0.02
        opts_xy@tmYLLabelFontHeightF    = 0.015
        opts_xy@PlotOrientation         = tc_plane@Orientation
        n_YMax = 400 ; domain 1
        n_XMax = 180; 260 ; id_d+1 ; domain 1
;        n_XMax =  260 ; id_d+1 ; domain 1
        n_XMin = 30  ; domain 1
;        n_XMin = 0  ; domain 1

        opts_xy@trYMaxF = n_YMax 
        opts_xy@trXMaxF = n_XMax 
        opts_xy@trXMinF = n_XMin 
;        opts_xy@gsnOrientation = "portrait"
  	opts_xy@gsnDraw  = True    ; Forces the plot to be drawn
  	opts_xy@gsnFrame = True    ; Frame advance
        opts_xy@tiMainString = "" ;chartostring(a->Times)
         opts_xy@vpHeightF        = 0.35
         opts_xy@vpWidthF         = 0.5
;        opts_xy@gsnRightString = chartostring(a->Times)
        

      ; Plotting options for CO2total
        opts_theta = opts_xy
        opts_theta@cnFillMode                =  "RasterFill"
        opts_theta@ContourParameters       = (/ 396., 406., 1. /)
        opts_theta@tmXBFormat =  "f"
        opts_theta@pmLabelBarOrthogonalPosF = -0.07
        opts_theta@cnFillOn                = True
        opts_theta@gsnSpreadColors       = True             ; use full range of colormap
;        opts_theta@cnFillColors            = (/"White","White","White", \
;                                            "White","Chartreuse","Green", \
;                                            "Green3","Green4", \
;                                            "ForestGreen","PaleGreen4"/)

      ; Plotting options for Temperature
        opts_tc = opts_xy
        opts_tc@cnInfoLabelOrthogonalPosF = 0.00
        opts_tc@ContourParameters  = (/ 5. /)


      ; Get the contour info for the theta and temp
        opts_theta@gsnMaximize = False 

       CO2total_plane!0 = "height"
       CO2total_plane!1 = "lat"
       CO2total_plane&height = fspan(0,zmax,dim(0) ) 
       CO2total_plane&lat    = X_plane  

       va_plane!0 = "height"
       va_plane!1 = "lat"
       va_plane&height = fspan(0,zmax,dim(0) ) 
       va_plane&lat    = X_plane  

       wa_plane!0 = "height"
       wa_plane!1 = "lat"
       wa_plane&height = fspan(0,zmax,dim(0) ) 
       wa_plane&lat    = X_plane  

   
        contour_theta = wrf_contour(a,wks,CO2total_plane,opts_theta)
     j_Altitude_GPS  = elevation_array(:,1)/(zmax*1000./dim(0)) 
     i_lat = (lat-min(X_plane))/((max(X_plane)-min(X_plane))/dimsizes(X_plane))
     print("max, min of X_plane"+ min(X_plane)+" "+max(X_plane)+"max, min of lat "+min(lat)+" "+max(lat)  )
     print("max, min of i_lat"+ min(i_lat)+" "+max(i_lat)  )
;      dummy1 = gsn_add_polymarker(wks,contour_theta,i_lat,j_Altitude_GPS,gsres)

  npts = dimsizes(lat)
;  print(contour_theta@contour)
;  print(contour_theta)
;  getvalues contour_theta@contour
;    "cnFillColors" : colors
;  end getvalues
;    levels = ispan(396, 406, 1)
    levels = ispan(80, 140, 5)
;     levels = ispan(46, 80, 2) ; for O3 overlay

  num_distinct_markers = dimsizes(levels)+1        ; number of distinct markers
    colors = round(fspan(2,201, num_distinct_markers),3) 
  print(colors)
  dims=dimsizes(CO2_obs) 
  y_new = new((/num_distinct_markers,dims/),float,-999)
  x_new = new((/num_distinct_markers,dims/),float,-999)
  do j = 0, num_distinct_markers-1
    if (j.eq.0) then
      indexes = ind(CO2_obs.lt.levels(0))
    end if
    if (j.eq.num_distinct_markers-1) then
      indexes = ind(CO2_obs.ge.max(levels))
;     print("color "+j+" indexes="+indexes)
    end if
    if (j.gt.0.and.j.lt.num_distinct_markers-1) then
      indexes = ind(CO2_obs.ge.levels(j-1).and.CO2_obs.lt.levels(j))
    end if

    if (.not.any(ismissing(indexes))) then
      npts_range = dimsizes(indexes)   ; # of points in this range.
      y_new(j,0:npts_range-1) = j_Altitude_GPS(indexes)
      x_new(j,0:npts_range-1) = i_lat(indexes)
     print("color "+j+"npts_range="+npts_range)
    end if
    delete(indexes)            ; Necessary b/c "indexes" may be a different
                               ; size next time.
  end do ; j = 0, num_distinct_markers-1
  gsres               = True
  gsres@gsMarkerIndex = 4          ; circle at first

  do j = 0, num_distinct_markers-1
    if (.not.ismissing(y_new(j,0)))
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.014
      gsres@gsMarkerColor      = "black"
    if (j.eq.0) then 
      dummy1 = gsn_add_polymarker(wks,contour_theta,i_lat,j_Altitude_GPS,gsres) 
        dummy = new(num_distinct_markers,typeof(dummy1))
    end if 
       gsres@gsMarkerIndex = 16          ; Use filled dots for markers.
      gsres@gsMarkerSizeF = 0.011
      gsres@gsMarkerColor      = colors(j)
;      if (j.eq.0) then 
;        dummy1 = gsn_add_polymarker(wks,contour_theta,x_new(j,:),y_new(j,:),gsres)
;         printVarSummary(dummy1)
;      end if 
      dummy(j) = gsn_add_polymarker(wks,contour_theta,x_new(j,:),y_new(j,:),gsres)
    end if
  end do

 delete(y_new)
 delete(x_new)
       txres               = True                         
       txres@txFontHeightF = 0.01 
      text = gsn_add_text(wks,contour_theta,flt2string(hours_since_midnight_UTC_simple(::1500)),i_lat(::1500), \
                                    j_Altitude_GPS(::1500)+12 ,txres)

        
  	resvector = True
         resvector@vpHeightF        = 0.35
         resvector@vpWidthF         = 0.5
  	resvector@gsnDraw  = True    ; Forces the plot to be drawn
  	resvector@gsnFrame = True    ; Frame advance
  	resvector@vcGlyphStyle = "CurlyVector" 
  	resvector@tmXBFormat =  "f"
;        resvector@trXMinF = 65.
        resvector@trYMaxF = n_YMax 
        resvector@trXMaxF = n_XMax 
        resvector@trXMinF = n_XMin 
        resvector@vcRefAnnoOn = True
   	resvector@vcRefMagnitudeF          = 8.             ; define vector ref mag
   	resvector@vcRefLengthF             = 0.025          ; define length of vec ref
   	resvector@vcRefAnnoOrthogonalPosF  = -1.            ; move ref vector
   	resvector@vcRefAnnoParallelPosF  = 0.95            ; move ref vector
   	resvector@vcMinDistanceF           = 0.025            ; larger means sparser
   	resvector@vcLineArrowHeadMaxSizeF  = 0.0075          ; default: 0.05 (LineArrow), 0.012 (CurlyVector)
        resvector@gsnMaximize = False 

  
        vector = wrf_vector(a,wks,va_plane,wa_plane*100,resvector)


;Contour terrain cross section
;      print(opts_ter)
    opts_ter@trXMaxF = X_plane(n_XMax) 
    opts_ter@trXMinF = X_plane(n_XMin) 
    opts_ter@trYMaxF = zcordinate(n_YMax)*1000 
         opts_ter@vpHeightF        = 0.35
         opts_ter@vpWidthF         = 0.5
  opts_ter@xyLineThicknessF = 4.0
  opts_ter@pmLegendSide           = "Top"               ; Change location of
  opts_ter@pmLegendParallelPosF   = .1084                  ; move units right
  opts_ter@pmLegendOrthogonalPosF = -0.286
  opts_ter@pmLegendWidthF         = 0.081                ; Change width and
  opts_ter@pmLegendHeightF        = 0.081                ; height of legend.
  opts_ter@pmLegendDisplayMode    = "Always"       ; create a legend
  opts_ter@lgPerimOn     = False
  opts_ter@lgLabelFontHeightF = 0.0148
  opts_ter@xyExplicitLegendLabels = (/"terrain","PBL top","lfc","lcl"/)
  opts_ter@xyLineColors       = (/"blue","black","blue"/)
  opts_ter@xyDashPatterns = (/0, 12, 16/)
  opts_ter@gsnYRefLine = 0.0
;  opts_ter@gsnAboveYRefLineColor = "black"

        ter_plane_mean_plusref = new ((/2,dimsizes(ter_plane)/),float)
        ter_plane_mean_plusref(0,:) = (/ter_plane/)
        ter_plane_mean_plusref(1,:) = (/ter_plane+PBLH_plane/)
;        contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
        contour_ter = gsn_csm_xy(wks,X_plane,ter_plane_mean_plusref,opts_ter)
;        contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
;        dummy1 = gsn_add_polymarker(wks,contour_ter,X_plane(id_d),0.08,gsres)

;     j_Altitude_GPS  = elevation_array(:,1)/(zmax*1000./dim(0)) 
;     i_lat = (lat-min(X_plane))/((max(X_plane)-min(X_plane))/dimsizes(X_plane))
;     print("max, min of X_plane"+ min(X_plane)+" "+max(X_plane)+"max, min of lat "+min(lat)+" "+max(lat)  )
;     print("max, min of i_lat"+ min(i_lat)+" "+max(i_lat)  )
;        dummy1 = gsn_add_polymarker(wks,contour_ter,lat,Altitude_GPS,gsres)


      ; MAKE PLOTS         
  pltres = True
  pltres@NoTitles = True
  pltres@CommonTitles = True; False 
  pltres@PanelPlot = True
  pltres@gsnPaperOrientation = "portrait"
  pltres@FramePlot = False

;        plot = wrf_overlays(a,wks,(/contour_theta,vector/),pltres)
        plot = wrf_overlays(a,wks,(/contour_theta,vector,contour_ter/),pltres)
      setvalues contour_theta
;         "vpHeightF"        : 0.35
;         "vpYF"             : 0.9-ifile_loop*0.39
;         "vpWidthF"         : 0.35
;         "vpXF"             : 0.15+ispecies*0.355
         "tiMainString"     : res@TimeLabel + " UTC";"Cross section"
         "tiMainOffsetYF"   : -0.01
;         "gsnRightString"   : "RString" ; res@TimeLabel
;         "gsnPaperOrientation" : "portrait"
       end setvalues
     draw(contour_theta)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

  delete(ter_plane)
  delete(wa_plane)
  delete(tc_plane)
  delete(theta_plane)
  delete(va_plane)
  delete(CO2total_plane);  = wrf_user_intrp3d(CO2total,z,"v",plane,angle,opts)
  delete(zcordinate)
  lbres                      = True
  lbres@vpHeightF            = 0.25
  lbres@vpWidthF             = 0.03
  lbres@lbAutoManage         = False
; No margins around labelbar.
  lbres@lbBottomMarginF      = 0.0
  lbres@lbLeftMarginF        = 0.0
  lbres@lbRightMarginF       = 0.0
  lbres@lbTopMarginF         = 0.0

; Turn various features on and off.
  lbres@lbLabelsOn           = True
  lbres@lbPerimOn            = False;True
  lbres@lbTitleOn            = True
  lbres@lbTitlePosition  = "Top"              ; title position
  lbres@lbTitleFontHeightF   = 0.01258                ; make title smaller
  lbres@lbLabelFontHeightF   = 0.01258
  lbres@lbTitleString    = "co~C~~V10~ppb"
  lbres@lbTitleDirection = "Across"             ; title direction
  lbres@lbTitleOffsetF = -0.015
  lbres@lbTitleJust = "CenterCenter"
  lbres@lbMonoFillPattern    = True
  lbres@lbOrientation        = "vertical";"Horizontal"
  lbres@lbBoxMinorExtentF    = 1.0

; Which point to position labelbar about.
  lbres@lbJustification      = "TopLeft"

  lbres@lbLabelAlignment     = "InteriorEdges"
  lbres@lbLabelStride    = 2
;  lbres@lbLabelJust          = "CenterLeft"
  lbres@lbLabelJust          = "CenterCenter"
;  lbres@lbLabelOffsetF       = -0.1

  labels = levels + ""
  lbres@lbFillColors         = colors
  ypos = 0.15
  lbid = gsn_create_labelbar_ndc(wks,dimsizes(levels)+1,labels, \
                                   0.89,0.305,lbres)
    draw(lbid)

  frame(wks)

;    system("convert -density 300 -rotate -90 -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png")
;    system("convert -density 300  -resize 1000x1000 -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.6.1/YSU/wrfNARR3dWSM6_CONUS_UCM_YSU.2011080700/"+figurename+".png")
 print("finish "+figurename+".eps")
end do ; ifiles
end