diff -r 000000000000 -r 0c6405ab2ff4 energy/08.monthly_avg.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/energy/08.monthly_avg.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,162 @@
+;************************************************
+;    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 = "Hyytiala"
+  file_name_add = "carboeurope"
+  year_start = 1996
+  year_end   = 2003
+  year_leap1 = 1996
+  year_leap2 = 2000
+  year_leap3 = 2004
+
+  lat_out = 61.8474
+  lon_out = 24.2948
+   
+; 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.or.year.eq.year_leap3) 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