undef("wrf_wps_dom_mod")
function wrf_wps_dom_mod(wks[1]:graphic,opt_args[1]:logical,lnres[1]:logical,txres[1]:logical) 

begin

  mpres = opt_args
  ;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 = mpres@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 = mpres@dx
  res@DY = mpres@dy
  res@LATINC = 0.0
  res@LONINC = 0.0
  if ( mpres@map_proj .eq. "lambert") then
    mpres@map_proj = 1
    res@MAP_PROJ = 1
  end if
  if ( mpres@map_proj .eq. "polar") then
    mpres@map_proj = 2
    res@MAP_PROJ = 2
  end if
  if ( mpres@map_proj .eq. "mercator") then
    mpres@map_proj = 3
    res@MAP_PROJ = 3
  end if
  if ( mpres@map_proj .eq. "lat-lon") then
    mpres@map_proj = 6
    res@MAP_PROJ = 6
    res@LATINC = mpres@dy
    res@LONINC = mpres@dx
  end if
  res@TRUELAT1 = mpres@truelat1
  res@TRUELAT2 = mpres@truelat2
  res@STAND_LON = mpres@stand_lon

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

  if ( isatt(mpres,"pole_lat") ) then
    res@POLE_LAT = mpres@pole_lat
  else
    res@POLE_LAT = 90.0
  end if
  if ( isatt(mpres,"pole_lon") ) then
    res@POLE_LON = mpres@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(mpres@e_we(0))
  ;yy = int2flt(mpres@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(mpres@e_we(0)-1) + adjust_grid
  yy = int2flt(mpres@e_sn(0)-1) + adjust_grid
  ;BPR END
  loc = wrf_ij_to_ll (xx,yy,res)
  end_lon = loc(0)
  end_lat = loc(1)

  mpres@start_lat   = start_lat
  mpres@start_lon   = start_lon
  mpres@end_lat     = end_lat
  mpres@end_lon     = end_lon
  mpres@mpGreatCircleLinesOn = True

  mp = wrf_wps_map(wks,mpres)

  if ( mpres@max_dom .gt. 1 ) then

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

    do idom = 1,mpres@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 ( mpres@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 = mpres@i_parent_start(idom)
        ;j_start = mpres@j_parent_start(idom)
        goffset = 0.5*(1-(1.0/mpres@parent_grid_ratio(idom)))
        i_start = mpres@i_parent_start(idom)-goffset
        j_start = mpres@j_parent_start(idom)-goffset
        ; end point
        ;Change to mass point
        ;i_end = (mpres@e_we(idom)-1)/mpres@parent_grid_ratio(idom) + i_start
        ;j_end = (mpres@e_sn(idom)-1)/mpres@parent_grid_ratio(idom) + j_start
        i_end = (mpres@e_we(idom)-2)/(1.0*mpres@parent_grid_ratio(idom)) + i_start
        j_end = (mpres@e_sn(idom)-2)/(1.0*mpres@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*mpres@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 ( mpres@parent_id(idom) .ge. 2) then
        ; corner value
        nd = mpres@parent_id(idom)
        ;BPR BEGIN
        ;i_points = ((mpres@e_we(idom)-1)/mpres@parent_grid_ratio(idom))
        ;j_points = ((mpres@e_sn(idom)-1)/mpres@parent_grid_ratio(idom))
        i_points = ((mpres@e_we(idom)-2)/(1.0*mpres@parent_grid_ratio(idom)))
        j_points = ((mpres@e_sn(idom)-2)/(1.0*mpres@parent_grid_ratio(idom)))
        goffset = 0.5*(1-(1.0/(1.0*mpres@parent_grid_ratio(idom))))
        ai_start = mpres@i_parent_start(idom)*1.0-goffset
        aj_start = mpres@j_parent_start(idom)*1.0-goffset
        ;ai_start = mpres@i_parent_start(idom)*1.0
        ;aj_start = mpres@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*mpres@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*mpres@parent_grid_ratio(nd-1))))
          ;ai_start = ai_start/mpres@parent_grid_ratio(nd-1) + mpres@i_parent_start(nd-1)
          ;aj_start = aj_start/mpres@parent_grid_ratio(nd-1) + mpres@j_parent_start(nd-1)
          ai_start = (ai_start-1)/(1.0*mpres@parent_grid_ratio(nd-1)) + mpres@i_parent_start(nd-1)-goffset
          aj_start = (aj_start-1)/(1.0*mpres@parent_grid_ratio(nd-1)) + mpres@j_parent_start(nd-1)-goffset
          ;i_points = (i_points/mpres@parent_grid_ratio(nd-1))
          ;j_points = (j_points/mpres@parent_grid_ratio(nd-1))
          i_points = (i_points/(1.0*mpres@parent_grid_ratio(nd-1)))
          j_points = (j_points/(1.0*mpres@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*mpres@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 = mpres@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(mp, xbox, ybox, x_out, y_out)
;;      gsn_polyline_ndc(wks, x_out, y_out, lnres)

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

      idd = idom + 1
      dom_text = "d0"+idd
      txstr = unique_string("text")
      if ( txres@txJust .eq. "BottomLeft" ) then
        mp@$txstr$ = gsn_add_text(wks,mp,dom_text,lon_NW,lat_NW,txres)
;;        gsn_text(wks,mp,dom_text,lon_NW,lat_NW,txres)
      else
        mp@$txstr$ = gsn_add_text(wks,mp,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(mp)

end