;This ncl script will plot the DEC PM2.5 day
;
;begin
;************************************
       f = "42-089-0001-SO2-1990-2020.Final.out"

       data = asciiread(f,(/11323,17/),"double")
        nrow = data
        x =  data(:,0)   ; Year/day in Julian decimal
        x@long_name = "Date/Time"
        year = x

        y = data(:,14)  ; SO2 concentration
        y@long_name = "SO2 concentration (ppb)"
       y@_FillValue=integertoshort(-999)
        y=where(y.lt.0,y@FillValue,y)

        y@_Fillvalue = -999
        ny = dimsizes(y)

;-----------------------------------------------------
; plot resources
;-----------------------------------------------------
  plot = new(1,"graphic")

  wks  = gsn_open_wks("png","manken")        ; send graphics to PNG file

  res                     = True                   ; plot mods desired
  res@gsnDraw             = False                  ; don't draw yet
  res@gsnFrame            = False                  ; don't advance frame yet

  res@xyDashPatterns      = 0                      ; solid line
  res@xyLineColors        = (/"black", "red", "blue"/)
  res@xyLineThicknesses   = (/1,3,1/)
  res@tiMainString        = "f"                   ; title

  res@vpHeightF           = 0.4                    ; change aspect ratio of plot
  res@vpWidthF            = 0.8
  res@vpXF                = 0.1                    ; syrt plot at x ndc coord

  res@trXMinF  = year(0)
  res@trXMaxF  = year(ny-1)+1
  res@trYMinF  = -0.6
  res@trYMaxF  =  0.6
  res@gsnYRefLine      = 0.0

  txres               =  True
  txres@txFontHeightF = 0.02
  txres@txJust        = "CenterCenter"
  txres@txFontThicknessF = 2.0      ; default=1.00
  txres@txFontHeightF    = 0.025    ; default=0.05

;----------------------------------------------------
; For each partition, calculate
; .   Simple linear regression line
; .   Mann-Kendall and Thiel-Sen estimates
;-----------------------------------------------------
; Do not plot immediately
;-----------------------------------------------------

  dplt = new ( (/ny,nrow/), typeof(y), y@_FillValue)


     rc   = regline(year(ntStrt:ntLast),y(ntStrt:ntLast))
     p    = trend_manken(y(ntStrt:ntLast), False, 0)

     dplt = y@_FillValue
     dplt(0,ntStrt:ntLast) = y(ntStrt:ntLast)
;     dplt(1,ntStrt:ntLast) = rc*(year(ntStrt:ntLast)-rc@xave) + rc@yave
;     dplt(2,ntStrt:ntLast) = p(1)*(year(ntStrt:ntLast)-rc@xave) + rc@yave

     plot(0)  = gsn_csm_xy (wks,year,dplt,res)           ; create plot

     text          = "p="+sprintf("%5.3f",p(0))+"  trend="+sprintf("%5.3f",p(1))
     plot@$unique_string("dum")$ = gsn_add_text(wks,plot(0),text, 2000, 0.30 ,txres)

  draw(plot(0))
  frame(wks)