energy/03.monthly_avg.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:60614259f6aa
       
     1 ;************************************************
       
     2 ;    Read half-hourly data and write monthly data
       
     3 ; input data is  : half-hourly 
       
     4 ;                                        
       
     5 ; output data is :  CO2_flux, RAD_FLUX, SH_FLUX, LH_FLUX
       
     6 ;                   date, lat, lon  
       
     7 ;************************************************
       
     8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
       
     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
       
    10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
       
    11 ;************************************************
       
    12 begin
       
    13 
       
    14   dir_root = "/fis/cgd/cseg/people/jeff/clamp_data/fluxnet/"
       
    15   station_name = "Lethbridge"
       
    16   file_name_add = "ameriflux"
       
    17   year_start = 1999
       
    18   year_end   = 2004
       
    19   year_leap1 = 2000
       
    20   year_leap2 = 2004
       
    21 
       
    22   lat_out = 49.709278
       
    23   lon_out = -112.940167 + 360.
       
    24    
       
    25 ; final output
       
    26   diro = dir_root + station_name + "/"
       
    27   filo = station_name+"_"+year_start+"-"+year_end+"_monthly.nc" 
       
    28   c = addfile(diro+filo,"c")
       
    29   print (filo)
       
    30 
       
    31 ; input dir
       
    32   diri = dir_root + station_name + "/"
       
    33 
       
    34   nyear = year_end - year_start + 1
       
    35   nmon = 12
       
    36   nlat = 1
       
    37   nlon = 1
       
    38 
       
    39 ; day_of_month  = (/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./)
       
    40   end_of_month1 = (/31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334.,365./)
       
    41   end_of_month2 = (/31.,60.,91.,121.,152.,182.,213.,244.,274.,305.,335.,366./)
       
    42 
       
    43   data_out1 = new((/nyear,nmon/),float)
       
    44   data_out2 = new((/nyear,nmon/),float)
       
    45   data_out3 = new((/nyear,nmon/),float)
       
    46   data_out4 = new((/nyear,nmon/),float)
       
    47   time_out  = new((/nyear,nmon/),integer)
       
    48 ; lat_out   = new((/nlat/),float)
       
    49 ; lon_out   = new((/nlon/),float)
       
    50 
       
    51   do m = 0,nyear-1
       
    52      year = year_start + m
       
    53 
       
    54 ;    if (mod(year,4) .eq. 0) then
       
    55      if (year .eq. year_leap1 .or. year .eq. year_leap2 ) then         
       
    56         hour_of_month = end_of_month2 * 24.
       
    57      else        
       
    58         hour_of_month = end_of_month1 * 24.
       
    59      end if
       
    60 ; print (hour_of_month) 
       
    61      
       
    62 ;    input file
       
    63      fili = file_name_add+"."+station_name+"."+year+".nc"
       
    64      b = addfile(diri+fili,"r")
       
    65      print (fili)
       
    66 
       
    67 ;    co2 flux unit: umol m-2 s-1
       
    68      data1 = b->FCO2M
       
    69 ;    net radiation flux unit: W m-2 
       
    70      data2 = b->RN
       
    71 ;    sensible heat flux unit: W m-2 
       
    72      data3 = b->H
       
    73 ;    latent   heat flux unit: W m-2 
       
    74      data4 = b->LE
       
    75 
       
    76 ;    time unit: hour since 01-01 (mm-dd)
       
    77      time = b->time
       
    78 
       
    79   do n= 0,nmon-1
       
    80      if (n.eq.0) then
       
    81         timeL = 0.
       
    82         timeR = hour_of_month(0)
       
    83      else
       
    84         timeL = hour_of_month(n-1)
       
    85         timeR = hour_of_month(n)
       
    86      end if
       
    87 ;    print (timeL)
       
    88 ;    print (timeR)
       
    89      
       
    90      i = ind(time.ge.timeL .and. time.lt.timeR)
       
    91 ;    print (i)
       
    92 
       
    93      data_out1(m,n) = avg(data1(i))
       
    94      data_out2(m,n) = avg(data2(i))
       
    95      data_out3(m,n) = avg(data3(i))
       
    96      data_out4(m,n) = avg(data4(i))
       
    97 
       
    98      time_out(m,n) = year*100 + n + 1
       
    99 ;    print (time_out(m,n))
       
   100 
       
   101      delete (i)
       
   102   end do
       
   103      if (m.lt.nyear-1) then  
       
   104         delete (data1)
       
   105         delete (data2)
       
   106         delete (data3)
       
   107         delete (data4)
       
   108         delete (time)
       
   109      end if
       
   110   end do
       
   111   
       
   112   data_out1!0  = "year"
       
   113   data_out1!1  = "month"
       
   114   data_out1@long_name  = data1@long_name
       
   115   data_out1@units      = data1@units
       
   116   data_out1@_FillValue = 1.e+36
       
   117 ; print (data_out1)
       
   118 
       
   119   data_out2!0  = "year"
       
   120   data_out2!1  = "month"
       
   121   data_out2@long_name  = data2@long_name
       
   122   data_out2@units      = data2@units
       
   123   data_out2@_FillValue = 1.e+36
       
   124 ; print (data_out2)
       
   125 
       
   126   data_out3!0  = "year"
       
   127   data_out3!1  = "month"
       
   128   data_out3@long_name  = data3@long_name
       
   129   data_out3@units      = data3@units
       
   130   data_out3@_FillValue = 1.e+36
       
   131 ; print (data_out3)
       
   132 
       
   133   data_out4!0  = "year"
       
   134   data_out4!1  = "month"
       
   135   data_out4@long_name  = data4@long_name
       
   136   data_out4@units      = data4@units
       
   137   data_out4@_FillValue = 1.e+36
       
   138 ; print (data_out4)
       
   139 
       
   140   time_out!0  = "year"
       
   141   time_out!1  = "month"
       
   142   time_out@long_name  = "current date as yyyymm"
       
   143   time_out@units      = "current date as yyyymm"
       
   144   print (time_out)
       
   145                  
       
   146   lat_out!0  = "lat"
       
   147   lon_out!0  = "lon"
       
   148   lat_out@units      = "degrees_north"
       
   149   lat_out@long_name  = "Latitude"
       
   150   lon_out@units      = "degrees_east"
       
   151   lon_out@long_name  = "Longitude"
       
   152   print (lat_out)
       
   153   print (lon_out)
       
   154            
       
   155   c->lat      = lat_out
       
   156   c->lon      = lon_out
       
   157   c->date     = time_out
       
   158   c->CO2_FLUX = data_out1
       
   159   c->RAD_FLUX = data_out2
       
   160   c->SH_FLUX  = data_out3
       
   161   c->LH_FLUX  = data_out4
       
   162   
       
   163 end