;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
;---Read data 
    dir    = "./"
   fnames = systemfunc("ls " + dir + "AOT_PM_2004-2011.nc")
    a      = addfile(fnames,"r")

  totc = a->pm_kunal

  lat  = a->lat_2
  lon  = a->lon_2

 ; printVarSummary(totc)

;---Open workstation 
  wks_type = "png"                       
  wks = gsn_open_wks(wks_type,"iiserb_campus_OND_CoD_PM")


  ;---Set a lat/lon box for extracting a subset.
  latv = (/18.3, 28.3/)
  lonv = (/72.1, 82.1/)

;---Calculate the indexes for this lat/lon box and use to extract data
  nm   = getind_latlon2d (lat, lon, latv, lonv)
  ilt1 = nm(1,0)  ; start lat index
  ilt2 = nm(0,0)  ; end lat index
  iln1 = nm(0,1)  ; start lon index
  iln2 = nm(1,1)  ; end lon index
;  print("ilt1 / ilt2 " + ilt1 + "/" + ilt2)
;  print("iln1 / iln2 " + iln1 + "/" + iln2)

  data_subset        = totc(:,ilt1:ilt2,iln1:iln2)
 ; printVarSummary(data_subset)

 ; print(data_subset)

 ;---Set a lat/lon box for extracting a subset.
  latz = (/23.285/)
  lonz = (/77.276/)

  mn = getind_latlon2d (lat,lon, latz, lonz)
;  printVarSummary(mn)
;  print(mn)

  do k=0,dimsizes(latz)-1
     n = mn(k,0)
     m = mn(k,1)
     print("nearest location: k="+k+"  n="+n+" m="+m+"  "+lat(n,m)+"  "+lon(n,m))
   ;  print("               totc="+totc(:,n,m))   ; nearest grid point at all time steps
     print("-----") 
  end do

   data_single = totc(:,n,m)
; print(data_single)
; printVarSummary(data_single)
   print("-----") 

 ;-----calculation of coefficient of divergence
    DATA_SINGLE_3D = conform(data_subset, data_single, 0)  
   x1 = (DATA_SINGLE_3D - data_subset)
   x2 = (DATA_SINGLE_3D + data_subset)
   CoD1   = (x1/where(x2.ne.0,x2,x2@_FillValue))^2 
  ; printVarSummary(CoD1)
  ; printMinMax(CoD1, 0)
   print("---")

   CoD2 = dim_avg_n_Wrap(CoD1, 0)    ; (nlat,mlon)
   CoD2@long_name = "Coef of Divergence"
  ; printVarSummary(CoD2)
  ; printMinMax(CoD2,0)

 ;  CoD2_avg = avg(CoD2)                        ; small area; no need to weight
                                               ; for plotting only
   
   CoD_final = sqrt(CoD2)
   printVarSummary(CoD_final)
   printMinMax(CoD_final,0)
   
   xAvg = SUM[CoD_final*w]/SUM[w]
   
   
end