load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  Model Ouputs ;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  DATADir = "./"
  FILES = systemfunc (" ls -1 " + DATADir + "wrfout_* ")
  numFILES = dimsizes(FILES)
  print("numFILES = " + numFILES)

  ifil = 0
  a = addfile(FILES(ifil)+".nc","r")

  lat = wrf_user_getvar(a,"XLAT",0)
  lon = wrf_user_getvar(a,"XLONG",0)

 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 ;;;Here we calculate the desired variables for every time available do loop;;;;;;
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

      it = 0 
      FirstTime = True

     ;========
     aJ = addfile(FILES(0)+".nc","r")
     bJ = addfile(FILES(1)+".nc","r")
     ;========
     
     pres = wrf_user_getvar(a,"pres",it) 
     tk = wrf_user_getvar(a,"tk",it) 

     Rspe=287.058
     density=pres/(Rspe*tk)

     lat    = wrf_user_getvar(a,"XLAT",it)
     lon    = wrf_user_getvar(a,"XLONG",it)  

     ;======
     dustJ1  = (/wrf_user_getvar(aJ,"DUST_1",it), wrf_user_getvar(bJ,"DUST_1",it)/) 
     dustJ2  = (/wrf_user_getvar(aJ,"DUST_2",it), wrf_user_getvar(bJ,"DUST_2",it)/) 
     dustJ3  = (/wrf_user_getvar(aJ,"DUST_3",it), wrf_user_getvar(bJ,"DUST_3",it)/) 
     dustJ4  = (/wrf_user_getvar(aJ,"DUST_4",it), wrf_user_getvar(bJ,"DUST_4",it)/) 
     dustJ5  = (/wrf_user_getvar(aJ,"DUST_5",it), wrf_user_getvar(bJ,"DUST_5",it)/) 

     dust_tot=(/dustJ1,dustJ2,dustJ3,dustJ4,dustJ5/)

     dustnew = dim_avg_n(dust_tot, (/0,1,2/)) 

     dimx = dimsizes(dustnew) 
     nlat = dimx(0)  
     mlon = dimx(1)
     npts = nlat*mlon

     if (FirstTime) then           ;NCL
     fname = "Dust_NEW.txt"
     data = new(npts,"string")
     end if

     npt=-1

     do nl =0,nlat-1
      do ml=0,mlon-1
       npt = npt+1
      
       data(npt) = sprinti("%0.5i ",(npt+1))
       data(npt) = data(npt) + sprintf("%7.3f ", lat(nl,ml))
       data(npt) = data(npt) + sprintf("%7.3f ", lon(nl,ml))

       data(npt) = data(npt) + sprintf("%7.5f ", dustnew(nl,ml))

      end do              ; do-loop lon
   end do                 ; do-loop lat

  asciiwrite (fname, data)