Page History

...

As an exercise to put all of this together, we will write a new macro to do the followingcompute the precipitation rate in mm per hour over a sub-area of the globe for each time step. The steps will be:

• compute the 'period' precipitation from precip.grib - this is what we already did earlier, so it's done! Just copy the code to your new macro and change the result variable name to precip_diff
• loop through the fields in the original fieldset and for each:
  • get the date and time of the forecast step
  • combine these into a Metview date variable
  • add it to a list
• use syntax similar to the line of code used to compute the 'period' precipitation to find the differences between the times of adjacent fields (ok, we know it's 3 hours, but in theory it could be anything)
• compute the mean value over a sub-area for each field in precip_diff (use the integrate() function)
• scale up from metres (as the data are stored) to mm by multiplying by 1000
• convert this into a rate, mm per hour, using the time differences computed earlier

...

XXXX remove once tested!

# read the original data file with accumulated precipitation
precip = read("precip.grib")

# compute the 'period' precipitation - the difference between steps
n = count(precip) # the number of fields in the fieldset
precip_diff = precip[2, n] - precip[1, n-1]

# extract the dates and times from the original fields and put into a list
dates = nil
for i = 1 to count(precip) do
    d  = grib_get_long(precip_diff[i], 'validityDate')
    t  = grib_get_long(precip_diff[i], 'validityTime') # ASSUME it's going to be in hours
    dt = date(d) + hour(t/100)
    dates = dates & [dt]
end for
print(dates)
 
# compute the differences between the valid times
date_diffs = dates[2, n] - dates[1, n-1]
print(date_diffs)

means = integrate(precip_diff, [60,-15,50,5]) # [N,W,S,E]
means = means * 1000 # remember to scale up from m to mm!
print(means)

precip_rates_per_hour = means / (date_diffs*24) # date_diffs*24 gives the number of hours in each interval
print(precip_rates_per_hour)

...