...

Hovmoeller?

Precipitation

Precipitation data provides an interesting challenge. Precipitation fields in MARS are stored as accumulated fields. Visualise the supplied precip.grib icon with the precip_shade visdef. The first field is empty (check using the Cursor Data). The first field has a step of 0, meaning that it contains the total precipitation accumulated between the run time and the run time plus step. Since these are the same, there is no accumulated precipitation! Subsequent steps show more and more precipitation (the amount accumulated over 3, 6, 9, etc hours).

...

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

...

at a point

As an exercise to put all of this together, we will write a new macro to compute the precipitation rate in mm per hour over a sub-area of the globe at a  particular location 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 (which was initialised to nil before the loop)
• 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 extract the mean point value over a sub-area for each field in precip_diff (use the integratenearest_gridpoint() function). Choose a location with some high precipitation
• 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

...

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

The integratenearest_gridpoint() function can be called in a number of ways, but we will use it like this:

 means values = integratenearest_gridpoint(precip_diffs, [N,W,S,E])lat, lon)

to compute the mean values over a sub-area of each field. The result is a list of values, a mean value 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).

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)