var  = "prate"
  yyyy = 1987
  
  diri = "./"
  fili = var+"."+yyyy+".nc"
  pthi = diri+fili
  print(pthi)
  
  f1    = addfiles (pthi, "r")   ; note the "s" of addfile
  ListSetType (f1, "cat")        ; "cat' is fefault
  x     = f1[:]->prate           ; read variable  from all files
  lon2d = f1[0]->lon             ; read NARR longitude from 1st file 
  lat2d = f1[0]->lat             ; read NARR latitude
  
  printVarSummary (x)
  printMinMax(x,0)

  if (var.eq."prate") then 
      x = where(x.lt.0.0, 0.0, x)   
      printMinMax(x,0)
  end if
  print("----------")
  printVarSummary (lat2d)         ; NARR grid
  printMinMax(lat2d,0)
  print("----------")
  printVarSummary (lon2d)
  printMinMax(lon2d,0)
  print("----------")

;---Calculate monthly values from daily values
  
  xMon = calculate_monthly_values (x, "avg", 0,True)
  printVarSummary(xMon)           
  printMinMax(xMon,0)
  print("----------")
  
;---Calculate JJA seasonal values

  xJJA = month_to_season (xMon, "JJA")
  xJJA@long_name = "JJA Seasonal Total"
  printVarSummary(xJJA)
  
;---Change units

  if (var.eq."prate") then       ; kg/m^2 => mm ; kg/m^2/s" ==> mm/s 
      xJJA = xJJA*86400*92       ; (mm/s)(s/day)(day)
      xJJA@units = "mm"          ; Seasonal Total
  end if
  
  printVarSummary(xJJA)
  printMinMax(xJJA,0)
  print("----------")

;---Examine distribution
  opt_sd = True
  opt_sd@PrintStat = True
  statJJA = stat_dispersion(xJJA, opt_sd )
  print("----------")

  nmsgJJA = num(ismissing(xJJA)) ; any missing values?
  print("nmsgJJA="+nmsgJJA)
  print("----------")
  
  yyyymm = cd_calendar(xJJA&time, -1)
;------------------------------------------
;---ESMF regrid
;------------------------------------------
  
;-- define the destination grid
  
  nlat = 190 ;-- number of latitudes
  nlon = 384 ;-- number of longitudes
  nlat@double = True
  nlon@double = True
  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
  lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east")
;;print(lat)
;;print("----------")
;;print(lon)
;;print("----------")

;---Create regrid options

  InterpMethod       = "bilinear"

  Opt                = True
  Opt@SrcTitle       = "NARR grid to Rectilinear: "+InterpMethod   ; optional
  Opt@WgtFileName    = "NARR_to_Rect.WgtFile_"+InterpMethod+"."+nlat+"x"+nlon+".nc"

;---Generate the names of the files containing the source and destination grid descriptions
;---Remove after Part A is complete

    Opt@SrcFileName    = "NARR.SCRIP_grid_description.nc"         ; Name of source and
    Opt@SrcRegional    = True
    
;---If source data contains missing values, set the 
;---special SrcMask2D option to indicate the missing values 

    if (nmsgJJA.gt.0) then
        Opt@SrcMask2D  = where(ismissing(xJJA),0,1)
    end if
    Opt@SrcGridType    = "curvilinear"
    Opt@SrcGridLat     = lat2d
    Opt@SrcGridLon     = lon2d
    Opt@SrcRegional    = True

    DstDirName         = "./"
    Opt@DstFileName    = DstDirName + "Rectilinear.SCRIP_grid_description.nc" 
    Opt@DstGridType    = "rectilinear"
    Opt@DstGridLat     = lat   
    Opt@DstGridLon     = lon
    Opt@DstRegional    = True

;---Specify other options  

    Opt@ForceOverwrite = True
    Opt@InterpMethod   = InterpMethod
    Opt@RemoveSrcFile  = True                  ; remove SCRIP grid destination files
    Opt@RemoveDstFile  = True 
    Opt@NoPETLog       = True                  ; 6.2.1 onward

;---Perform the regrid: NARR ==> rectilinear (_reclin)

    xJJA_regrid = ESMF_regrid(xJJA, Opt)       ; NARR to rectilinear
    printVarSummary(xJJA_regrid)
    printMinMax(xJJA_regrid,0)

;************************************************
; create plots
;************************************************  
    wks = gsn_open_wks("png","NARR_JJA")      ; send graphics to PNG file
    plot = new(2,graphic)                      ; create a plot array
    
    res                      = True
    res@gsnDraw              = False           ; don't draw
    res@gsnFrame             = False           ; don't advance frame
    res@gsnAddCyclic         = False           ; regional data
  
    res@cnInfoLabelOn        = False           ; turn off cn info label
    res@cnFillOn             = True            ; turn on color
   ;res@cnFillMode           = "RasterFill"
    res@cnLinesOn            = False
    res@cnLineLabelsOn       = False
    res@cnLevelSelectionMode = "ExplicitLevels"              
    res@cnLevels             = (/25,50,75,100,150,200,250,300,350,400,450,500,750,1000/)
    res@cnFillPalette        = "WhiteBlueGreenYellowRed"
  ;;res@cnFillPalette        = "precip2_17lev"
  ;;res@cnFillPalette        = (/"Snow","PaleTurquoise","PaleGreen"\ ; 11 contour colors
  ;;                          ,"SeaGreen3" ,"Yellow","Orange"    \
  ;;                          ,"HotPink","Orange","HotPink","Red"\
  ;;                          ,"Violet", "Purple", "Red"       /)       

    res@lbLabelBarOn         = False           ; turn off individual cb's
  
    res@mpMinLatF            = min(lat2d)      ; range to zoom in on
    res@mpMaxLatF            = max(lat2d)
    res@mpMinLonF            = min(lon2d)
    res@mpMaxLonF            = max(lon2d)
    res@mpCenterLonF         = -107.0
    res@mpFillOn             = False
   ;res@mpGridAndLimbOn      = True
   ;res@mpGridLineDashPattern= 10                 ; lat/lon lines as dashed

    res@sfXArray             = lon2d              ; NARR grid
    res@sfYArray             = lat2d  
    
    res@gsnLeftString = ""
    res@gsnRightString= ""
      
    nt   = 0  ; first xJJA
    yyyy = yyyymm(nt)/100

    res@gsnCenterString= "Source NARR"
    plot(0) = gsn_csm_contour_map(wks,xJJA(nt,:,:)     ,res) 

    delete([/res@sfXArray, res@sfYArray/])

    res@gsnCenterString= "NARR: Rectilinear Grid: "+InterpMethod+": "+nlat+"x"+nlon  
    plot(1) = gsn_csm_contour_map(wks,xJJA_regrid(nt,:,:) ,res) 

;************************************************
; create panel
;************************************************
    resP                     = True                ; modify the panel plot
    resP@gsnPanelMainString = "JJA ("+yyyy+"): Total Precipitation (mm)"
    resP@gsnPanelLabelBar    = True                ; add common colorbar
    resP@gsnMaximize         = True  
    resP@pmLabelBarWidthF    = 0.75 
    resP@lbLabelFontHeightF  = 0.0125
    gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot