Even Sphere

It is frequently useful to generate a series of synthetic observations located at roughly evenly-distributed locations [1] on a sphere.

There are three methods described here.

1. A Matlab script and the standard DART observation generation utilities.

2. A csh script with all the parts of 1. (Not available for all models).

3. A stand-alone Fortran program.

The Fortran program does not generate the nice plots that the Matlab process does, but it may be faster and easier to automate for generating a large number of obs.

Matlab Scripts Plus Standard DART Observation Executables

This involves multiple steps:

1. determine how many locations are needed

2. determine the vertical levels needed

3. run the MATLAB function even_sphere.m to create the text file containing the input for create_obs_sequence

4. run program create_obs_sequence to create an observation sequence (usually set_def.out, although it is possible to create obs_seq.out files directly if you don’t really care about the observation values).

5. if desired, run program create_fixed_network_seq to create a longer observation sequence file.

6. run program perfect_model_obs to harvest the synthetic observations from a chosen model.

This directory contains a MATLAB function (even_sphere.m) that generates input for the program create_obs_sequence . After executing create_obs_sequence, the resulting observation sequence file will have a template for ‘RADIOSONDE_TEMPERATURE’,’RADIOSONDE_U_WIND_COMPONENT’,and ‘RADIOSONDE_V_WIND_COMPONENT’ observations at specified pressure levels and roughly evenly-distributed locations across the entire globe. Optionally, bogus observation values may also inserted; which may be useful in certain circumstances.

even_sphere.m

even_sphere.m has many optional arguments to tailor its behavior. It has exactly 1 required argument - the number of horizontal locations desired.

• It will create a text file even_create_input to be used as input to create_obs_sequence

• The choice of pressure levels is described here.

• The default number of pressure levels is 21. Argument nlevels specifies how many levels to use from the beginning of the levels list.

• The default observation error variances for each observation type are level-dependent and are consistent with DART/observations/obs_converters/obs_error/ncep_obs_err_mod.f90

• The default is to create ‘empty’ observation sequences - i.e. they have no actual observation values and are suitable to be used with perfect_model_obs

• The default date of the observations is 2017-12-25 00:00:00

• A plot of the locations will be created. The number of gridlines is configurable but defaults to 288 in longitude and 192 in latitude.

• All the defaults can be changed by specifying ‘variable-value’ pairs of options, as described below. Examples of some options are also available via the normal MATLAB help facility. (Documenting all of them in the help makes the help page too long.)

Note that the number of observations will be the number of locations * the number of vertical levels * the number of variables (i.e. 3) even_sphere.m also takes observation error variances and includes them in the observation sequences.

Optional Argument Variable-Value pairs

The optional variable-value pairs can appear in any order.

optional variable	example value	Description
‘nlevels’	5	number of pressure levels to use. May be less than the length of the ‘levels’ array, but cannot be more.
levels	[1000 500 300 200 100]	pressure levels desired. see Levels section for discussion.
T_error_var	[1.44 0.64 0.81 1.44 0.64]	level-specific Temperature error variances. see Levels section for discussion.
W_error_var	[1.96 4.41 9.00 7.29 4.41]	level-specific error variances for both U, V wind components. see Levels section for discussion.
‘YMD’	‘2017-12-25’	Date required for create_obs_sequence. If create_fixed_network_seq is run, this time is replaced.
fill_obs	false	‘true’ inserts a bogus observation value of 1.0 and a bogus QC value of 0.’false’ does not insert bogus values and essentially creates an empty obs sequence file (typically set_def.out)
‘nlon’	288	number of longitude grid lines in plot
‘nlat’	192	number of latitude grid lines in plot

Examples

1. 30 horizontal locations at 6 pressure levels:

nprofiles   = 30;
levels      = [1000  850  500  300  200  100];
T_error_var = [1.44 0.64 0.64 0.81 1.44 0.64];
W_error_var = [1.96 2.25 4.41 9.00 7.29 4.41];
even_sphere(nprofiles, 'levels', levels, ...
           'T_error_var', T_error_var, 'W_error_var', W_error_var)

2. 30 horizontal locations at 3 pressure levels. Note that the nlevels argument specifies that only the first 3 pressure levels are used even though there are 6 potential pressure levels. Similarly, only the matching error variances are used.

nprofiles   = 30;
nlevels     = 3 ;
levels      = [1000  850  500  300  200  100];
T_error_var = [1.44 0.64 0.64 0.81 1.44 0.64];
W_error_var = [1.96 2.25 4.41 9.00 7.29 4.41];
even_sphere(nprofiles, 'nlevels', nlevels, 'levels', levels, ...
           'T_error_var', T_error_var, 'W_error_var', W_error_var)

Levels

Attention

If you need realistic error variances attached to your observations, be careful to align your levels and variances.

The default levels that this program generates are the mandatory pressure levels defined in the AMS glossary. The corresponding error variances are from ncep_obs_err_mod. See Obs Error. Levels at the top can be excluded by setting nprofiles < 21 (size(levels)).

levels      = [1000  925  850  700  500  400  300   250  200  150  100   70   50   30   20   10    7    5    3    2    1];
T_error_var = [1.44 1.00 0.64 0.64 0.64 0.64 0.81  1.44 1.44 1.00 0.64 0.64 0.81 1.00 1.69 2.25 2.25 2.25 2.25 2.25 2.25];
W_error_var = [1.96 2.25 2.25 2.56 4.41 6.76 9.00 10.24 7.29 5.76 4.41 4.41 4.41 4.41 4.41 4.41 4.41 4.41 4.41 4.41 4.41];

Here’s an example of replacing the AMS levels with a set that has more levels near the surface and none above 150 hPa. Note that the error variances should change to be consistent with the levels.

levels      = [1000  950  900  850  800  750  700  650  600  550  500  400  300  200  150];
T_error_var = [1.44 1.21 0.81 0.64 0.64 0.64 0.64 0.64 0.64 0.64 0.64 0.64 0.81 1.44 1.00];
W_error_var = [1.96 2.25 2.25 2.25 2.56 2.56 2.56 3.24 3.61 4.00 4.41 6.76 9.00 7.29 5.76];

Running Matlab in Batch Mode

If you would prefer to run even_sphere.m in batch mode (i.e. from within a shell script), here is an example syntax that worked for me. The script ran in the same directory as even_sphere.m. There are many ways to construct the input, naturally - but you don’t have to explicitly edit even_sphere.m this way.

#!/bin/csh

\rm -rf matlab_input.m

cat >> matlab_input.m << EndOfInput

   nprofiles   = 30;
   levels      = [1000  850  500  300  200  100];
   T_error_var = [1.44 0.64 0.64 0.81 1.44 0.64];
   W_error_var = [1.96 2.25 4.41 9.00 7.29 4.41];
   even_sphere(nprofiles, 'levels', levels, ...
              'T_error_var', T_error_var, 'W_error_var', W_error_var)
   fname = sprintf('even_sphere_%d_profiles',nprofiles);
   orient landscape
   print(fname,'-dpdf')

EndOfInput

matlab -nosplash -nodesktop -r "try; cd $PWD; matlab_input; catch; end; exit";

Automation Scripts

Here there are also scripts (run_fixed_network_*.csh) which use the output from create_obs_sequence and the program create_fixed_network_seq to generate a series of observation sequence files.

run_fixed_network_seq.csh

Calls create_fixed_network_seq to create a separate file for each time period. By default, it makes 2 files/day, 12 hours apart, single time per file. The intervals and dates can be changed by editing the script. It assumes that create_fixed_network has any model-specific files it needs in this directory. It requires a set_def.out file (usually created by create_obs_sequence).

run_fixed_network_daily.csh

Calls create_fixed_network_seq to create a separate file for each time period. By default it makes 1 file/day, single time (noon) per file. The dates and time can be changed by editing the script. It assumes that create_fixed_network has any model-specific files it needs in this directory. It requires a set_def.out file (usually created by create_obs_sequence).

The process, end to end:

MATLAB:

Set the number of profiles, the levels, etc. and run even_sphere.m in MATLAB. It creates the necessary text file even_create_input for the next step. It will also make a plot - which you can save.

DART:

Then you have a choice about building and running the create_obs_sequence and create_fixed_network_seq programs:

1. building them in the models/template/work directory

2. using the ones which were built in models/your_model/work directory by quickbuild.sh.

Choice A uses programs which have no model specific file dependencies, but may involve more separate steps than B.

A

1. Build the programs in template/work

2. Link (or copy) these files to the directory in which you want to create obs_seq files.

./even_create_sequence
./run_fixed_network_{seq or daily}.csh
models/template/work/create_fixed_network_seq
models/template/work/create_obs_sequence
models/template/work/input.nml

3. In your obs_seq directory, run create_obs_sequence, which creates a set_def.out file.

./create_obs_sequence < even_create_input > /dev/null

4. Edit and run your choice of run_fixed_network_*.csh for the desired dates. These call create_fixed_network_seq, which creates an obs_seq.in file for each specified date.

B

This choice may involve fewer steps, if there is a model specific script which combines the steps in A). See the cam-fv example (models/cam-fv/shell_scripts/synth_obs_locs_to_seqs.csh). If there is not a script like that for your model, you can follow the steps in A), substituting your model name for the “template” in the pathnames. NOTE: you may need to link any additional input files which your model requires into the directory where you will run the programs. These typically contain grid information and are found in your_model/work. For example, cam-fv needs a caminput.nc and cam_phis.nc.

Fortran program for generating obs directly

cd into the work directory and run quickbuild.sh.

This builds the create_even_sphere executable. Edit the input.nml to set the number of obs to generate and the date in the namelist. Run the program and the output file will be generated.

DETAILS of generating points evenly-distributed on a sphere

This is the algorithm that’s being used [1]:

N     := the number of profiles you want
dlong := pi*(3-sqrt(5))  /* ~2.39996323 */
dy    := 2.0/N
phi   := 0
y     := 1 - dy/2

for k := 0 .. N-1
    r       := sqrt(1-y*y)
    node[k] := (cos(phi)*r, sin(phi)*r, y)
    y       := y - dy
    phi     := phi + dlong

For the geometric and visually minded:

1. Picture a unit sphere in cartesion space (x,y,z).

2. Choose a value -1 < y < 1, which defines an x-z plane. That plane intersects with the unit sphere to form a circle whose center is on the y axis. (The circle radius is small near y = +/-1 and is 1 at y=0.)

3. Choose an angle (“phi”) and draw a ray from the center of the circle to a point on the circle using this angle relative to the x positive direction. Where the ray intersects the circle (and sphere) is one of the evenly distributed points on the sphere which we want.

4. Its x and z coordinates can then be combined with the already defined y coordinate to define the cartesian location of the point.

5. The choice of the y and angle for each point is where the magic enters the algorithm. They are derived from the Fibonacci or Golden Spiral formula (derived elsewhere).

[1] (1,2)

A python example of the Golden Section spiral algorithm can be found in https://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere See the contribution from Fab von Bellinghousen.