Improvements of the Driver of ORCHIDEE

These improvements will concern first the specificities of the driver (its capacity to read forcing files), not its performances that will be improved hopefully by the use of the IOServer. The focus is mainly on the use of forcing files at a 6 or 3-hourly time resolution (such as re-analysis), for which the driver interpolates the meteorological fields to calculate half-hourly fields.

Reminder

The driver of ORCHIDEE is expected to read forcing files with a very specific temporal structure. All time information are defined on UTC time. The fields read at a certain time step t0 in the NetCDF file correspond to:

• the instantaneous value at t0+dt of scalar variables : Temperature, Wind Speed, Air Pressure and Specific Humidity.
• the mean value between t0 and t0+dt for fluxes : Long and Short Wave incoming Radiation and the Precipitation.

Where dt is the time step of the forcing file.

AD (26/5/16): Note that, in the above case, the time at which we read the forcing is t0+dt. Let's introduce an index itau of the forcing records: itau=1 is the first record in the nc file for all variables, etc. Imagine we are currently reading the record itau_n, coming juste after the record itau_nm1. These are the notation of dim2driver/read2dim, and we have itau_n=itau_nm1+1. The "time_stamps" corresponding to these two records are separated by dt (the dt between two records, and not dt_sechiba=dt': dt=split*dt'). To make the link with the above notations of time, the time stamp of the record itau_n is t0+dt, and the one of of itau_nm1 is t0.

Now we need to remember where we are in the code, and what we need the forcing for. When we read itau_n, with values corresponding to t0+dt, sechiba hasn't yet run between t0 and t0+dt. It has run (or been initialized) up to t0, with forcing controlled by the forcing records itau_nm1 and before. So when we read the record itau_n corresponding to the time stamp t0+dt, we know the state of sechiba at t0, the forcing values at t0 (from itau_nm1), and we want to correctly define the forcing variables over each dt' between t0 and t0+dt (dt=split*dt').

For me, these are general principles which need to be true whichever the forcing if we decide to read together all the variables that have the same time stamp. The fact that fluxes are averages over a period of length dt after or before the time stamp can be dealt with in the above general framework (but the latter is an opinion).

Interpolation

• Instantaneous fields and the Long-Wave Radiation are then linearly interpolated for each model time-step t' (typically 1/2 hour) based on the values available at t0 and t0+dt with t0 < t' < t0+dt. Should we not change that for long-wave radiation so that we have one consistent procedure for scalars and another for fluxes ?
• The Precipitation field is downscaled to model's temporal resolution using a step distribution function with one parameter (nb_spread, that is the number of half-hourly time steps over which we spread the precipitation). nb_spread can vary from 1 (all precipitations over a forcing time period (3 ou 6 hours) are put on the first model time step) to SPLIT_DT (number of integration time steps with a forcing time period, in that case, the precipitation are distributed uniformly). By default, nb_spread is set to 1. This default value is not appropriate for all climates and it would make sense to have typical values for various regions of the world (see thesis of T. d'Orgeval.).
• The Short-Wave radiation is interpolated using a function distribution that corresponds to the solar angle distribution over a forcing time period and conserving the average flux for the [t0,t0+dt] interval. Under low incidence angles this interpolation gives at times strange values, this needs to be checked and improved.

AD (26/5/2016): I agree on the above general guidelines for interpolation, but :

• LW should not be interpolated linearly: the mean value over the appropriate interval ([t0,t0+dt] in dim2driver) is the value read at t0+dt ; if we linearly interpolate LW between t0 and t0+dt, the mean value becomes LW(t0)+LW(t0+dt), which is usually different from LW(t0+dt). So if we do a linear interpolation, we change the amount of energy that is received by the surface compared to what is defined in the forcing file.
• Apparently, the default behavior in dim2driver is to linearly interpolate LW (using the same interpolation scheme used for the instantaneous variables like Tair). This results from inter_lin=F by default, which induces netrad_cons=F and leads to use a linear interpolation, in contradiction with inter_lin=F. This needs to be double checked by someone, and reported in a ticket if confirmed.
• The simplest conservative interpolation is to use a uniform function, but it creates "steps". If we wanted to smooth this and remain conservative, we could use a function to constrain the time variations, as it is done with the zenith angle for SW. Based on the physics of LW and the data we have, we could define T4_mean as the mean of Tair⁴ over the 6 (split) dt' in the interpolation period. Then the value of LW(t0+n*dt') = LW(itau_n) * Tair(t0+n*dt')⁴ / T4_mean. I think this is conservative, but it needs to checked.
• Shall we understand that the default spred_prec is one in orchidee_driver ? We recently made a collective choice to take spred_prec=split/2, is it overlooked in orchideedriver?
• For the SW, what are the low incidence angle problems? It must be kept in mind that there is no particular reason why the diurnal cycle of SW should be either very smooth or well centered at noon. We do not deal with extra-terrestrial radiation, but incident SW at the surface, and the noise induced by clouds can be very large.
• Finally, the graphics below remain mysterious to me and I can't see any clear link with the above principles.

Jan's understanding of the working of dim2driver (Jan 26/05/2016)

It is important to understand that readim2.f90 is written in time steps and not physical time. There are 2 time step counters in the driver :

• itau_n : is the last read index of the time axis in the forcing file and readdim2.f90. itau_nm1 is the previously read value and will be used with itau_n for the interpolation. Note that itau_nm1 = itau_n - 1.
• itauin : is the time stepping of ORCHIDEE and takes 'split' values between itau_nm1' and itau_n.

The dates at each of the time steps is given by t(itau_n) or t(itauin) for ORCHIDEE and we can write t(itau_nm1)=t0 et t(itau_n)=t0+dt.

The readdim2.f90 will then interpolate using values at itau_nm1 and itau_n to obtain the forcing at itauin, which is the current time step of ORCHIDEE. Once itauin has stepped through split values, the forcing values of itau_n will be copied to itau_nm1 and the new values for itau_n will be read from the netCDF file. Because of this time-step approach the physical time chosen at itau_n=1 is very important.

The time-stepping is done independently from the physical time. Thus it is the choice of the physical time at itau_n=1 which will determine at which dates the values of the forcing at itau_nm1 and itau_n are valid. This poses no problem for instantaneous variables in the forcing file as then the indexing space (itau_n) and the physical time are identical.

The issue becomes more complicated when we consider fluxes which are averaged over an interval. The only valid assumption is that the dates corresponding to index itau_n is within the time interval over which the average was performed. In readim2.f90 there is no information if the date of n is at the start, end or anywhere else within the averaging period.

The graphic below explains probably better the above description.

The upper panel of the figure is meant to illustrate that the interpolation for fluxes in all generality (here for fluxes valid over the forcing time step before the instantaneous values) will not produce an interpolated value at the same time as the scalar field. readdim2.f90 assumes implicitly that the interval of validity of the fluxes at index itau_n are between itau_n-1/2 and itau_n+1/2 else the flux and scalar interpolations give results at different physical times.

The following explanations will probably help understand the meaning of the interpolated values in physical time. Let us suppose itauin is exactly the middle of the interval itau_nm1 and itau_n. From the upper panel of the figure we can deduce the following :

• Tair(itauin) = (Tair(itau_nm1)+Tair(itau_n))/2 will be valid at t(itauin) = (t(itau_n-1)+t(itau_n))/2 thus the interpolated value is exactly in the middle of the physical time interval.
• LWdown(itauin) = (LWdown(itau_nm1)+LWdown(itau_n))/2 will be valid at t(itauin) = (t(itau_n-2)+t(itau_n))/2 = t(itau_nm1). This is thus different from the time at which Tair(itauin) is valid. Exactly half a forcing time step away.

The lower panel will show that this error does not occur only if, the fluxes are averages centred on the time corresponding at which the scalars are valid. Thus the old driver is only correct if the fluxes at time step itau_n are exactly valid over the time interval t(itau_n-1/2) and t(itau_n+1/2).

AD(27/05/2016): I disagree with the above statement that motivates the development of a new interpolation method : "The time-stepping is done independently from the physical time. [...] The issue becomes more complicated when we consider fluxes which are averaged over an interval. The only valid assumption is that the dates corresponding to index itau_n is within the time interval over which the average was performed. In readim2.f90 there is no information if the date of n is at the start, end or anywhere else within the averaging period."

First, the time is defined without ambiguity from the initial time and the nb of records that were read. Second, when we read a flux at the record itau_n, we perfectly know to which interval it corresponds. For WFDEI for instance, we know from ​http://www.eu-watch.org/gfx_content/documents/README-WFDEI%20(v2016).pdf that the record at itau_n holds the average flux over the 3-hourly period that precedes the "time stamp" of itau_n (t0+dt according to the above notations). This is perfectly clear, perfectly coherent with the first panel of the above graphic, and tractable by dim2driver with the good values of inter_lin and netrad_cons.

It further means that, whichever the interpolation we perform, for WFDEI case, the average of the interpolated flux over the 3h that precede t0+dt must equal the value read in record itau_n. It thus implies that any interpolation that does not preserve the mean over [t0,t0+dt] is inappropriate for such a flux.

Looking at LWdown at the bottom of the 2nd panel above, the mean between [itau_nm1,itau_n]=[t0,t0+dt] is not equal to the value read at itau_n in any of these two cases: (a) uniform values over the green and brown intervals, (b) linear interpolation betwen the centre of the green and the centre of the brown intervals. Maybe another solution is coded, but then, it would be nice to explain it.

Add an explicit time information within the forcing files

As described above, the driver of ORCHIDEE is relatively rigid because it makes strict assumptions on temporal specification for the forcing files. At first, we do not expect to develop a more flexible driver but we would like to add explicit time information in the driver in order to avoid misuse. We suggest to add a time_bounds attribute to the time variables stored in the Netcdf forcing file according to the CM convention : ​http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html

Each file will contain at least 2 time axes, one for the scalar variables and one for the flux variables. Each of the time axes will have its "time_bounds" variable in order to describe the interval each value describes. A time_op attribute (average, point) for each variable could be also added. The time specification of a forcing file (based on time_bounds and time_op attributes) should be read by the driver, in order to check for consistency.

Questions that remain:

• How to implement in detail such a solution without perturbing too much the day-to-day usage of ORCHIDEE.
• It would be most easily managed if we build a library of forcing files available to all at IPSL and which contains only forcing data sets with all the information.
• Should we try to support old and incomplete forcing files and for how long ?
• The improved driver should use all the extra time information in the forcing files but this might require some changes to IOIPSL.

Restart the last forcing time step read

As the driver is built now, the first time step read (for the Temperature for instance) is valid for t+dt. This means that for a yearly file at a 6-hourly time resolution, the first temperature read is for "1-JAN 06:00".
For a simulation starting from scratch (no restart), the driver uses the first value for the second day ("2-JAN 00:00") in order to interpolate with the first one readed (and to calculate half-hourly values between "1-JAN 00:00" and "1-JAN 06:00". At the end of a simulation, the last forcing file read is not written in the restart file of the driver.
When we perform monthly runs using a yearly forcing file, it is managed within the driver to use the appropriate time step in the forcing file but when we perform yearly runs using a yearly file, we process as for a run starting from scratch (while we could use the last values read during the last run).
We suggest to restart the last time step read of the meteorological variables in order to be make possible this.

Once the improved driver knows about the time interval over which the forcing is valid, then we can test if the restart date is correct and we can dump the next restart at the end of the time interval.

Remove the reading of the watchout files

In ORCHIDEE, there is the possibility to write in an output file, the meteorological fields used during the simulation. It is a feature useful when doing a coupled simulation (ORCHIDEE coupled to LMDz) in order to be capable to re-use this 'forcing' in off-line mode. These output files are named WATCHOUT files. These WATCHOUT files are not based on the same template than the standard forcing files read by ORCHIDEE.
Consequently, there are specific reading for these WATCHOUT files.
We suggest to :

• either generate the WATCHOUT files in ORCHIDEE under the format than a standard forcing file
• or, develop a specific tool that convert the WATCHOUT file into the same format than a standard forcing file. This tool will be launch as a pre-treatment.