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/csm/shea_util.ncl"
begin
;.......................................................................
;This NCL script converts hdf5 data into netCDF for the selected variables
;Note that line numbers 16, 67 must be changed for different variables
;on line 67 the appropriate emission coefficient variable name from lines 23-25 must be changed
;glob_mn variable printed at the end is for total global emissions averaged for the 18 years. To be compared with website tables.
;................................................................................................................................
nyear = 18
st_yr = 1997
tot_mon = 12*18

aeroType    = "SO2"            ;change this for different aerosol type.
outDir      = "/HOME/rmahmood/GFED4s/ncfiles/"   ;........ output dir  
dirName     = "/HOME/rmahmood/GFED4s/"           ;........ input dir
outfileName = aeroType+"_GFED4.0s_1997-2014.nc"

;................ emission factors...............................
;--------- SAVE , BORF,TEMF,  DEFO,PEAT,AGRI
 ef_OC  = (/2.620,9.600,9.600,4.710,6.020,2.300/)
 ef_BC  = (/0.370,0.500,0.500,0.520,0.040,0.750/)
 ef_SO2 = (/0.480,1.100,1.100,0.400,1.100,0.400/)
;.....................................................................
 mnt = (/"01","02","03","04","05","06","07","08","09","10","11","12"/)
 srcs= (/"DM_SAVA","DM_BORF","DM_TEMF","DM_DEFO","DM_PEAT","DM_AGRI"/)

 vr_ini = "/emissions/"
 vr_end = "/partitioning/"
 dm     ="/DM"

  in_tmp = addfile(dirName+"GFED4.0s_"+st_yr+".hdf5","r")
  bm2    = in_tmp->$"/emissions/01/partitioning/DM_AGRI"$
  mon_tot_emi  = new((/tot_mon,720,1440/),typeof(bm2))    ;arrays to hold monthly data

;.................................... for testing output ........
  grid_area = in_tmp->$"/ancill/grid_cell_area"$           ; 
  tst_sum_yr   = new((/18,720,1440/),typeof(bm2))          ;
  tst_sum_src   = new((/12,720,1440/),typeof(bm2))         ;
;................................................................
  delete(bm2)
  delete(in_tmp)

 do ny=0,nyear-1
   year = ny+st_yr
 
   print("                       processing data for the year =   "+year)
   mult = 12*ny

  fileName = "GFED4.0s_"+year+".hdf5"
  in1 = addfile(dirName+fileName,"r")

  vars_all_src = new((/6,720,1440/),float)
 
  do n=0,dimsizes(mnt)-1
  
        var_dum1  =  vr_ini+mnt(n)+dm
        dm_emi    =  in1->$var_dum1$

        do m=0, dimsizes(srcs)-1
      
            var_dum2  = vr_ini+mnt(n)+vr_end+srcs(m)
            vars_prt  = in1->$var_dum2$
          
            vars_all_src(m,:,:) = dm_emi*vars_prt*ef_SO2(m)    ;change ef_* for different species...... calculate emissions for each biomass type.

          end do
    
        mon_tot_emi(n+mult,:,:) = dim_sum_n(vars_all_src,0)              ;sum emissions for each month
        
        tst_sum_src(n,:,:)      = dim_sum_n(vars_all_src,0)              ; for testing output 
   

      end do
       
         tst_sum_yr(ny,:,:)     = dim_sum_n(tst_sum_src,0)               ; for testing output
         

   delete(vars_all_src) 

   end do

;--------------------------- for testing output ---------        
         mean_all_years = dim_avg_n_Wrap(tst_sum_yr,0)
         tgb   = grid_area*mean_all_years
         glob_mn = (sum(tgb))/1e9 
     print(glob_mn)
;.......................................................

;***************************************** coordinate and attribute assignment *******************************  
  
 lat = in1->lat
 lon = in1->lon

    mon_tot_emi = mon_tot_emi(:,::-1,:)    ;............. make it -90 to 90 instead of 90 to -90    

     mon_tot_emi!0    = "time"
     mon_tot_emi!1    = "lat"
     mon_tot_emi!2    = "lon"
     mon_tot_emi&time = ispan(0,tot_mon-1,1)
     mon_tot_emi&lat  = (/lat(::-1,0)/)   ;................... gives 1D lats ........       
     mon_tot_emi&lon  = (/lon(360,:)/)    ;................... gives 1D lons .......

 mon_tot_emi&lat@units = "degrees_north"
 mon_tot_emi&lon@units = "degrees_east"
 mon_tot_emi@units     = "g/m2"

;--------------------------- write netcdf -----------
in31   =addfile(outDir+outfileName,"c")
in31->$aeroType$ = mon_tot_emi
;........................................

end