Page History

...

• 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_perioddiff
• 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)extract a sub-area from the computed precipitation data
• compute the mean value over that a sub-area for each field in precip_diff (use the averageintegrate() 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

...

You can get the date and time of a field as numbers like this:

field_valid_dated = grib_get_stringlong(precip_3hdiff[i], 'validityDate')
field_valid_timet = grib_get_stringlong(precip_3hdiff[i], 'validityTime')

then combine those numbers into proper date variables.

...

dates = nil
for i = .... do
    ddt = .....  # construct a date/time variable
    dates = dates & [ddt]
end for

The integrate() function can be called in a number of ways, but we will use it like thisTo extract a sub-area, use the following syntax (derived from using the GRIB Filter icon's Area  parameter):

 f means = read(data: fintegrate(precip_diffs, area:[N,W,S,E])

The average() function returns to compute the mean values over a sub-area of each field. The result is a list of values, a mean for each field. You can directly multiply a list variable by a number to obtain a new list where each element has been multiplied.

The final calculation requires converting the time intervals into hours (because if the time difference between two steps is 7 hours, then the rate of precip per hour is the mean precip value divided by 7).

 XXXXXXXX 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_perioddiff = 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_valid_date  = grib_get_stringlong(precip_3hdiff[i], 'validityDate')
    t  = grib_get_stringlong(precip_3hdiff[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]
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)