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

;    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 = "Lethbridge"

  file_name_add = "ameriflux"

  year_start = 1999

  year_end   = 2004

  year_leap1 = 2000

  year_leap2 = 2004

  lat_out = 49.709278

  lon_out = -112.940167 + 360.

; 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)

; lat_out   = new((/nlat/),float)

; lon_out   = new((/nlon/),float)

  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