energy/05.monthly_avg.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:74889f3e69f5
       
     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 = "Vielsalm"
       
    16   file_name_add = "carboeurope"
       
    17   year_start = 1998
       
    18   year_end   = 2003
       
    19   year_leap1 = 1996
       
    20   year_leap2 = 2000
       
    21 
       
    22   lat_out = 50.3089
       
    23   lon_out = 5.99861
       
    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 
       
    49   do m = 0,nyear-1
       
    50      year = year_start + m
       
    51 
       
    52 ;    if (mod(year,4) .eq. 0) then
       
    53      if (year .eq. year_leap1 .or. year .eq. year_leap2 ) then         
       
    54         hour_of_month = end_of_month2 * 24.
       
    55      else        
       
    56         hour_of_month = end_of_month1 * 24.
       
    57      end if
       
    58 ; print (hour_of_month) 
       
    59      
       
    60 ;    input file
       
    61      fili = file_name_add+"."+station_name+"."+year+".nc"
       
    62      b = addfile(diri+fili,"r")
       
    63      print (fili)
       
    64 
       
    65 ;    co2 flux unit: umol m-2 s-1
       
    66      data1 = b->FCO2M
       
    67 ;    net radiation flux unit: W m-2 
       
    68      data2 = b->RN
       
    69 ;    sensible heat flux unit: W m-2 
       
    70      data3 = b->H
       
    71 ;    latent   heat flux unit: W m-2 
       
    72      data4 = b->LE
       
    73 
       
    74 ;    time unit: hour since 01-01 (mm-dd)
       
    75      time = b->time
       
    76 
       
    77   do n= 0,nmon-1
       
    78      if (n.eq.0) then
       
    79         timeL = 0.
       
    80         timeR = hour_of_month(0)
       
    81      else
       
    82         timeL = hour_of_month(n-1)
       
    83         timeR = hour_of_month(n)
       
    84      end if
       
    85 ;    print (timeL)
       
    86 ;    print (timeR)
       
    87      
       
    88      i = ind(time.ge.timeL .and. time.lt.timeR)
       
    89 ;    print (i)
       
    90 
       
    91      data_out1(m,n) = avg(data1(i))
       
    92      data_out2(m,n) = avg(data2(i))
       
    93      data_out3(m,n) = avg(data3(i))
       
    94      data_out4(m,n) = avg(data4(i))
       
    95 
       
    96      time_out(m,n) = year*100 + n + 1
       
    97 ;    print (time_out(m,n))
       
    98 
       
    99      delete (i)
       
   100   end do
       
   101      if (m.lt.nyear-1) then  
       
   102         delete (data1)
       
   103         delete (data2)
       
   104         delete (data3)
       
   105         delete (data4)
       
   106         delete (time)
       
   107      end if
       
   108   end do
       
   109   
       
   110   data_out1!0  = "year"
       
   111   data_out1!1  = "month"
       
   112   data_out1@long_name  = data1@long_name
       
   113   data_out1@units      = data1@units
       
   114   data_out1@_FillValue = 1.e+36
       
   115 ; print (data_out1)
       
   116 
       
   117   data_out2!0  = "year"
       
   118   data_out2!1  = "month"
       
   119   data_out2@long_name  = data2@long_name
       
   120   data_out2@units      = data2@units
       
   121   data_out2@_FillValue = 1.e+36
       
   122 ; print (data_out2)
       
   123 
       
   124   data_out3!0  = "year"
       
   125   data_out3!1  = "month"
       
   126   data_out3@long_name  = data3@long_name
       
   127   data_out3@units      = data3@units
       
   128   data_out3@_FillValue = 1.e+36
       
   129 ; print (data_out3)
       
   130 
       
   131   data_out4!0  = "year"
       
   132   data_out4!1  = "month"
       
   133   data_out4@long_name  = data4@long_name
       
   134   data_out4@units      = data4@units
       
   135   data_out4@_FillValue = 1.e+36
       
   136 ; print (data_out4)
       
   137 
       
   138   time_out!0  = "year"
       
   139   time_out!1  = "month"
       
   140   time_out@long_name  = "current date as yyyymm"
       
   141   time_out@units      = "current date as yyyymm"
       
   142   print (time_out)
       
   143                  
       
   144   lat_out!0  = "lat"
       
   145   lon_out!0  = "lon"
       
   146   lat_out@units      = "degrees_north"
       
   147   lat_out@long_name  = "Latitude"
       
   148   lon_out@units      = "degrees_east"
       
   149   lon_out@long_name  = "Longitude"
       
   150   print (lat_out)
       
   151   print (lon_out)
       
   152            
       
   153   c->lat      = lat_out
       
   154   c->lon      = lon_out
       
   155   c->date     = time_out
       
   156   c->CO2_FLUX = data_out1
       
   157   c->RAD_FLUX = data_out2
       
   158   c->SH_FLUX  = data_out3
       
   159   c->LH_FLUX  = data_out4
       
   160   
       
   161 end