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

begin
;************************************************
; define experiment parameters
;************************************************
  EXPNAME   = getenv("expname")
  filepre   = getenv("filepre")
  DATA_dir  = getenv("DATA_DIR")+"/"+EXPNAME+"/"
  FIG_dir   = getenv("FIG_DIR")+"/"
  ix=stringtointeger(getenv("ix"))
  iy=stringtointeger(getenv("iy"))
  ib=stringtointeger(getenv("ib"))
  cfield    = getenv("cfield")
  region    = getenv("region")
  ffpre     = getenv("filefigpre")
  ffpost    = getenv("filefigpost")
;************************************************
; define averaging parameters
;************************************************
  ;yyyymmdd = yyyymmdd_time(2014, 2016, "integer")
  ;yyyymmdd_request=yyyymmdd({20140401:20160630})
  yyyymmdd_request=asciiread("dates_to_process_fig",-1,"string")
  nday=dimsizes(yyyymmdd_request)
  daterange=yyyymmdd_request(0)+"_"+yyyymmdd_request(nday-1)

;************************************************
; read in data
;************************************************
tmp1D = new((/nday,ix*iy/),float)
tmp1D@_FillValue=-999.
do id=0,nday-1
  filename = filepre+yyyymmdd_request(id)+"_fanl"
  if (fileexists(DATA_dir+filename+".txt")) then
     tmp1D(id,:) =  asciiread(DATA_dir+filename+".txt",(/ix*iy/),"float")
  else
     tmp1D(id,:) = tmp1D@_FillValue
  end if
end do
if (cfield.eq."APCP") then
  tmp1D=where(tmp1D.lt.0.,tmp1D@_FillValue,tmp1D)
else
  tmp1D=where(tmp1D.le.0.,tmp1D@_FillValue,tmp1D)
  tmp1D=where(tmp1D.gt.500.,tmp1D@_FillValue,tmp1D)
end if
data1D_fanl = dim_avg_n(tmp1D,0)

tmp1D = new((/nday,ix*iy/),float)
tmp1D@_FillValue=-999.
do id=0,nday-1
  filename = filepre+yyyymmdd_request(id)+"_fcst"
  if (fileexists(DATA_dir+filename+".txt")) then
     tmp1D(id,:) =  asciiread(DATA_dir+filename+".txt",(/ix*iy/),"float") 
  else
     tmp1D(id,:) = tmp1D@_FillValue
  end if
end do
if (cfield.eq."APCP") then
  tmp1D=where(tmp1D.lt.0.,tmp1D@_FillValue,tmp1D)
else
  tmp1D=where(tmp1D.le.0.,tmp1D@_FillValue,tmp1D)
  tmp1D=where(tmp1D.gt.500.,tmp1D@_FillValue,tmp1D)
end if

data1D_fcst = dim_avg_n(tmp1D,0)

tmp1D = new((/nday,ix*iy/),float)
tmp1D@_FillValue=-99.99
do id=0,nday-1
  filename = filepre+yyyymmdd_request(id)+"_rpss"
  if (fileexists(DATA_dir+filename+".txt")) then
     tmp1D(id,:) =  asciiread(DATA_dir+filename+".txt",(/ix*iy/),"float")
  else
     tmp1D(id,:) = tmp1D@_FillValue
  end if
end do
data1D_rpss = dim_avg_n(tmp1D,0)


tmp1D = new((/nday,ix*iy/),float)
tmp1D@_FillValue=-999
do id=0,nday-1
  filename = filepre+yyyymmdd_request(id)+"_mask"
  if (fileexists(DATA_dir+filename+".txt")) then
     tmp1D(id,:) =  asciiread(DATA_dir+filename+".txt",(/ix*iy/),"float")
  else
     tmp1D(id,:) = tmp1D@_FillValue
  end if
end do
data1D_mask = dim_avg_n(tmp1D,0)

printMinMax(data1D_fcst,False)
printMinMax(data1D_fanl,False)


; apply missing value to fanl and fcst
;data1D_fanl=where(ismissing(data1D_rpss),-999.,data1D_fanl)
;data1D_fcst=where(ismissing(data1D_rpss),-999.,data1D_fcst)

printVarSummary(data1D_mask)


;setfileoption ("bin", "ReadByteOrder", "BigEndian")
;data=fbinrecread("/scratch2/NCEPDEV/ensemble/save/Eric.Sinsky/landmask/rfcgrid1deg_shift.bin",0,-1,"float")
;data=fbinrecread("/scratch2/NCEPDEV/ensemble/save/Eric.Sinsky/landmask/rfcgrid1deg.bin",0,-1,"float")
;printVarSummary(data)

;Landmask from Chris's txt file
;data1D_mask=where(data1D_mask.eq.0.0,tmp1D@_FillValue,data1D_mask)





;************************************************
; data 1D --> data 2D
;************************************************
data2D_fanl = new((/iy,ix/),"float")
data2D_fcst = new((/iy,ix/),"float")
data2D_rpss = new((/iy,ix/),"float")
data2D_fanl@_FillValue = -999.
data2D_fcst@_FillValue = -999.
data2D_rpss@_FillValue = -999.
nxy=0
do j = 0, iy-1 
  do i = 0, ix-1
;Yan Luo's landmask binary file
;     if (data(nxy).lt.1.0)then
;Chris's landmask text file
;    if (ismissing(data1D_mask(nxy))) then
    if (ismissing(data1D_rpss(nxy))) then
      data2D_fanl(j,i) = data2D_fanl@_FillValue
      data2D_fcst(j,i) = data2D_fcst@_FillValue 
      data2D_rpss(j,i) = data2D_rpss@_FillValue
    else
      data2D_fanl(j,i) = data1D_fanl(nxy)
      data2D_fcst(j,i) = data1D_fcst(nxy)
      data2D_rpss(j,i) = data1D_rpss(nxy)
    end if
    nxy=nxy+1
  end do
end do

printMinMax(data2D_fcst,False)
printMinMax(data2D_fanl,False)



;exit


;************************************************
; calculate bias
;************************************************
data2D_bias = data2D_fcst-data2D_fanl


;************************************************
; define coordinate variables
;************************************************
lat = fspan(90.,-90.,iy)
if (cfield.eq."APCP") then
  lon = fspan(0.75,359.75,ix)
else
  lon = fspan(0.,357.5,ix)
end if
data2D_bias!0 = "lat";
data2D_bias!1 = "lon";
data2D_bias&lat = lat;
data2D_bias&lon = lon;
data2D_bias&lat@units="degrees_north"
data2D_bias&lon@units="degrees_east"
data2D_bias@_FillValue = -999
copy_VarMeta(data2D_bias,data2D_fanl)
copy_VarMeta(data2D_bias,data2D_fcst)
copy_VarMeta(data2D_bias,data2D_rpss)

;************************************************
; calculate spatial average
;************************************************

rad    = 4.0*atan(1.0)/180.0
clat   = cos(lat*rad)

data2D_rpss_spavg=wgt_areaave_Wrap(data2D_rpss,clat,1.0,0)
data2D_bias_spavg=wgt_areaave_Wrap(data2D_bias,clat,1.0,0)



rpss_spavg_str="RPSS="+sprintf("%5.2f",data2D_rpss_spavg)
bias_spavg_str="Bias="+sprintf("%5.2f",data2D_bias_spavg)


txres               = True                            
txres@txFontHeightF = 20.00             ; Set the font height
xstr=26.0
ystr=289.0



;************************************************
; create plots
;************************************************

;************************************************
  filename_FIG = ffpre+"_rpss_"+ffpost+"_"+daterange
  wks = gsn_open_wks("png",FIG_dir+"/"+filename_FIG)          ; send graphics to PNG file
  
  res                   = True
  res@mpProjection      = "CylindricalEqualArea"       ; choose projection
  res@mpGridAndLimbOn   = True              ; turn on lat/lon lines
  res@mpPerimOn         = False             ; turn off box around plot
  res@mpGridLatSpacingF = 30.               ; spacing for lat lines
  res@mpGridLonSpacingF = 30.               ; spacing for lon lines
  res@mpFillOn          = False
  res@gsnDraw = False
  res@gsnFrame = False
  res@mpGeophysicalLineThicknessF = 2.0
  if (region.eq."_NA") then
    res@mpPerimOn         = True
    res@mpLimitMode       = "LatLon"
    res@mpMinLatF         = 20.
    res@mpMaxLatF         = 60.
    res@mpMinLonF         = -140.
    res@mpMaxLonF         = -50.
    res@mpOutlineBoundarySets     = "geophysicalandusstates"; turn on states
    res@mpDataBaseVersion         = "mediumres"             ; select database
    res@mpDataSetName             = "Earth..2"
  end if
  if (region.eq."_WEST") then
    res@mpPerimOn         = True
    ;res@mpProjection      = "CylindricalEquidistant"
    res@mpLimitMode       = "LatLon"
    res@mpMinLatF         = 30.
    res@mpMaxLatF         = 52.5
    res@mpMinLonF         = -127.5
    res@mpMaxLonF         = -107.5
    
    res@mpOutlineBoundarySets     = "geophysicalandusstates"; turn on states
    res@mpDataBaseVersion         = "mediumres"             ; select database
    res@mpDataSetName             = "Earth..2"
  end if
  if (region.eq."_CONUS") then
    res@mpPerimOn         = True
    ;res@mpProjection      = "CylindricalEquidistant"
    res@mpLimitMode       = "LatLon"
    res@mpMinLatF         = 25.
   res@mpMaxLatF         = 52.5
    res@mpMinLonF         = -125.
    res@mpMaxLonF         = -65.

    res@mpOutlineBoundarySets     = "geophysicalandusstates"; turn on states
    res@mpDataBaseVersion         = "mediumres"             ; select database
    res@mpDataSetName             = "Earth..2"
  end if

  res@gsnAddCyclic = True 
  res@cnFillOn          = True              ; color plot desired
  cmap = read_colormap_file("hkss_cmap.rgb")
  res@cnFillPalette     = cmap              ; set color map
  res@cnLineLabelsOn    = False             ; turn off contour lines
  res@cnLinesOn         = False
  res@txFontHeightF     = 0.015 
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF    = -0.5
  res@cnMaxLevelValF    = 1.0
  res@cnLevelSpacingF   = 0.05

  res@mpShapeMode   = "FreeAspect"           ; make plot bigger
  res@vpWidthF      = 0.8
  res@vpHeightF     = 0.4
  res@vpXF          = 0.1
  res@vpYF          = 0.9

  res@lbLabelFontHeightF  = 0.015           ; label bar font height
  res@tiMainString       = filename_FIG  ; add a title
  res@tiMainFontHeightF  = .018                               ; font height

  contour = gsn_csm_contour_map(wks,data2D_rpss,res)  ; create the plot
  text = gsn_add_text(wks,contour,rpss_spavg_str,ystr,xstr,txres)
  draw(contour)
  frame(wks)
;************************************************

;************************************************
  filename_FIG = ffpre+"_bias_"+ffpost+"_"+daterange
  wks = gsn_open_wks("png",FIG_dir+"/"+filename_FIG)          ; send graphics to PNG file
  res@tiMainString       = filename_FIG

  delete(cmap)
  cmap = read_colormap_file("nrl_sirkes")
  delete(res@cnFillPalette)
  res@cnFillPalette     = cmap
  if (cfield.eq."APCP") then
    res@cnMinLevelValF    = -50
    res@cnMaxLevelValF    = 50
    res@cnLevelSpacingF   = 5
  else
    res@cnMinLevelValF    = -5
    res@cnMaxLevelValF    = 5
    res@cnLevelSpacingF   = 0.5
  end if
  contour = gsn_csm_contour_map(wks,data2D_bias,res)  ; create the plot 
  text = gsn_add_text(wks,contour,bias_spavg_str,ystr,xstr,txres)
  draw(contour)
  frame(wks)

;************************************************

;************************************************
  filename_FIG = ffpre+"_fcst_"+ffpost+"_"+daterange
  wks = gsn_open_wks("png",FIG_dir+"/"+filename_FIG)          ; send graphics to PNG file
  res@tiMainString       = filename_FIG

  delete(cmap)
  delete(res@cnFillPalette)
  if (cfield.eq."APCP") then
    cmap = read_colormap_file("precip4_11lev")
    res@cnMinLevelValF    = 0.
    res@cnMaxLevelValF    = 125.
    res@cnLevelSpacingF   = 12.5
  else
    cmap = read_colormap_file("GMT_haxby")
    res@cnMinLevelValF    = 260.
    res@cnMaxLevelValF    = 320.
    res@cnLevelSpacingF   = 5.
  end if
  res@cnFillPalette     = cmap
  contour = gsn_csm_contour_map(wks,data2D_fcst,res)  ; create the plot
  draw(contour)
  frame(wks)
;************************************************

;************************************************
  filename_FIG = ffpre+"_fanl_"+ffpost+"_"+daterange
  wks = gsn_open_wks("png",FIG_dir+"/"+filename_FIG)          ; send graphics to PNG file
  res@tiMainString       = filename_FIG

  delete(cmap)
  delete(res@cnFillPalette)
  if (cfield.eq."APCP") then
    cmap = read_colormap_file("precip4_11lev")
    res@cnMinLevelValF    = 0.
    res@cnMaxLevelValF    = 125.
    res@cnLevelSpacingF   = 12.5
  else
    cmap = read_colormap_file("GMT_haxby")
    res@cnMinLevelValF    = 260.
    res@cnMaxLevelValF    = 320.
    res@cnLevelSpacingF   = 5.
  end if
  res@cnFillPalette     = cmap
  contour = gsn_csm_contour_map(wks,data2D_fanl,res)  ; create the plot
  draw(contour)
  frame(wks)
;************************************************

end