;************************************************

 ;    Read half-hourly data and write monthly data

 ; input data is  : half-hourly 

 ;                                        

 ; output data is :  CO2_flux, RAD_FLUX, SH_FLUX, LH_FLUX

 ;                   date, lat, lon  

 ;************************************************

 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"   

 ;************************************************

 begin

   dir_root = "/fis/cgd/cseg/people/jeff/clamp_data/fluxnet/"

   station_name = "Vielsalm"

   file_name_add = "carboeurope"

   year_start = 1998

   year_end   = 2003

   year_leap1 = 1996

   year_leap2 = 2000

   lat_out = 50.3089

   lon_out = 5.99861

 ; final output

   diro = dir_root + station_name + "/"

   filo = station_name+"_"+year_start+"-"+year_end+"_monthly.nc" 

   c = addfile(diro+filo,"c")

   print (filo)

 ; input dir

   diri = dir_root + station_name + "/"

   nyear = year_end - year_start + 1

   nmon = 12

   nlat = 1

   nlon = 1

 ; day_of_month  = (/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./)

   end_of_month1 = (/31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334.,365./)

   end_of_month2 = (/31.,60.,91.,121.,152.,182.,213.,244.,274.,305.,335.,366./)

   data_out1 = new((/nyear,nmon/),float)

   data_out2 = new((/nyear,nmon/),float)

   data_out3 = new((/nyear,nmon/),float)

   data_out4 = new((/nyear,nmon/),float)

   time_out  = new((/nyear,nmon/),integer)

   do m = 0,nyear-1

      year = year_start + m

 ;    if (mod(year,4) .eq. 0) then

      if (year .eq. year_leap1 .or. year .eq. year_leap2 ) then         

         hour_of_month = end_of_month2 * 24.

      else        

         hour_of_month = end_of_month1 * 24.

      end if

 ; print (hour_of_month) 

 ;    input file

      fili = file_name_add+"."+station_name+"."+year+".nc"

      b = addfile(diri+fili,"r")

      print (fili)

 ;    co2 flux unit: umol m-2 s-1

      data1 = b->FCO2M

 ;    net radiation flux unit: W m-2 

      data2 = b->RN

 ;    sensible heat flux unit: W m-2 

      data3 = b->H

 ;    latent   heat flux unit: W m-2 

      data4 = b->LE

 ;    time unit: hour since 01-01 (mm-dd)

      time = b->time

   do n= 0,nmon-1

      if (n.eq.0) then

         timeL = 0.

         timeR = hour_of_month(0)

      else

         timeL = hour_of_month(n-1)

         timeR = hour_of_month(n)

      end if

 ;    print (timeL)

 ;    print (timeR)

      i = ind(time.ge.timeL .and. time.lt.timeR)

 ;    print (i)

      data_out1(m,n) = avg(data1(i))

      data_out2(m,n) = avg(data2(i))

      data_out3(m,n) = avg(data3(i))

      data_out4(m,n) = avg(data4(i))

      time_out(m,n) = year*100 + n + 1

 ;    print (time_out(m,n))

      delete (i)

   end do

      if (m.lt.nyear-1) then  

         delete (data1)

         delete (data2)

         delete (data3)

         delete (data4)

         delete (time)

      end if

   end do

   data_out1!0  = "year"

   data_out1!1  = "month"

   data_out1@long_name  = data1@long_name

   data_out1@units      = data1@units

   data_out1@_FillValue = 1.e+36

 ; print (data_out1)

   data_out2!0  = "year"

   data_out2!1  = "month"

   data_out2@long_name  = data2@long_name

   data_out2@units      = data2@units

   data_out2@_FillValue = 1.e+36

 ; print (data_out2)

   data_out3!0  = "year"

   data_out3!1  = "month"

   data_out3@long_name  = data3@long_name

   data_out3@units      = data3@units

   data_out3@_FillValue = 1.e+36

 ; print (data_out3)

   data_out4!0  = "year"

   data_out4!1  = "month"

   data_out4@long_name  = data4@long_name

   data_out4@units      = data4@units

   data_out4@_FillValue = 1.e+36

 ; print (data_out4)

   time_out!0  = "year"

   time_out!1  = "month"

   time_out@long_name  = "current date as yyyymm"

   time_out@units      = "current date as yyyymm"

   print (time_out)

   lat_out!0  = "lat"

   lon_out!0  = "lon"

   lat_out@units      = "degrees_north"

   lat_out@long_name  = "Latitude"

   lon_out@units      = "degrees_east"

   lon_out@long_name  = "Longitude"

   print (lat_out)

   print (lon_out)

   c->lat      = lat_out

   c->lon      = lon_out

   c->date     = time_out

   c->CO2_FLUX = data_out1

   c->RAD_FLUX = data_out2

   c->SH_FLUX  = data_out3

   c->LH_FLUX  = data_out4

 end