load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "tools/trimPNG.ncl"

begin
  ;; HOW TO USE
  ;; set the details in the namelist: "../av-tool/mynamelists/cp.ens-list"
  ;; set the input conditions below

  ;; Start of input conditions
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  var = "pr"
  xcm = "CPM"
  frq = "day"
  grd = "rcm"
  obs = "obs-merge"
  cond = "true"
  cgol = "false"
  ; true/false [CPM + cond only]
  ths = "0-100"
  syr = "all"      ; all/sel
  ; all/sel
  mapd = "manual"  ; "manual" or "auto"
  focus = "av"     ; "av" or "ccds/pav-near/far" 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; End of input conditions

  ; mask list
  subregs = (/      "CFR",      "SFR",      "CRS",    "SGR",      "CH", \
                 "NIT",   "CIT",   "SIT",    "EAL",  "EAC"/)
  obsregs = (/"COMEPHORE","COMEPHORE","COMEPHORE","RADKLIM","RdisaggH", \
              "EURO4M","GRIPHO","GRIPHO","EURO4M","EURO4M"/)
  numregs = dimsizes(subregs)

  if ( focus.ne."av" ) then
    dsdel = str_split(focus,"-")
    dstyp = dsdel(0)
    dstim = dsdel(1)
    delete(dsdel)
  end if    

  if ( syr .eq. "all" ) then
    syd = ""
    yrs = "all-years"
  else if ( syr .eq. "sel" ) then
    syd = "/selyr"
    yrs = "????-????"
  end if
  end if 
  if ( xcm .eq. "RCM" ) then
    fn = "1-3"
  else if ( xcm .eq. "GCM" ) then
    fn = "1-2"
  end if
  end if

  cdir = ""
  cnam = "nocond"
  if ( cond .eq. "true" ) then
    if ( xcm .eq. "CPM" ) then
      if ( cgol .eq. "false" ) then
        cdir = "cond-cro/"
        cnam = "cond-cro"
        ctitl = "CPM-RCM"
      else 
        cdir = "cond-cgo/"
        cnam = "cond-cgo"
        ctitl = "CPM-GCM"
      end if
    else  
      cdir = "cond-rgo/"
      cnam = "cond-rgo"
      ctitl = "RCM-GCM"
    end if
  end if

  if ( focus.ne."av" ) then
    cdir=""
    if ( xcm.eq."CPM" ) then
      if ( cgol.eq."true" ) then
        cnam = "cgo"
        ctitl = "CPM-GCM"
      else
        cnam = "cro"
        ctitl = "CPM-RCM"
      end if
    else
      cnam = "rgo"
      ctitl = "RCM-GCM"
    end if
  end if

  tdir = ""
  if ( ths .ne. "0-100" ) then
    tdir = "th-"+ths+"/"
  end if

  nlist = "../av-tool/mynamelists/cp.ens-list"
  bnum = 18
  bnum2 = 33 
  system("cat "+nlist+" | grep num_mod > nums.txt" )
  nmt = asciiread("nums.txt", (/1/), "string")
  nms = str_split(nmt,"'")
  nm = nms(1)
  nm := stringtointeger(nm)
  system("rm nums.txt") 

  names = new((/nm/),string)

;---Start the graphics
  system("mkdir -p images")
  outim = "images/box-"+focus+"-"+var+"-"+frq+"-"+xcm+"-"+grd+"-"+cnam+"-"+ths
  formo = "png"
  wks = gsn_open_wks (formo, outim)

  system("cat "+nlist+" | grep workdir > wd.txt" )
  wd = str_split(asciiread("wd.txt", (/1/), "string"),"'")
  wd := wd(1)
  system("rm wd.txt")
  if ( focus.ne."av" ) then
    if ( dstim.eq."near" ) then
      jsrch = "job_near"
    end if
    if ( dstim.eq."far" ) then
      jsrch = "job_far"
    end if
  else
    jsrch = "job_master"
  end if
  system("cat "+nlist+" | grep "+jsrch+" > jm.txt" )
  jm = str_split(asciiread("jm.txt", (/1/), "string"),"'")
  jm := jm(1)
  system("rm jm.txt")
  hdir = wd+"/"+jm+"/"
  if ( focus.eq."av" ) then
    remaskd = "/remap_"+grd+"/mask_"+obs
    avdsd = "av/"
  else
    remaskd = "/remap_"+grd
    if ( dstyp.eq."pav" ) then
      avdsd = "ccds/ads/pav/"
    else
      avdsd = "ccds/"
    end if
  end if
  idir = hdir+xcm+"/"+frq+syd+remaskd+"/pdfs/"+cdir+tdir+"rdif/"+avdsd

  true_nm = nm
  c = 0
  do n=1, nm
    nline = bnum+n
    mline := asciiread(nlist, (/nline/), "string")
    mline := mline(bnum+n-1)
    modn := str_split(mline,"'")
    modn := modn(1)
    ;print("modn = "+modn)
    modn_sub = str_split(modn,"_")
    cpmn = modn_sub(3)
    rcmn = modn_sub(2)
    if ( cgol .eq. "true" ) then
      rcmn = "none"
    end if
    gcmn = modn_sub(0)+" ("+modn_sub(1)+")"
    gcmo = modn_sub(0)
    ensn = modn_sub(1)
    if ( xcm .eq. "RCM" .or. xcm .eq. "GCM" ) then
      system("echo "+modn+" | cut -d_ -f"+fn+" > modnam.txt" )
      modn2 := asciiread("modnam.txt", (/1/), "string")
      modn = modn2
      system("rm modnam.txt")
    end if
    if ( xcm .eq. "CPM" ) then
      modn = gcmo+"_"+ensn+"_"+rcmn+"_"+cpmn+"-CP" ;modn+"-CP"
    end if
    print("modn = "+modn)

    system("ls "+idir+var+"_"+modn+"_"+frq+"_"+yrs+".nc 2> err.t | wc -l > num.t")
    ncheck = asciiread("num.t", (/1/), "string")
    system("rm num.t err.t")
    if ( ncheck .eq. "0" ) then
      print("skipped")
      true_nm = true_nm -1
      delete(ncheck)
      continue
    end if
    if ( xcm .eq. "CPM" .and. cgol .eq. "false" .and. rcmn .eq. "none" ) then
      print("skipped")
      true_nm = true_nm -1
      delete(ncheck)
      continue
    end if
    c = c+1
 
    inf = idir+var+"_"+modn+"_"+frq+"_"+yrs+".nc"
    system("ls "+inf+" > plot-tmp.txt") 
    ff = asciiread("plot-tmp.txt", (/1/), "string")
    system("rm plot-tmp.txt")
 
    fid   = addfile (ff , "r")
    v    = fid->$var$

    v@lat2d = fid->lat
    v@lon2d = fid->lon
    pv = v(0,:,:)
    lat2d = fid->lat
    lon2d = fid->lon
    nlat = dimsizes(lat2d(:,0))
    nlon = dimsizes(lat2d(0,:))

    if ( c.eq.1 )
      pvs = new((/nm,nlat,nlon/),float)
    end if

    dimsz = dimsizes(v)
    ndims = dimsizes(dimsz)

    pvs(n-1,:,:) = pv(:,:)
    delete(pv)
  end do

  xpos = tofloat(ispan(1,numregs,1))
  ypos = new((/numregs,5/),float)
  ymax = new((/numregs/),float)
  yavg = ymax
  ymin = ymax

  ; apply masks
  remaskd = "/"+frq+"/remap_"+grd
  do s=0, numregs-1
    sr  = subregs(s)
    obs = obsregs(s)
    print("Mask: "+sr)
    inf = hdir+"OBS/"+obs+remaskd+"/sub_mask_"+sr+"/"+sr+"_mask.nc"
    fid = addfile(inf, "r")
    m  := fid->mask
    mdim = dimsizes(m)
    mmdim = dimsizes(mdim)
    if ( mmdim.gt.2 ) then
      m := m(0,:,:)
    end if
   
    ; mask data
    do n=0 , nm-1
      mpv = where(m.gt.0,pvs(n,:,:),pvs@_FillValue)
    end do
    statp = stat_dispersion(mpv,False)
    ymin(s)   = statp(2)  ; minimum
    ypos(s,0) = statp(23) ; lower 1%       - 1th percentile
    ypos(s,1) = statp(6)  ; lower quartile - 25th percentile
    ypos(s,2) = statp(8)  ; median
    yavg(s)   = statp(0)  ; mean
    ypos(s,3) = statp(10) ; upper quartile - 75th percentile
    ypos(s,4) = statp(26) ; upper 1%       - 99th percentile
    ymax(s)   = statp(14) ; maximum
    delete(statp)
    delete(mpv)
    delete(m)
  end do

;---Set some plotting options
  res = True
; res@gsnMaximize         = True
  res@vpWidthF      = 0.8
  res@vpHeightF     = 0.4

; res@tiYAxisOn    = True
; res@tiMainOn    = True
; res@tiMainString = focus+" "+ctitl+" ("+frq+" "+ths+") N="+c
  res@tiMainFontHeightF = 0.015
  res@tmXBLabels = subregs
  res@tmXBLabelFontHeightF = 0.015
  if ( focus .eq. "av" ) then
    ylab = "A.V."
  else 
    if ( dstyp .eq. "ccds" ) then
      ylab = "C.C.D.S."
    else if ( dstyp .eq. "pav" ) then
      ylab = "P.A.V."
    end if
    end if
  end if 
  res@tiYAxisString = ylab
  res@tmYLFormat    = "@f.1" 
  res@trYMinF = -.5
  res@trYMaxF =  .5

  res@tfPolyDrawOrder = "PostDraw"
; res@tiXAxisOn    = False
; res@tmXTLabelsOn = False
; res@tmXBLabelsOn = True
; res@tmYLLabelsOn = True
; res@tmXTOn = False
; res@tmXBOn = True  ; Turn off tickmarks and labels
; res@tmYLOn = True
; res@tmYROn = False
; res@pmTickMarkDisplayMode = "Always"            ; turn on built-in tickmarks
; res@tmYLLabelsOn = True
; res@tmYLOn = True
; res@tmYROn = False

  boxOpts = True
  boxOpts@gsLineDashPattern = 0
; boxOpts@boxColors = "green"

  linres = True
  linres@gsLineDashPattern = 0

  plot = boxplot(wks,xpos,ypos,boxOpts,res,linres)

  uptfill = NhlNewMarker(wks,"u",34,0.0,0.0,1.0,1.0,0.0)
  dotfill = NhlNewMarker(wks,"u",34,0.0,0.0,1.0,1.0,180.0)
  mxres = True
  mxres@gsMarkerIndex = dotfill
  mxres@gsMarkerSizeF = 9.0
  mxres@gsMarkerColor = "black"
  maxmk = gsn_add_polymarker(wks,plot,xpos,ymax,mxres)

  mnres = True 
  mnres@gsMarkerIndex = uptfill  
  mnres@gsMarkerSizeF = 9.0
  mnres@gsMarkerColor = "black"
  minmk = gsn_add_polymarker(wks,plot,xpos,ymin,mnres)

  avgln = new((/numregs/),graphic)
  avres = True
  avres@gsLineColor = "blue"
  avres@gsLineThicknessF = 1.5
  avres@tfPolyDrawOrder = "PostDraw"
  do r=0 , numregs-1
    dx = 0.15
    xl = (/(xpos(r)-dx),(xpos(r)+dx)/)
    yl = (/     yavg(r),     yavg(r)/) 
    avgln(r) = gsn_add_polyline(wks,plot,xl,yl,avres)
  end do

  resp = True
  resp@gsFillIndex  = 0
  resp@gsFillColor  = "grey"
  resp@tfPolyDrawOrder = "PreDraw"
  xbox = (/1.5,2.,2.,1.5 /)
  ybox = (/.2,.2,0.,0. /)
  pbox = gsn_add_polygon(wks,plot,xbox,ybox,resp)

;  resP = True
;  maximize_output(wks,resP) 
  draw(plot)
  frame(wks)
  
  iform2 = "tif"
  print("saved.") ; "+outim+"."+formo)
  trimPNG(outim+"."+formo)

end