pro gfed4
;
  ; External functions
  FORWARD_FUNCTION CTM_Grid,  CTM_Make_DataInfo, CTM_Type, $
                   MFindFile, Extract_FileName

  ; Input Grid parameters
  InType = CTM_Type( 'generic', res=0.25, halfpolar=0, center180=0 )
  InGrid = CTM_Grid( InType, /No_Vertical )

; write to bpch or ncdf?
bpch = 0
ncdf = 1

; types to write
drymat = 0
daily  = 1
hourly = 1

; path to source files
gfeddir = './hdf5'

;
ftype = ['AGRI', 'BORF', 'DEFO', 'SAVA', 'PEAT', 'TEMF', 'TOTL']
month = ['01/', '02/', '03/', '04/', '05/', '06/', $
         '07/', '08/', '09/', '10/', '11/', '12/']
ig    = [0, 1, 3, 2, 2, 5, 4, 3]
;

for year = 2011,2012 do begin
date=long(year)
dates=strlowcase(strcompress(string(date,format='(a20)'),/remove_all))
gfedfile = gfeddir + '/' + 'GFED4.1s_' + dates + '.hdf5' 
print, gfedfile
fileid=h5f_open(gfedfile)
;

if ( (year eq 1996) || (year eq 2000) || (year eq 2004) || (year eq 2008) || (year eq 2012) ) then begin
 IsLeapYear = 1
endif else begin
 IsLeapYear = 0
endelse

; FOR LEAP YEARS
if IsLeapYear eq 1 then begin
 secsinmonth=[2.6784e6, 2.5056e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.5920e6, 2.6784e6]
; FOR REGULAR YEARS
endif else begin
 secsinmonth=[2.6784e6, 2.4192e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.6784e6, 2.5920e6, 2.6784e6, 2.5920e6, 2.6784e6]
endelse
;
for im = 1, 12 do begin
;
 date=100*long(year)+long(im)
 nmonth = im+1
 nyear = year
 if nmonth eq 13 then begin
    nmonth = 1
    nyear = year+1
 endif
 ndate=100*long(nyear)+long(nmonth)
 YYYYMMDD  = Long( date  ) * 100L + 01L
 YYYYMMDD1  = Long( ndate  ) * 100L + 01L
 Tau0      = Nymd2Tau( YYYYMMDD  )
 Tau1      = Nymd2Tau( YYYYMMDD1 )

; define lon & lat
 lon=(FINDGEN(1440)*0.25)-179.875
 lat=(FINDGEN(720)*0.25)-89.875
 dates=strlowcase(strcompress(string(date,format='(a20)'),/remove_all))

; dry matter (if selected)
 if drymat eq 1 then begin
 OutFilename=strarr(7)
 OutFilename[0]='GFED4_DM_AGRI_' + dates +'.generic.025x025'
 OutFilename[1]='GFED4_DM_BORF_' + dates +'.generic.025x025'
 OutFilename[2]='GFED4_DM_DEFO_' + dates +'.generic.025x025'
 OutFilename[3]='GFED4_DM_SAVA_' + dates +'.generic.025x025'
 OutFilename[4]='GFED4_DM_PEAT_' + dates +'.generic.025x025'
 OutFilename[5]='GFED4_DM_TEMP_' + dates +'.generic.025x025'
 OutFilename[6]='GFED4_DM_TOTL_' + dates +'.generic.025x025'
;
 OutVar=strarr(7)
 OutVar[0]='DM_AGRI'
 OutVar[1]='DM_BORF'
 OutVar[2]='DM_DEFO'
 OutVar[3]='DM_SAVA'
 OutVar[4]='DM_PEAT'
 OutVar[5]='DM_TEMP'
 OutVar[6]='DM_TOTL'
;
 dataall=fltarr(1440,720)
 datafxg=fltarr(1440,720,8)
;
;
 for itype=0,6 do begin
;
  igx=ig[itype]
;
  if itype eq 0 then begin
   dataset = '/emissions/' + month[im-1] + '/DM'
  endif else begin
   dataset = '/emissions/' + month[im-1] +  '/partitioning/DM_' + ftype[itype-1]
  endelse

  dataid=h5d_open(fileid, dataset)
;
  ; Read data (g DM/m^2/month for TOTL or fractions for partitions)
  data=h5d_read(dataid)
;
  dataf=fltarr(1440,720)
;
  ; Flip j so that j=1 is at south pole
;
  for j=0,719 do begin
    dataf[*,719-j] = data[*,j]
  endfor
;
  if itype eq 0 then begin
   dataall[0:1439, 0:719] = dataf[0:1439, 0:719]
   dataf[0:1439, 0:719] = 1.0
  endif 
;
  datafx=fltarr(1440,720)
;
  ; Convert all data to g DM/m^2/month
  datafx = dataf * dataall
;
  for j=0,719 do begin
  for i=0,1439 do begin
      datafxg[i,j,itype]=datafx[i,j]
  endfor
  endfor
;
 endfor
;
 ; Ensure that total DM burned (datafxg[*,*,7]) is sum of components
 for igx=1,6 do begin
  for j=0,719 do begin
  for i=0,1439 do begin
   datafxg[i,j,7]=datafxg[i,j,7]+datafxg[i,j,igx]
  endfor
  endfor
 endfor
;
 for igx=1,7 do begin

; -- Write to bpch --
  if bpch eq 1 then begin

   ; Convert from kg/m2/month to kg/m2/s
   InData = fltarr(1440, 720)
   ;InData[0:1439,0:719]=datafxg[0:1439,0:719,igx]*1.e-3/secsinmonth[im-1]
   InData[0:1439,0:719]=datafxg[0:1439,0:719,igx]/secsinmonth[im-1]
;
   Tracer = 90L + Long(igx)
 
   ; Make a DATAINFO structure
   Success = CTM_Make_DataInfo( InData,                   $
                                ThisDataInfo,             $
                                ThisFileInfo,             $
                                ModelInfo=InType,         $
                                GridInfo=InGrid,          $
                                DiagN='GFED4-BB',         $
                                Tracer=Tracer,            $
                                Tau0=Tau0,                $
                                Tau1=Tau1,                $
                                Unit='kg/m2/s',           $
                                Dim=[InGrid.IMX,          $
                                     InGrid.JMX, 0, 0],   $
                                First=[1L, 1L, 1L],       $
                                /No_Global )
 
   ; Error check
   if ( not Success ) then Message, 'Cannot make DATAINFO!'
 
   ; Write output file
   CTM_WriteBpch, ThisDataInfo, ThisFileInfo, FileName=OutFileName[igx-1]

   ; Undefine stuff
   UnDefine, ThisDataInfo
   UnDefine, ThisFileInfo
  endif

; -- Write to netCDF --
  if ncdf eq 1 then begin
 
   ; Convert from kg/m2/month to kg/m2/s
   InData = fltarr(1440, 720, 1)
   ;InData[0:1439,0:719,0]=datafxg[0:1439,0:719,igx]*1.e-3/secsinmonth[im-1]
   InData[0:1439,0:719,0]=datafxg[0:1439,0:719,igx]/secsinmonth[im-1]

   ThisFileName = './ncdf/' + OutFilename[igx-1] + '.nc'
   ThisVar      = OutVar[igx-1]

   ; open/create file
   id = NCDF_CREATE(ThisFileName, /CLOBBER)

   ; begin definition section
   NCDF_CONTROL, id, /FILL

   xid = NCDF_DIMDEF(id, 'lon', 1440)
   yid = NCDF_DIMDEF(id, 'lat',  720)
   tid = NCDF_DIMDEF(id, 'time', /UNLIMITED)

   tmid = NCDF_VARDEF(id, 'time', [tid], /FLOAT)
   NCDF_ATTPUT, id, tmid, 'units', 'hours since 1985-01-01 00:00:00'

   lnid = NCDF_VARDEF(id, 'lon', [xid], /FLOAT)
   NCDF_ATTPUT, id, lnid, 'units', 'degrees_east'

   ltid = NCDF_VARDEF(id, 'lat', [yid], /FLOAT)
   NCDF_ATTPUT, id, ltid, 'units', 'degrees_north'

   did = NCDF_VARDEF(id, ThisVar, [xid,yid,tid], /FLOAT)
   NCDF_ATTPUT, id, did, 'units', 'kg/m2/s'

   ; pass data
   NCDF_CONTROL, id, /ENDEF
   NCDF_VARPUT,  id, lnid, lon
   NCDF_VARPUT,  id, ltid, lat
   NCDF_VARPUT,  id, tmid, Tau0
   NCDF_VARPUT,  id, did,  InData

   ; close
   NCDF_CLOSE, id
  endif
;
 endfor ;type
 endif  ;drymat

; -- Extract daily data --
 if year gt 2002 then begin 

  if daily eq 1 then begin
  daysinmonth = secsinmonth[im-1] / 24.0 / 60.0 / 60.0
  tdim = FIX(daysinmonth)

  ; average fraction per day
  avgfrac = 1.0 / daysinmonth

  dataf=fltarr(1440,720,tdim)
  dataf[0:1439,0:719,0:(tdim-1)]=0.0

  for iday=1,tdim do begin
    dataset = '/emissions/' + month[im-1] + 'daily_fraction/day_' + strtrim(iday,2)
    dataid=h5d_open(fileid, dataset)
    data=h5d_read(dataid)

    ; Flip j so that j=1 is at south pole. Also normalize by avgfrac, so that dataf
    ; contains the emission scale factor relative to a monthly mean value.
    for j=0,719 do begin
      dataf[*,719-j,iday-1] = data[*,j] / avgfrac
    endfor
  endfor

  ; Now that we have all data, set scale factors to 1.0 if there are no scale factors defined at all in a given grid box
  for j=0,719 do begin
  for i=0,1439 do begin
    sum = TOTAL(dataf[i,j,*])
    if sum eq 0.0 then begin 
      dataf[i,j,*]=1.0
    endif
  endfor
  endfor

  ; Define time array
  times=fltarr(tdim)
  times[0] = Tau0
  for j=1,(tdim-1) do begin
    times[j] = times[j-1] + 24.0
  endfor

   ThisFileName = './ncdf/GFED4_dailyfrac_gen.025x025.' + dates + '.nc'
   ThisVar      = "GFED_FRACDAY"

   ; open/create file
   id = NCDF_CREATE(ThisFileName, /CLOBBER)

   ; begin definition section
   NCDF_CONTROL, id, /FILL

   xid = NCDF_DIMDEF(id, 'lon', 1440)
   yid = NCDF_DIMDEF(id, 'lat',  720)
   tid = NCDF_DIMDEF(id, 'time', /UNLIMITED)

   tmid = NCDF_VARDEF(id, 'time', [tid], /FLOAT)
   NCDF_ATTPUT, id, tmid, 'units', 'hours since 1985-01-01 00:00:00'

   lnid = NCDF_VARDEF(id, 'lon', [xid], /FLOAT)
   NCDF_ATTPUT, id, lnid, 'units', 'degrees_east'

   ltid = NCDF_VARDEF(id, 'lat', [yid], /FLOAT)
   NCDF_ATTPUT, id, ltid, 'units', 'degrees_north'

   did = NCDF_VARDEF(id, ThisVar, [xid,yid,tid], /FLOAT)
   NCDF_ATTPUT, id, did, 'units', '1'

   ; pass data
   NCDF_CONTROL, id, /ENDEF
   NCDF_VARPUT,  id, lnid, lon
   NCDF_VARPUT,  id, ltid, lat
   NCDF_VARPUT,  id, tmid, times
   NCDF_VARPUT,  id, did,  dataf

   ; close
   NCDF_CLOSE, id
   endif ;daily

; -- Extract hourly data --
  if hourly eq 1 then begin
  tdim = 8
  dataf=fltarr(1440,720,tdim)
  dataf[0:1439,0:719,0:(tdim-1)]=1.0

  ; average fraction per 3-hour bin
  avgfrac = 1.0 / 8.0

  dc_name = ['UTC_0-3h', 'UTC_3-6h', 'UTC_6-9h', 'UTC_9-12h', 'UTC_12-15h', 'UTC_15-18h', 'UTC_18-21h', 'UTC_21-24h' ]
  for ihr=1,tdim do begin
    dataset = '/emissions/' + month[im-1] + 'diurnal_cycle/' + dc_name[ihr-1] 
    dataid=h5d_open(fileid, dataset)
    data=h5d_read(dataid)

    ; Flip j so that j=1 is at south pole. Also normalize by avgfrac, so that dataf
    ; contains the emission scale factor relative to a daily mean value.
    for j=0,719 do begin
      dataf[*,719-j,ihr-1] = data[*,j] / avgfrac
    endfor
  endfor

  ; Now that we have all data, set scale factors to 1.0 if there are no scale factors defined at all in a given grid box
  for j=0,719 do begin
  for i=0,1439 do begin
    sum = TOTAL(dataf[i,j,*])
    if sum eq 0.0 then begin 
      dataf[i,j,*]=1.0
    endif
  endfor
  endfor

  ; Define time array
  times=fltarr(tdim)
  times[0] = Tau0
  for j=1,(tdim-1) do begin
    times[j] = times[j-1] + 3.0
  endfor

   ThisFileName = './ncdf/GFED4_3hrfrac_gen.025x025.' + dates + '.nc'
   ThisVar      = "GFED_FRAC3HR"

   ; open/create file
   id = NCDF_CREATE(ThisFileName, /CLOBBER)

   ; begin definition section
   NCDF_CONTROL, id, /FILL

   xid = NCDF_DIMDEF(id, 'lon', 1440)
   yid = NCDF_DIMDEF(id, 'lat',  720)
   tid = NCDF_DIMDEF(id, 'time', /UNLIMITED)

   tmid = NCDF_VARDEF(id, 'time', [tid], /FLOAT)
   NCDF_ATTPUT, id, tmid, 'units', 'hours since 1985-01-01 00:00:00'

   lnid = NCDF_VARDEF(id, 'lon', [xid], /FLOAT)
   NCDF_ATTPUT, id, lnid, 'units', 'degrees_east'

   ltid = NCDF_VARDEF(id, 'lat', [yid], /FLOAT)
   NCDF_ATTPUT, id, ltid, 'units', 'degrees_north'

   did = NCDF_VARDEF(id, ThisVar, [xid,yid,tid], /FLOAT)
   NCDF_ATTPUT, id, did, 'units', '1'

   ; pass data
   NCDF_CONTROL, id, /ENDEF
   NCDF_VARPUT,  id, lnid, lon
   NCDF_VARPUT,  id, ltid, lat
   NCDF_VARPUT,  id, tmid, times
   NCDF_VARPUT,  id, did,  dataf

   ; close
   NCDF_CLOSE, id
   endif ;hourly

 endif ; year>=2003

;
endfor ;month
endfor ;year
;
end