...
Note that when passing numeric dates such as 20150105
to other modules, such as the MARS Retrieval module, these do not need to be converted into date variables. However, MARS treats Date and Time as separate parameters, so a date variable would need to be split into these components.
Date arithmetic
When dealing with dates, the number 1 represents one day. So the expression d1 + 1 gives a date one day later than day 1. To compute the difference, in days, between two dates, it's simply:
...
Looping through dates
Three examples (no need to type these in), to get a feel for it:
Code Block | ||
---|---|---|
| ||
for d = 2015-01-01 to 2015-03-01 do print(d) ... # each step is 1 day end for for d = 2015-01-01 to 2015-03-01 by 2 do print(d) ... # each step is 2 days end for for d = 2015-01-01 to 2015-03-01 by hour(6) do print(d) # each step is 6 hours end for |
Computing mean precipitation rate over an area
As an exercise to put all of this together, write a new macro to do the following:
- 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_period
- 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 area for each field (use the
average()
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
Here are some hints to help.
You can get the date and time of a field as numbers like this:
Code Block | ||
---|---|---|
| ||
field_valid_date = grib_get_string(precip_3h[i], 'validityDate')
field_valid_time = grib_get_string(precip_3h[i], 'validityTime') |
then combine those numbers into proper date variables.
A list is built up like this:
Code Block | ||
---|---|---|
| ||
dates = nil for i = .... do d = ..... # construct a date variable dates = dates & [d] end for |
To extract a sub-area, use the following syntax (derived from using the GRIB Filter icon's Area parameter):
Code Block | ||
---|---|---|
| ||
f = read(data: f, area:[N,W,S,E]) |
The average()
function returns 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).
Code Block | ||
---|---|---|
| ||
XXXX remove once tested! precip = read("precip.grib") n = count(precip) # the number of fields in the fieldset precip_period = precip[2, n] - precip[1, n-1] for i = 1 to count(precip) do d_valid_date = grib_get_string(precip_3h[i], 'validityDate') t = grib_get_string(precip_3h[i], 'validityTime') step is 6 hours end for |