;--------------------------------------------------------------------------------
; This is a modified version of wrf_wps_dom which combines code from
; plotgrids.ncl and the original wrf_wps_dom function to add WRF domain
; boxes to an existing map created with wrf_map_overlays.
;
; To use this function, you need to do something like this:
;  . . .
;  contour = wrf_contour(a,wks,ter,res)
;  pltres = True                              ; Set plot options
;  mpres = True                               ; Set map options
;  pltres@PanelPlot = True                    ; Don't remove overlay plot(s) afterwards
;  plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)  ; Plot the data over a map background
;  . . .
;  plot = wrf_add_wps_dom(wks,plot)           ;  Add domain boxes to map
;  draw(plot)     ; This will draw contours, map, domain lines, and text
;  frame(wks)
;----------------------------------------------------------------------
; You can modify the lnres and txres resource settings below to
; customize the look of the domain boxes and text.
;----------------------------------------------------------------------
undef("wrf_add_wps_dom")
function wrf_add_wps_dom(wks[1]:graphic,map[1]:graphic)
begin

;---Need to set this resource so lines added to map are straight and not curved.
  setvalues map
    "mpGreatCircleLinesOn" : True
  end setvalues

; read the following namelist file
  filename = "./namelist.wps"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Do not change anything between the ";;;;;" lines

  dom_opts = True

  dom_opts@max_dom   = stringtoint  (systemfunc("grep max_dom    " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  dom_opts@dx        = stringtofloat(systemfunc("grep dx         " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  dom_opts@dy        = stringtofloat(systemfunc("grep dy         " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  dom_opts@ref_lat   = stringtofloat(systemfunc("grep ref_lat    " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  dom_opts@ref_lon   = stringtofloat(systemfunc("grep ref_lon    " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  test = systemfunc("grep truelat1 " +filename ) 
  if ( .not. ismissing(test) )
    dom_opts@truelat1  = stringtofloat(systemfunc("grep truelat1   " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  else
    dom_opts@truelat1  = 0.0
  end if
  test = systemfunc("grep truelat2 " +filename ) 
  if ( .not. ismissing(test) )
    dom_opts@truelat2  = stringtofloat(systemfunc("grep truelat2   " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  else
    dom_opts@truelat2  = 0.0
  end if
  dom_opts@stand_lon = stringtofloat(systemfunc("grep stand_lon  " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )

  test = systemfunc("grep lambert " +filename ) 
  if ( .not. ismissing(test) )
    dom_opts@map_proj          = "lambert"
  end if
  test = systemfunc("grep mercator " +filename ) 
  if ( .not. ismissing(test) )
    dom_opts@map_proj          = "mercator"
  end if
  test = systemfunc("grep polar " +filename ) 
  if ( .not. ismissing(test) )
    dom_opts@map_proj          = "polar"
  end if
  testa = systemfunc("grep 'lat-lon' " +filename ) 
  if ( .not. ismissing(testa) )
    dom_opts@map_proj          = "lat-lon"
    dom_opts@pole_lat = stringtofloat(systemfunc("grep pole_lat  " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
    dom_opts@pole_lon = stringtofloat(systemfunc("grep pole_lon  " +filename+ " | cut -f2 -d'=' | cut -f1 -d','" ) )
  end if

  parent_id         = new ( dom_opts@max_dom , integer )
  parent_grid_ratio = new ( dom_opts@max_dom , integer )
  i_parent_start    = new ( dom_opts@max_dom , integer )
  j_parent_start    = new ( dom_opts@max_dom , integer )
  e_we              = new ( dom_opts@max_dom , integer )
  e_sn              = new ( dom_opts@max_dom , integer )

  do n = 1, dom_opts@max_dom 

    n0 = n - 1
    parent_id(n0)         = stringtoint(systemfunc("grep parent_id         " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
    parent_grid_ratio(n0) = stringtoint(systemfunc("grep parent_grid_ratio " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
    i_parent_start(n0)    = stringtoint(systemfunc("grep i_parent_start    " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
    j_parent_start(n0)    = stringtoint(systemfunc("grep j_parent_start    " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
    e_we(n0)              = stringtoint(systemfunc("grep e_we              " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
    e_sn(n0)              = stringtoint(systemfunc("grep e_sn              " +filename+ " | cut -f2 -d'=' | cut -f"+n+" -d','" ) )
 
  end do 

  dom_opts@parent_id         = parent_id
  dom_opts@parent_grid_ratio = parent_grid_ratio
  dom_opts@i_parent_start    = i_parent_start
  dom_opts@j_parent_start    = j_parent_start
  dom_opts@e_we              = e_we
  dom_opts@e_sn              = e_sn
  

  ;BPR BEGIN
  ;grid_to_plot = 0 => plot using the corner mass grid points
  ;grid_to_plot = 1 => plot using the edges of the corner grid cells
  ;  This uses the locations of the corner mass grid points +/- 0.5 grid cells
  ;  to get from the mass grid point to the edge of the grid cell
  grid_to_plot = 1
  ;grid_to_plot = dom_opts@grid_to_plot
  if(grid_to_plot.eq.0) then
  ; print("Plotting using corner mass grid points")
  else if(grid_to_plot.eq.1) then
  ; print("Plotting using edges of the corner grid cells")
  else
   print("ERROR: Invalid value for grid_to_plot = "+grid_to_plot)
  end if
  end if
  ;BPR END

  res = True
  res@DX = dom_opts@dx
  res@DY = dom_opts@dy
  res@LATINC = 0.0
  res@LONINC = 0.0
  if ( dom_opts@map_proj .eq. "lambert") then
    dom_opts@map_proj = 1
    res@MAP_PROJ = 1
  end if
  if ( dom_opts@map_proj .eq. "polar") then
    dom_opts@map_proj = 2
    res@MAP_PROJ = 2
  end if
  if ( dom_opts@map_proj .eq. "mercator") then
    dom_opts@map_proj = 3
    res@MAP_PROJ = 3
  end if
  if ( dom_opts@map_proj .eq. "lat-lon") then
    dom_opts@map_proj = 6
    res@MAP_PROJ = 6
    res@LATINC = dom_opts@dy
    res@LONINC = dom_opts@dx
  end if
  res@TRUELAT1 = dom_opts@truelat1
  res@TRUELAT2 = dom_opts@truelat2
  res@STAND_LON = dom_opts@stand_lon

  res@REF_LAT = dom_opts@ref_lat
  res@REF_LON = dom_opts@ref_lon
  if ( isatt(dom_opts,"ref_x") ) then
    res@KNOWNI = dom_opts@ref_x
  else
    res@KNOWNI = int2flt(dom_opts@e_we(0))/2.
  end if
  if ( isatt(dom_opts,"ref_y") ) then
    res@KNOWNJ = dom_opts@ref_y
  else
    res@KNOWNJ = int2flt(dom_opts@e_sn(0))/2.
  end if

  if ( isatt(dom_opts,"pole_lat") ) then
    res@POLE_LAT = dom_opts@pole_lat
  else
    res@POLE_LAT = 90.0
  end if
  if ( isatt(dom_opts,"pole_lon") ) then
    res@POLE_LON = dom_opts@pole_lon
  else
    res@POLE_LON = 0.0
  end if

  ;BPR BEGIN
  ;Determine adjustment needed to convert from mass grid to chosen grid
  if(grid_to_plot.eq.0) then
   adjust_grid = 0.0 
  else if(grid_to_plot.eq.1) then
   adjust_grid = 0.5 
  else 
   print("ERROR: Invalid value for grid_to_plot = "+grid_to_plot)
   adjust_grid = 0.0 
  end if
  end if

  xx = 1.0 - adjust_grid
  yy = 1.0 - adjust_grid
  ;xx = 1.0
  ;yy = 1.0
  ;BPR END
  loc = wrf_ij_to_ll (xx,yy,res)
  start_lon = loc(0)
  start_lat = loc(1)
  ;BPR BEGIN
  ;e_we is the largest U grid point and e_sn the largest V gridpoint
  ;xx = int2flt(dom_opts@e_we(0))
  ;yy = int2flt(dom_opts@e_sn(0))
  ;Change it so it is in terms of mass grid points since wrf_ij_to_ll is
  ;in terms of mass grid points 
  xx = int2flt(dom_opts@e_we(0)-1) + adjust_grid
  yy = int2flt(dom_opts@e_sn(0)-1) + adjust_grid
  ;BPR END
  loc = wrf_ij_to_ll (xx,yy,res)
  end_lon = loc(0)
  end_lat = loc(1)

  dom_opts@start_lat   = start_lat
  dom_opts@start_lon   = start_lon
  dom_opts@end_lat     = end_lat
  dom_opts@end_lon     = end_lon

;---Set line and text resources
  lnres = True 
  lnres@gsLineThicknessF = 4.0   ; was 2.5
  lnres@domLineColors    = (/ "white", "Red" , "Red" , "Blue" /)

  txres = True
  txres@txFont = "helvetica-bold"
  ;txres@txJust = "BottomLeft"
  txres@txJust = "TopLeft"
  txres@txPerimOn = False
  txres@txFontHeightF = 0.015


  if ( dom_opts@max_dom .gt. 1 ) then

    numLineColors = 0
    if ( isatt(lnres,"domLineColors") ) then
       numLineColors = dimsizes(lnres@domLineColors)
    end if

    do idom = 1,dom_opts@max_dom-1

      if ( numLineColors .gt. 0 ) then
        if ( numLineColors .ge. idom ) then
          lnres@gsLineColor = lnres@domLineColors(idom-1)
          txres@txFontColor = lnres@domLineColors(idom-1)
        else
          lnres@gsLineColor = lnres@domLineColors(numLineColors-1)
          txres@txFontColor = lnres@domLineColors(numLineColors-1)
        end if
      end if


      ; nest start and end points in large domain space
      if ( dom_opts@parent_id(idom) .eq. 1) then
        ; corner value
        ;BPR BEGIN
        ;Due to the alignment of nests we need goffset in order to
        ;find the location of (1,1) in the fine domain in coarse domain
        ;coordinates
        ;i_start = dom_opts@i_parent_start(idom)
        ;j_start = dom_opts@j_parent_start(idom)
        goffset = 0.5*(1-(1.0/dom_opts@parent_grid_ratio(idom)))
        i_start = dom_opts@i_parent_start(idom)-goffset
        j_start = dom_opts@j_parent_start(idom)-goffset
        ; end point
        ;Change to mass point
        ;i_end = (dom_opts@e_we(idom)-1)/dom_opts@parent_grid_ratio(idom) + i_start
        ;j_end = (dom_opts@e_sn(idom)-1)/dom_opts@parent_grid_ratio(idom) + j_start
        i_end = (dom_opts@e_we(idom)-2)/(1.0*dom_opts@parent_grid_ratio(idom)) + i_start
        j_end = (dom_opts@e_sn(idom)-2)/(1.0*dom_opts@parent_grid_ratio(idom)) + j_start
        
        if(grid_to_plot.eq.0) then
         adjust_grid = 0.0 
        else if(grid_to_plot.eq.1) then
         adjust_grid = 0.5/(1.0*dom_opts@parent_grid_ratio(idom))
        else 
         print("ERROR: Invalid value for grid_to_plot = "+grid_to_plot)
         adjust_grid = 0.0 
        end if
        end if
        
        ;BPR END
      end if
      if ( dom_opts@parent_id(idom) .ge. 2) then
        ; corner value
        nd = dom_opts@parent_id(idom)
        ;BPR BEGIN
        ;i_points = ((dom_opts@e_we(idom)-1)/dom_opts@parent_grid_ratio(idom))
        ;j_points = ((dom_opts@e_sn(idom)-1)/dom_opts@parent_grid_ratio(idom))
        i_points = ((dom_opts@e_we(idom)-2)/(1.0*dom_opts@parent_grid_ratio(idom)))
        j_points = ((dom_opts@e_sn(idom)-2)/(1.0*dom_opts@parent_grid_ratio(idom)))
        goffset = 0.5*(1-(1.0/(1.0*dom_opts@parent_grid_ratio(idom))))
        ai_start = dom_opts@i_parent_start(idom)*1.0-goffset
        aj_start = dom_opts@j_parent_start(idom)*1.0-goffset
        ;ai_start = dom_opts@i_parent_start(idom)*1.0
        ;aj_start = dom_opts@j_parent_start(idom)*1.0

        if(grid_to_plot.eq.0) then
         adjust_grid = 0.0 
        else if(grid_to_plot.eq.1) then
         adjust_grid = 0.5/(1.0*dom_opts@parent_grid_ratio(idom)) 
        else 
         print("ERROR: Invalid value for grid_to_plot = "+grid_to_plot)
         adjust_grid = 0.0 
        end if
        end if

        do while ( nd .gt. 1)
          ;Note that nd-1 is used in the following because the WPS namelist is
          ;one-based but arrays in NCL are zero-based
          goffset = 0.5*(1-(1.0/(1.0*dom_opts@parent_grid_ratio(nd-1))))
          ;ai_start = ai_start/dom_opts@parent_grid_ratio(nd-1) + dom_opts@i_parent_start(nd-1)
          ;aj_start = aj_start/dom_opts@parent_grid_ratio(nd-1) + dom_opts@j_parent_start(nd-1)
          ai_start = (ai_start-1)/(1.0*dom_opts@parent_grid_ratio(nd-1)) + dom_opts@i_parent_start(nd-1)-goffset
          aj_start = (aj_start-1)/(1.0*dom_opts@parent_grid_ratio(nd-1)) + dom_opts@j_parent_start(nd-1)-goffset
          ;i_points = (i_points/dom_opts@parent_grid_ratio(nd-1))
          ;j_points = (j_points/dom_opts@parent_grid_ratio(nd-1))
          i_points = (i_points/(1.0*dom_opts@parent_grid_ratio(nd-1)))
          j_points = (j_points/(1.0*dom_opts@parent_grid_ratio(nd-1)))
          if(grid_to_plot.eq.0) then
           adjust_grid = 0.0 
          else if(grid_to_plot.eq.1) then
           adjust_grid = adjust_grid/(1.0*dom_opts@parent_grid_ratio(nd-1))
          else 
           print("ERROR: Invalid value for grid_to_plot = "+grid_to_plot)
           adjust_grid = 0.0 
          end if
          end if

          ;nd = nd - 1
          nd = dom_opts@parent_id(nd-1)
        end do
        ;i_start = tointeger(ai_start + .5 )
        ;j_start = tointeger(aj_start + .5 )
        i_start = ai_start
        j_start = aj_start
        ; end point
        ;i_end = i_points + i_start + 1
        ;j_end = j_points + j_start + 1
        i_end = i_points + i_start
        j_end = j_points + j_start
        ;BPR END
      end if

      ; get the four corners

      xx = i_start - adjust_grid
      yy = j_start - adjust_grid
      ;xx = int2flt(i_start)
      ;yy = int2flt(j_start)
      loc = wrf_ij_to_ll (xx,yy,res)
      lon_SW = loc(0)
      lat_SW = loc(1)

      xx = i_end + adjust_grid
      yy = j_start - adjust_grid
      ;xx = int2flt(i_end)
      ;yy = int2flt(j_start)
      loc = wrf_ij_to_ll (xx,yy,res)
      lon_SE = loc(0)
      lat_SE = loc(1)

      xx = i_start - adjust_grid
      yy = j_end + adjust_grid
      ;xx = int2flt(i_start)
      ;yy = int2flt(j_end)
      loc = wrf_ij_to_ll (xx,yy,res)
      lon_NW = loc(0)
      lat_NW = loc(1)

      ;xx = int2flt(i_end)
      ;yy = int2flt(j_end)
      xx = i_end + adjust_grid
      yy = j_end + adjust_grid
      ;BPR END
      loc = wrf_ij_to_ll (xx,yy,res)
      lon_NE = loc(0)
      lat_NE = loc(1)

      xbox = (/lon_SW, lon_SE, lon_NE, lon_NW, lon_SW /)
      ybox = (/lat_SW, lat_SE, lat_NE, lat_NW, lat_SW /)

;;      x_out = new(dimsizes(xbox),typeof(xbox))
;;      y_out = new(dimsizes(ybox),typeof(ybox))
;;      datatondc(map, xbox, ybox, x_out, y_out)
;;      gsn_polyline_ndc(wks, x_out, y_out, lnres)

      lnstr = unique_string("line")
      map@$lnstr$ = gsn_add_polyline(wks, map, xbox, ybox, lnres)

      idd = idom + 1
      dom_text = "d0"+idd
      txstr = unique_string("text")
      if ( txres@txJust .eq. "BottomLeft" ) then
        map@$txstr$ = gsn_add_text(wks,map,dom_text,lon_NW,lat_NW,txres)
;;        gsn_text(wks,map,dom_text,lon_NW,lat_NW,txres)
      else
        map@$txstr$ = gsn_add_text(wks,map,dom_text,xbox(3)+0.5,ybox(3)-0.5,txres)
;;        gsn_text_ndc(wks,dom_text,x_out(3)+0.01,y_out(3)-0.01,txres)
      end if

    end do

  end if

  return(map)

end