;*************************************************
; WRF: DUST_5
;************************************************
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/wrf/WRF_contributed.ncl"

begin

f = addfile ("wrfout_d01_2011-06-22_00:00:00", "r")

Times = f->Times                    ; Times(Time, DateStrLen)    (type character)
Time = wrf_times_c( Times, 0 )      ; "hours since" initial time on file   (double); units recognized by cd_calendar
a = f->TAUAER2
a&Time = Time                       ; associate the 'Time' coordinate with the variable using  standard NCL & syntax
printVarSummary(Time)
print(Time)
    
ymd = cd_calendar(Time, -2)
print(ymd)   

YMD = 20110628                       ; desired date
mt   = ind(YMD.eq.ymd)               ; time index corresponding to the requested data




wks = gsn_open_wks("pdf" ,"WRF_aod_Daily-Average") ; ps,pdf,x11,ncgm,eps 
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map



res = True    ; plot mods desired
res@cnConstFEnableFill =  True
res@gsnMaximize = True
res@gsnSpreadColors = True
res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
WRF_map_c(f, res, 0)
res@tfDoNDCOverlay = True
res@pmTickMarkDisplayMode = "Always"

print("-----")
opt = True

a2 = f->TAUAER2
aDay = calculate_daily_values (a2, "avg", 0, opt)
   printVarSummary(aDay)
   printMinMax (aDay,1)

 b3 = f->TAUAER3
 bDay = calculate_daily_values (b3, "avg", 0, opt)
  printVarSummary(bDay)
     printMinMax (bDay,1)
;----------------------------------------------------------------------
;-----------------------------------------------------------------------
angstrom_exponent = -(log(aDay)-log(bDay))/log(400/600)
 printVarSummary(angstrom_exponent)

AOD550_3D = aDay * ((400/550)^angstrom_exponent)
printVarSummary(AOD550_3D)

AOD550_2D = dim_sum_n_Wrap(AOD550_3D,0)

if (any(isnan_ieee(AOD550_2D))) then
	value = 9.96921e+36
	replace_ieeenan (AOD550_2D,value,0)
	AOD550_2D@_FillValue = value
end if

 printVarSummary(AOD550_2D)
   printMinMax (AOD550_2D,1)

res@tiMainString = "WRF-CHEM (aod_550) " + times(mt)
res@gsnLeftString = a@description

plot = gsn_csm_contour_map(wks,AOD550_2D(:,:),res)

end