;************************************************
; READ WRF GRID
;************************************************
  inDir            = "./"
  inFile           = "geo_em.d01.markus_mingel.nc"
  inPath           = inDir + inFile
  in               = addfile(inPath,"r")
  WRF_lon          = rm_single_dims(in->XLONG_M)
  WRF_lat          = rm_single_dims(in->XLAT_M)
  printVarSummary(WRF_lat)
  printMinMax(WRF_lat,0)
  print("---")
  printVarSummary(WRF_lon)
  printMinMax(WRF_lon,0)
  print("---")
                   ; use the following for plotting and filtering
  WRF_lon_min = min(WRF_lon)
  WRF_lon_max = max(WRF_lon)
  WRF_lat_min = min(WRF_lat)
  WRF_lat_max = max(WRF_lat)

;************************************************
; READ OMI DATA; EXAMINE DATA DISTRIBUTION
;************************************************
  varname            = "tropospheric_no2_vertical_column"

  inFile             = "QA4ECV_L2_NO2_OMI_20150101T062700_o55658_fitB_v1.nc"
  inPath             = inDir + inFile
  in                 = addfile(inPath,"r")

  OMI_lon            = rm_single_dims(in->longitude)
  OMI_lat            = rm_single_dims(in->latitude)
  OMI_NO2            = rm_single_dims(in->$varname$) 
  printVarSummary(OMI_lat)  
  printMinMax(OMI_lat,0)  
  print("---")
  printVarSummary(OMI_lon) 
  printMinMax(OMI_lon,0)  
  print("---")
  printVarSummary(OMI_NO2) 

  opt_omi  = True
  opt_omi@PrintStat = True
  stat_omi = stat_dispersion(OMI_NO2, opt_omi )
  print("---")

;************************************************
; Remove negative values; Make sure associated lat/lon are also missing  
; Examine distribution after removing negative values
;************************************************
  OMI_NO2  = where(OMI_NO2.le.0, OMI_NO2@_FillValue, OMI_NO2)
  OMI_lat  = where(ismissing(OMI_NO2), OMI_lat@_FillValue, OMI_lat)
  OMI_lon  = where(ismissing(OMI_NO2), OMI_lon@_FillValue, OMI_lon)

  stat_omi = stat_dispersion(OMI_NO2, opt_omi )
  print("---")
  
;************************************************
; Extract indices of OMI *within* the WRF region              
;************************************************
  OMI_NO2_1d = ndtooned(OMI_NO2)
  OMI_lat_1d = ndtooned(OMI_lat)
  OMI_lon_1d = ndtooned(OMI_lon)

  OMI_ind   := ind(.not.ismissing(OMI_NO2_1d) .and. OMI_lat_1d.ge.WRF_lat_min \
                                              .and. OMI_lat_1d.le.WRF_lat_max \
                                              .and. OMI_lon_1d.ge.WRF_lon_min \
                                              .and. OMI_lon_1d.le.WRF_lon_max )
  if (ismissing(OMI_ind(0))) then
      print("No OMI data are with the target WRF grid")
      exit
  end if 

  n_OMI_NO2  = dimsizes(OMI_ind)
  print("n_OMI_NO2="+n_OMI_NO2)
  print("---")

  OMI_lat_1d   :=  OMI_lat_1d(OMI_ind)   ; := is reassignment operator
  OMI_lon_1d   :=  OMI_lon_1d(OMI_ind)
  OMI_NO2_1d   :=  OMI_NO2_1d(OMI_ind)

                                         ; eliminate rightmost swath path
  OMI_NO2_1d    = where(OMI_lon_1d.ge.92, OMI_NO2_1d@_FillValue, OMI_NO2_1d)
  OMI_lat_1d    = where(ismissing(OMI_NO2_1d), OMI_lat_1d@_FillValue, OMI_lat_1d)
  OMI_lon_1d    = where(ismissing(OMI_NO2_1d), OMI_lon_1d@_FillValue, OMI_lon_1d)

  OMI_ind      := ind(.not.ismissing(OMI_NO2_1d))
  OMI_lat_1d   :=  OMI_lat_1d(OMI_ind)   ; := is reassignment operator
  OMI_lon_1d   :=  OMI_lon_1d(OMI_ind)
  OMI_NO2_1d   :=  OMI_NO2_1d(OMI_ind)
 
;************************************************
; ESMF Interpolation
;************************************************
  Opt                   = True
  Opt@SrcFileName       = "SCRIP.nc"
  Opt@DstFileName       = "ESMF.nc"
  Opt@WgtFileName       = "Weight.nc"
  Opt@ForceOverwrite    = True
 ;Opt@InterpMethod      = "neareststod"
  Opt@InterpMethod      = "bilinear"
  Opt@PrintTimings      = True
  Opt@Debug             = False
  Opt@DstGridLat        = WRF_lat
  Opt@DstGridLon        = WRF_lon
  Opt@DstGridType       = "curvilinear"
  Opt@SrcGridLat        = OMI_lat_1d
  Opt@SrcGridLon        = OMI_lon_1d
  Opt@SrcGridType       = "unstructured"

  Opt@SrcRegional       = True
  Opt@DstRegional       = True

  OMI_NO2_regrid     = ESMF_regrid(OMI_NO2_1d, Opt) 

  copy_VarCoords(WRF_lat, OMI_NO2_regrid)
 ;OMI_NO2_regrid@long_name   = OMI_NO2@long_name
  OMI_NO2_regrid@long_name   = "tropospheric NO2"
  OMI_NO2_regrid@units       =  OMI_NO2@units
  OMI_NO2_regrid@coordinates = "XLAT_M XLONG_M"
  OMI_NO2_regrid@coordinates = "XLAT_M XLONG_M"
;;delete( [/OMI_NO2_regrid@lat2d, OMI_NO2_regrid@lon2d/] )  ; use lat2d, lon2d for plotting

  printVarSummary(OMI_NO2_regrid)
  printMinMax(OMI_NO2_regrid,0)
  stat_omi_regrid = stat_dispersion(OMI_NO2_regrid, opt_omi )
  print("---")

;----------------------------------------------------------------------
; Plotting section
;----------------------------------------------------------------------

    scale = 1e14                             
    OMI_NO2_regrid = OMI_NO2_regrid/scale  ; nicer plots
    OMI_NO2_regrid@units       =  OMI_NO2@units+" [1e-14]"
    wks = gsn_open_wks("png","ESMF_regrid_OMI")  ; send graphics to PNG file

;---Resources to share between both plots
    res                      = True              ; Plot mods desired.
    res@gsnMaximize          = True              ; Maximize plot
    res@gsnAddCyclic         = False

    res@cnFillOn             = True              ; color plot desired
    res@cnLinesOn            = False             ; turn off contour lines
    res@cnLineLabelsOn       = False             ; turn off contour labels
    res@cnFillMode           = "RasterFill"      ; turn raster on      

   ;res@cnLevelSelectionMode = "ManualLevels"
   ;res@cnMinLevelValF       =  1.0
   ;res@cnMaxLevelValF       = 16.0
   ;res@cnLevelSpacingF      =  0.5

    res@mpFillOn             = True              ; True is default
    res@mpMinLatF            = WRF_lat_min
    res@mpMaxLatF            = WRF_lat_max
    res@mpMinLonF            = WRF_lon_min
    res@mpMaxLonF            = WRF_lon_max

  ;;res@tiMainString         = ...         

    plot_regrid = gsn_csm_contour_map(wks,OMI_NO2_regrid,res)