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