1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/energy/01.monthly_avg.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,161 @@
1.4 +;************************************************
1.5 +; Read half-hourly data and write monthly data
1.6 +; input data is : half-hourly
1.7 +;
1.8 +; output data is : CO2_flux, RAD_FLUX, SH_FLUX, LH_FLUX
1.9 +; date, lat, lon
1.10 +;************************************************
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.14 +;************************************************
1.15 +begin
1.16 +
1.17 + dir_root = "/fis/cgd/cseg/people/jeff/clamp_data/fluxnet/"
1.18 + station_name = "Kaamanen"
1.19 + file_name_add = "carboeurope"
1.20 + year_start = 2000
1.21 + year_end = 2003
1.22 +
1.23 +; final output
1.24 + diro = dir_root + station_name + "/"
1.25 + filo = station_name+"_"+year_start+"-"+year_end+"_monthly.nc"
1.26 + c = addfile(diro+filo,"c")
1.27 + print (filo)
1.28 +
1.29 +; input dir
1.30 + diri = dir_root + station_name + "/"
1.31 +
1.32 + nyear = year_end - year_start + 1
1.33 + nmon = 12
1.34 + nlat = 1
1.35 + nlon = 1
1.36 +
1.37 +; day_of_month = (/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./)
1.38 + end_of_month1 = (/31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334.,365./)
1.39 + end_of_month2 = (/31.,60.,91.,121.,152.,182.,213.,244.,274.,305.,335.,366./)
1.40 +
1.41 + data_out1 = new((/nyear,nmon/),float)
1.42 + data_out2 = new((/nyear,nmon/),float)
1.43 + data_out3 = new((/nyear,nmon/),float)
1.44 + data_out4 = new((/nyear,nmon/),float)
1.45 + time_out = new((/nyear,nmon/),integer)
1.46 + lat_out = new((/nlat/),float)
1.47 + lon_out = new((/nlon/),float)
1.48 +
1.49 + do m = 0,nyear-1
1.50 + year = year_start + m
1.51 +
1.52 +; if (mod(year,4) .eq. 0) then
1.53 + if (year .eq. 2000) then
1.54 + hour_of_month = end_of_month2 * 24.
1.55 + else
1.56 + hour_of_month = end_of_month1 * 24.
1.57 + end if
1.58 +; print (hour_of_month)
1.59 +
1.60 +; input file
1.61 + fili = file_name_add+"."+station_name+"."+year+".nc"
1.62 + b = addfile(diri+fili,"r")
1.63 + print (fili)
1.64 +
1.65 +; co2 flux unit: umol m-2 s-1
1.66 + data1 = b->FCO2M
1.67 +; net radiation flux unit: W m-2
1.68 + data2 = b->RN
1.69 +; sensible heat flux unit: W m-2
1.70 + data3 = b->H
1.71 +; latent heat flux unit: W m-2
1.72 + data4 = b->LE
1.73 +
1.74 +; time unit: hour since 01-01 (mm-dd)
1.75 + time = b->time
1.76 +
1.77 + do n= 0,nmon-1
1.78 + if (n.eq.0) then
1.79 + timeL = 0.
1.80 + timeR = hour_of_month(0)
1.81 + else
1.82 + timeL = hour_of_month(n-1)
1.83 + timeR = hour_of_month(n)
1.84 + end if
1.85 +; print (timeL)
1.86 +; print (timeR)
1.87 +
1.88 + i = ind(time.ge.timeL .and. time.lt.timeR)
1.89 +; print (i)
1.90 +
1.91 + data_out1(m,n) = avg(data1(i))
1.92 + data_out2(m,n) = avg(data2(i))
1.93 + data_out3(m,n) = avg(data3(i))
1.94 + data_out4(m,n) = avg(data4(i))
1.95 +
1.96 + time_out(m,n) = year*100 + n + 1
1.97 +; print (time_out(m,n))
1.98 +
1.99 + delete (i)
1.100 + end do
1.101 + if (m.lt.nyear-1) then
1.102 + delete (data1)
1.103 + delete (data2)
1.104 + delete (data3)
1.105 + delete (data4)
1.106 + delete (time)
1.107 + end if
1.108 + end do
1.109 +
1.110 + lat_out = 69.1407
1.111 + lon_out = 27.2950
1.112 +
1.113 + data_out1!0 = "year"
1.114 + data_out1!1 = "month"
1.115 + data_out1@long_name = data1@long_name
1.116 + data_out1@units = data1@units
1.117 + data_out1@_FillValue = 1.e+36
1.118 +; print (data_out1)
1.119 +
1.120 + data_out2!0 = "year"
1.121 + data_out2!1 = "month"
1.122 + data_out2@long_name = data2@long_name
1.123 + data_out2@units = data2@units
1.124 + data_out2@_FillValue = 1.e+36
1.125 +; print (data_out2)
1.126 +
1.127 + data_out3!0 = "year"
1.128 + data_out3!1 = "month"
1.129 + data_out3@long_name = data3@long_name
1.130 + data_out3@units = data3@units
1.131 + data_out3@_FillValue = 1.e+36
1.132 +; print (data_out3)
1.133 +
1.134 + data_out4!0 = "year"
1.135 + data_out4!1 = "month"
1.136 + data_out4@long_name = data4@long_name
1.137 + data_out4@units = data4@units
1.138 + data_out4@_FillValue = 1.e+36
1.139 +; print (data_out4)
1.140 +
1.141 + time_out!0 = "year"
1.142 + time_out!1 = "month"
1.143 + time_out@long_name = "current date as yyyymm"
1.144 + time_out@units = "current date as yyyymm"
1.145 + print (time_out)
1.146 +
1.147 + lat_out!0 = "lat"
1.148 + lon_out!0 = "lon"
1.149 + lat_out@units = "degrees_north"
1.150 + lat_out@long_name = "Latitude"
1.151 + lon_out@units = "degrees_east"
1.152 + lon_out@long_name = "Longitude"
1.153 + print (lat_out)
1.154 + print (lon_out)
1.155 +
1.156 + c->lat = lat_out
1.157 + c->lon = lon_out
1.158 + c->date = time_out
1.159 + c->CO2_FLUX = data_out1
1.160 + c->RAD_FLUX = data_out2
1.161 + c->SH_FLUX = data_out3
1.162 + c->LH_FLUX = data_out4
1.163 +
1.164 +end