Subversion Repositories lagranto.um

\documentclass[a4paper,10pt]{article}

\usepackage[utf8x]{inputenc}

\usepackage{graphicx}

\graphicspath{{/home/sprenger/lagranto/docu/tutorial/}}

\textwidth16cm

\textheight22.5cm

\oddsidemargin0.5cm

\evensidemargin0cm

\topmargin0.cm

\headsep0cm

\topskip-0.5cm

\title{Lagranto - Tutorial}

\author{Michael Sprenger and Heini Wernli}

\begin{document}

\maketitle

\section{Definition of the Problem}

In this tutorial we consider a simple way how to find airstreams which transport air from the surface into the upper troposphere and lower stratosphere, i.e. the UTLS region. More specifically, we intend

\begin{itemize}

\item[1)] to initialize starting positions over the North Atlantic and Europe (80 W to 20 E, 40 N to 80 N) for 00 UTC, 20 October 1989. The starting positions are horizontally equidistant with 80 km horizontal spacing and are set 100 hPa above ground level.

\item[2)] to calculate trajectories several hours forward in time.

\item[3)] to select trajectories which ascend to levels above 400 hPa within 48 hours and are found at starting time below 700 hPa.

\item[4)] to trace several meteorological fields along the trajectories: potential vorticity (PV), potential temperature (TH), relative (RH) and specific humidity (Q).

\item[4] to select subsamples of trajectories: a) those reaching the stratosphere; b) those travelling at most 2000 km; ... 

\item[5)] to show the densities of the trajectories on a geographical map.

\end{itemize}

\section{Meteorological Data}

The following dates and netCDF files are needed for a Lagranto calculation covering the time period from 07\,UTC 14 February 2011 to 11\,UTC 14 February 2011. 

\begin{verbatim}

> datelist stdout -create 20110214_07 20110214_11 interval 1

 20110214_07

 20110214_08

 20110214_09

 20110214_10

 20110214_11

\end{verbatim}

\noindent

There are two different files involved, the P and the S files. Lagranto expects them to be in the running directory:

\begin{verbatim}

> ls -1

 P20110214_07 P20110214_08 P20110214_09 P20110214_10 P20110214_11

 S20110214_07 S20110214_08 S20110214_09 S20110214_10 S20110214_11

\end{verbatim}

\noindent

The meteorological fields on the primary P files are at least: zonal wind (U, in m/s); meridional wind (V, in m/s); vertical wind (OMEGA, in Pa/s); 

surface pressure (PS, in hPa). Additional fields might be available on the P files, e.g. temperature (T), specific

humidity (Q),... Secondary fields can be saved in the S files, which must have the same grid structure as the P files.

In the example the following fields are saved on the S files: potential temperature (TH, in K); Ertel potential vorticity (PV, in pvu); relative humidity (in \%). 

Furthermore, the surface pressure (PS) is also saved on the S files; it must be exactly identical to the one in the P files.\\

\noindent

If surface pressure is not saved on the S files, it must be copied from the corresponding P files. This can most readily be done using the NCO tools (http://nco.sourceforge.net/):

\begin{verbatim}

> ncks -A -v PS P20110214_07 S$20110214_07

\end{verbatim}

\noindent

If the P and S files are stored at another place, they might be linked with a simple Shell script (in csh), usimg the command {\em datelist}:

\begin{verbatim}

> foreach date ( `datelist stdout -create 20110214_07 20110214_11 -interval 1` )

>   ln -s {SOURCE DIR}/P${date} {DEST DIR}/P${date}

>   ln -s {SOURCE DIR}/S${date} {DEST DIR}/S${date}

> end

\end{verbatim}

Note that the command {\em datelist} offers several options how to work with date list - creating, stepping through, comparing.

\section{Starting Positions}

In a first step the starting positions must be specified. To this aim a file {\em regionf} must be created with the definition of the region. Note that the UM uses a rotated coordinated system, whereas the starting region must be specified in true longitude and latitude in {\em regionf}. The position of the rotated north pole can easily be extracted with

\begin{verbatim}

> rot2geo -file 20110214_07

112 157.5

\end{verbatim}

\noindent

where 112.0 and 157.5 are longitude and latitude of the rotated north pole. The file {\em regionf} is in the same directory as the meteorological data (section 1):

\begin{verbatim}

> more regionf 

"1 -69.0 -67.0 -68.0 -66.0"

\end{verbatim}

\noindent

The first number specifies a region ID (here 1) and the other values are: west boundary (69 W), east boundary (67 W), southern boundary (68 S) and northern boundary (66 S). It is easily checked whether the lat/lon points are within the rotated domain. To this aim, use for instance

\begin{verbatim}

> geo2rot -69.0 -68.0 -file P20110214_07

  -0.3746045      -0.5030226

\end{verbatim}

\noindent

where now the output give the rotated longitude and latitude. These coordinates can be compared to the ones (domxmin,domxmax,domymin,domymax) on the netCDF files:

\begin{verbatim}

> ncdump -h P20110214_09

...

/ global attributes:

                :domxmin = 359.3602f ;

                :domxmax = 362.9647f ;

                :domymin = -1.5665f ;

                :domymax = 2.038f ;

                :domzmin = 1000.f ;

                :domzmax = 50.f ;

...

\end{verbatim}

\noindent

The starting positions are then created with

\begin{verbatim}

> create_startf 20110214_07 startf.2 'region.eqd(1,10) ...

                             ... @ level(100) @ hPa,agl' -changet

\end{verbatim}

\noindent

The starting positions are written to {\em startf.2}, i.e. in format 2, and cover the region 1 specified in {\em regionf}. The horizontal start points are equidistantly distributed with 10 km spacing, and they are all at 100 hPa above ground. All points refer to the starting date 07 UTC, 14 February 2011.

In total, 3750 starting positions are written to {\em startf.2}. The first few lines of the file look as follows:

\begin{verbatim}

> head -10 startf.2

Reference date 20110214_0700 / Time range       0 min

  time      lon     lat     p     level

---------------------------------------

   0.00   -68.996  -67.928   908   100.000

   0.00   -68.756  -67.928   908   100.000

   0.00   -68.516  -67.928   908   100.000

   0.00   -68.276  -67.928   908   100.000

   0.00   -68.037  -67.928   908   100.000

\end{verbatim}

\noindent

The different columns are: time, longitude, latitude, pressure (in hPa) and level (in hPa,agl). Note that the 'same' starting file could have been created without a region file. In this case, the command would have been:

\begin{verbatim}

> create_startf 19891020_00 startf.2 'box.eqd(-69,-67,-68,-66,10) ...

                             ... @ level(100) @ hPa,agl' -changet

\end{verbatim}

\noindent

However, note that the starting positions are not exactly the same as with the previous command: it is only guaranteed that the starting points are equidistantly distributed within the region.

\section{Trajectory Calculation}

In a next step, the trajectories are calculated 72 h forward in time, with starting date 00 UTC, 20 October 1989. The command is:

\begin{verbatim}

> caltra 19891020_00 19891023_00 startf.2 traj.4 -j

\end{verbatim}

\noindent

The taring positions are taken from{\em startf.2} and the output is written to {\em traj.4}. Furthermore, the jumping flag {\em -j} is set, i.e. if trajectories run into the ground they are lifted a little and allowed to move on.\\

\noindent

Note that the output file {\em traj.4} is not in ASCII format. To look at the file, use the command {\em trainfo}, for instance:

\begin{verbatim}

> trainfo traj.4 vars

time lon lat p 

> trainfo traj.4 dim

3750           13            4

> trainfo trai.4 startdate

19891020_0000

\end{verbatim}

\noindent

It is also possible to convert the trajectory file into ASCII format with the command {\em reformat}:

\begin{verbatim}

> reformat traj.4 traj.1

> more traj.1

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p

-----------------------------

   0.00   -79.61   40.45   862

   6.00   -80.57   43.23   791

  12.00   -82.23   45.89   782

  18.00   -84.94   47.07   744

\end{verbatim}

\noindent

The command {\em trafinfo} cna also be used to look at the trajectory tables, i.e. without a conversion to ASCII format. To this aim, use

\begin{verbatim}

> trainfo traj.4 list

\end{verbatim}

\noindent

Whereas the ASCII format is most convenient for visual inspection, it is the least compact format. In particular, if the output of {\em caltra} should be further processed, e.g. with {\em trace} or {\em select}, the binary formats should be used (see documentation for {\em reformat}).

\section{Pre-Selection of Trajectories}

Quite often, the position of the air parcels is sufficient to select trajectories. It is then most efficient to pre-select these trajectories and do the tracing of additional fields along the trajectories only on the pre-selected ones. In the example, airstreams should be identified which ascend from below 700 hPa at initial time to levels above 300 hPa. The ascent has to take place within 48 h. This selection can be achieved with the command:

\begin{verbatim}

> select traj.4 wcb.1 'GT:p:700:FIRST & LT:p(MIN):400:0 to 48'

\end{verbatim}

\noindent

The first criterion selects all trajectories for which the pressure (p) at the initial time (FIRST) is greater than (GT) 700 hPa. This criterison is logically AND-combined with the second criteron: consider all times between 0 and 48 h and take the minimum pressure, i.e. p(MIN), over this time interval; if the minium pressure is less than (LT) 400 hPa, the trajectory is selected. A sample trajectory looks like:

\begin{verbatim}

> trainfo wcb.1 list 

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p

-----------------------------

   0.00   -72.98   40.45   918

   6.00   -76.45   43.14   879

  12.00   -78.53   46.69   808

  18.00   -80.08   48.70   770

  24.00   -84.49   48.71   563

  30.00   -87.89   43.32   377

  36.00   -80.69   37.24   396

  42.00   -73.05   39.00   477

  48.00   -67.62   47.21   488

  54.00   -63.53   54.61   455

  60.00   -53.79   58.53   447

  66.00   -38.79   59.08   452

  72.00   -27.72   55.51   493

\end{verbatim}

\noindent

In total, 99 trajectories are selected. Further Lagrangian selection criteria might be reasonable. For instance, it could be of interest whether the air parcels are far away from their initial position after 48 h:

\begin{verbatim}

> select traj.4 wcb.1 'GT:p:700:FIRST & LT:p(MIN):400:0 to 48 & GT:DIST0:5000:48'

\end{verbatim}

\noindent

The last criterion test whether the spherical distance from the initial position exceeds 5000 km at 48 h (only met by 4 trajectories). Note that the field {\em DIST0} is not available on the trajectory file {\em traj.4}, but is implicitely calculated. \\

\noindent

Similarly, it can be tested whether a trajectory passes through a target region (e.g. 20E-30E,50N-60N). Such a region might be defined in the region file {\em regionf}:

\begin{verbatim}

> more regionf 

# Starting positions

"1 -80  20 40 80"

# Target region

"2 20 30 50 60"

\end{verbatim}

\noindent

The call to {\em select} now looks as follows:

\begin{verbatim}

> select traj.4 wcb.1 'GT:p:700:FIRST & LT:p(MIN):400:0 to 48 & ...

                       ... TRUE:INREGION:2:42 to 54(ANY)'

\end{verbatim}

\noindent

The last criterion is interpreted in the following way: consider the times from 42 h to 48 h (42 to 54) and check whether a trajectory is at any of these times (ANY) in the traget region 2, as specified in {\em regionf}. A sample trajectory is given below:

\begin{verbatim}

> more wcb.1

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p

-----------------------------

   0.00   -44.56   40.45   914

   6.00   -41.22   38.95   885

  12.00   -37.50   37.59   849

  18.00   -33.54   36.76   823

  24.00   -29.45   36.45   770

  30.00   -24.56   37.44   622

  36.00   -18.89   41.17   429

  42.00   -10.48   49.48   349

  48.00     8.75   57.02   355

  54.00    28.74   56.97   354

  60.00    37.14   53.71   320

  66.00    38.99   49.75   334

  72.00    37.57   45.94   349

\end{verbatim}

\section{Tracing Meteorological Fields}

Meteorological fields can be traced along the trajectories with the command {\em trace}. Most often, a list of fields to trace will be listed in a file {\em tracevars}:

\begin{verbatim}

> more tracevars 

PS    1.  0 P

Q  1000.  0 P

TH    1.  0 S

RH    1.  1 *

\end{verbatim}

\noindent

The following fields are to be traced: surface pressure (PS), specific humidity (Q), potential temperature (TH) and relative humidity (RH). PS and Q are available on the P files and need not to be calculated; Q will be scaled by a factor 1000 to convert from Kg/kg to g/kg. TH is found on the S file and need not to be calculated. Finally, RH is found neither on the P nor on the S file and must becomputed - the flag 1 in the third column.\\

\noindent

With the {\em tracevars} file ready, the tracing is started with:

\begin{verbatim}

> trace wcb.1 wcb.1

\end{verbatim}

\noindent

Note that the input and output file are allowed to have the same name. The following table shows the first few lines of the new trajectory file:

\begin{verbatim}

> more wcb.1 

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p        PS         Q        TH        RH

---------------------------------------------------------------------

   0.00   -44.56   40.45   914  1014.093     9.664   294.921    86.667

   6.00   -41.22   38.95   885  1012.965     8.380   296.549    77.992

  12.00   -37.50   37.59   849  1014.056     8.565   299.122    81.321

  18.00   -33.54   36.76   823  1012.074     8.720   300.798    85.407

  24.00   -29.45   36.45   770  1012.254     7.666   303.487    85.262

\end{verbatim}

\noindent

If it is later found that additional fields should be traced, this can be done with a new {\em tracevars} file or for a single field with (the second number being the scaling factor):

\begin{verbatim}

> trace wcb.1 wcb.1 -f PV 1.

\end{verbatim}

\noindent

Note that several fields are allowed for online computation, i.e. with the computation flag set in {\em tracevars}. This is convenient for interactive mode and for few trajectories. However, if tracing is needed for many trajectories, a pre-calculation and saving on the S files is much more efficient! A list of fields for online computation is found in the reference guide for {\em trace}. \\

\noindent

It is also possible to trace the surrounding of a trajectory position, i.e. to get for instance not the temperature at the air parcel's position, but at 50 hPa above or below it. This is done with:

\begin{verbatim}

> trace wcb.1 wcb.1 -f T:-50HPA 1.

> trace wcb.1 wcb.1 -f T:+50HPA 1.

\end{verbatim}

\noindent

Finally, if it is decided that a field is no longer needed in the trajectory file, or if it has to be corrected, it is possible to remove columns from the trajectory file. This can be achieved with {\em extract} - for instance if only PS, TH and RH should be kept:

\begin{verbatim}

> extract wcb.1 wcb.1 -var PS TH RH

\end{verbatim}

\noindent

Note that {\em extract} can also be used to extract different times or starting positions (see the reference documentation). 

\section{Final Selection of Trajectories}

In this section the selection of trajectories should be refined, i.e. it is not only based on positional information but also on further meteorological parameters. We look at several questions:\\

\begin{itemize}

\item[a)] Is there a trajectory which reaches saturation ($RH>99$\%)? The trajectories should be saved in a new trajectory file.

\begin{verbatim}

> select wcb.1 sat.1 'GT:RH:99:0 to 72(ANY)'

> more sat.1

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p        PS         Q        TH        RH        PV   

---------------------------------------------------------------------------------

   0.00   -72.98   40.45   918  1018.161     9.503   292.020   100.722     0.920    

   6.00   -76.45   43.14   879   986.319     8.723   294.933    92.837     1.101    

  12.00   -78.53   46.69   808   972.550     7.621   297.875    97.737     0.794    

  18.00   -80.08   48.70   770   973.957     6.147   297.914    97.912     1.078    

  24.00   -84.49   48.71   563   962.279     2.327   307.548    87.923     1.034    

  30.00   -87.89   43.32   377   977.415     0.319   314.210    65.759     0.108    

  36.00   -80.69   37.24   396   939.606     0.303   312.705    52.845     0.323    

  42.00   -73.05   39.00   477  1013.693     0.298   314.614    16.248     0.309    

  48.00   -67.62   47.21   488   970.025     0.442   312.975    23.890     0.463    

  54.00   -63.53   54.61   455   950.011     0.386   313.047    30.182     0.479    

  60.00   -53.79   58.53   447  1007.039     0.392   311.951    36.578     0.487    

  66.00   -38.79   59.08   452  1006.532     0.319   311.316    29.286     0.443    

  72.00   -27.72   55.51   493  1009.871     0.279   311.428    15.950     0.513    

\end{verbatim}

\item[b)] Get a list of all trajectories which pass through a circle around 20\,W/40\,N and radius 500\,km.

\begin{verbatim}

> select wcb.1 indlist 'TRUE:INCIRCLE:-20,40,500:ALL(ANY)' -index

> more indlist

            4

            5

            6

           11

           12

           13

           14

           19

           20

           21

           22

           47

\end{verbatim}

Hence, the trajectories 4,5,... pass through the circle. The trajectories themselves can be extracted in a second step with

\begin{verbatim}

> extract wcb.1 pass.1 -index indlist

\end{verbatim}

where now the selected trajectories are written to the trajectory file {\em pass.1}.

\item[c)] Select all trajectories which are in the stratosphere after 48 h. The dynamical tropopause is defined as the 2-PVU isosurface?

\begin{verbatim}

> select wcb.1 out 'GT:PV:2:48 to 72(ALL) & LT:P:500:48 to 72(ALL)' 

\end{verbatim}

The second criterion guarantees that the air parcel is at a height above 500 hPa; indeed, low-level high-PV regions might mimick a stratosphere, although they are of diabatic origin.  

\item[d)] Select all trajectories which are within 2000 km distance of their starting position after 72 h.

\begin{verbatim}

> select wcb.1 sel.1 'LT:DIST0:2000:LAST'

\end{verbatim}

Note that the fields {\em DIST0} needs not to be available on the trajectory file - it is calculated during the selection. {\em DIST0} refers to the spherical distance (in km) from the strting position. If the path length (in km) is needed, {\em DIST} can be used instead.

\item[e)] How many trajectories ascend more than 550 hPa between 12 h and 54 h? We are only interested in the number of selected trajectories.

\begin{verbatim}

> select wcb.1 count 'GT:P(DIFF):550:12,54' -count

> more count

3

\end{verbatim}

\item[f)] We would like to select all trajectories which reach potential vorticity (PV) greather than 2 PVU at levels above 500 hPa. In a first attempt, this might be accomplished with the criterion

\begin{verbatim}

> select wcb.1 wcb.1 'GT:PV:2:ALL(ANY) & LT:p:500:ALL(ANY)' 

\end{verbatim} 

But note that this is not exactly what we want - the first criterion might be fulfilled at a time 48 h, for instance, whereas the second criterion is fulfilled at another time, say 72 h. Hence they are not both fulfilled at the same time! A way around this problem is possible if a TRIGGER column is used to mark the two events. The original trajectory file looks as follows:

\begin{verbatim}

> more wcb.1

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p        PS          RH        PV

------------------------------------------------------------

   0.00   -19.56   46.94   905  1005.242     83.514     0.291

   6.00   -14.72   48.17   892   999.182     88.325     0.242

  12.00   -10.58   50.53   862   993.145     97.718     0.293

  18.00    -7.22   53.02   792   972.076     99.216     0.738

  24.00    -3.71   55.89   724   956.135     93.218     1.076

  30.00    -0.19   58.87   629   971.334     70.088     1.076

  36.00     1.46   61.62   452   966.406     66.056     0.558

  42.00     0.01   62.49   328   977.209     65.319     1.754

  48.00    -1.54   63.41   313   983.930     56.822     2.727

  54.00    -3.59   64.77   322   984.627     58.328     1.874

  60.00    -9.91   66.07   323   988.185     57.894     2.052

  66.00   -20.91   66.02   316   976.560     57.989     2.565

  72.00   -28.89   66.19   319  1007.175     54.477     2.693

\end{verbatim}

Then we mark the two events with a TRIGGER:

\begin{verbatim}

> select wcb.1 wcb.1 'GT:PV:2:1(TRIGGER) & LT:p:500:2(TRIGGER)' -trigger

\end{verbatim}

The first criterion (PV) gets the trigger 1 (in binary system 01), the second one (pressure) get the trigger 2 (in binary system 10). If both criteria are fulfilled, the trigger column becomes 3, which corresponds in the binary system to 11 - i.e. each flag corresponds to a bit in the trigger value. With the option '-trigger' the trigger column is written to the output trajectory file:

\begin{verbatim}

> more wcb.1

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p        PS         RH        PV     TRIGGER

-------------------------------------------------------------------------

   0.00   -19.56   46.94   905  1005.242    83.514     0.291      0.000

   6.00   -14.72   48.17   892   999.182    88.325     0.242      0.000

  12.00   -10.58   50.53   862   993.145    97.718     0.293      0.000

  18.00    -7.22   53.02   792   972.076    99.216     0.738      0.000

  24.00    -3.71   55.89   724   956.135    93.218     1.076      0.000

  30.00    -0.19   58.87   629   971.334    70.088     1.076      0.000

  36.00     1.46   61.62   452   966.406    66.056     0.558      2.000

  42.00     0.01   62.49   328   977.200    65.319     1.754      2.000

  48.00    -1.54   63.41   313   983.930    56.822     2.727      3.000

  54.00    -3.59   64.77   322   984.627    58.328     1.874      2.000

  60.00    -9.91   66.07   323   988.185    57.894     2.052      3.000

  66.00   -20.91   66.02   316   976.560    57.989     2.565      3.000

  72.00   -28.89   66.19   319  1007.175    54.477     2.693      3.000

\end{verbatim}

Now the selection can be achieved by refering to the TRIGGER column:

\begin{verbatim}

> select wcb.1 wcb.1 'ALL:TRIGGER:1,2:ALL(ANY)'

\end{verbatim}

This selection means that the trigger values 1 and 2 must be set - the operator ALL (first term) guaranteeing that all selected triggers are set. The time specification ALL(ANY) is as before, i.e. the check is performed for all times and he criterion must be fulfilled at any of these times.\\

\item[g)] Select all trajectories which pass at time 60 h over Switzerland! The coordinates of the Swiss boundary are listed in a file {\em borders.dat}:

\begin{verbatim}

> more borders.dat

8.55 47.45

7.942863 46.002075

7.949024 46.001195

7.956945 46.000022

7.984226 46.000022

7.989800 46.001489

8.000068 46.007356

8.011508 46.018503

...

\end{verbatim}

The first line is a point (longitude, latitude) within Switzerland (Zurich), the other lines define the boundary of Switzerland (as 1373 points). With this polygon file, the selection command becomes

\begin{verbatim}

> select wcb.1 out.1 'TRUE:INPOLYGON:borders.dat:60'

\end{verbatim} 

Note that in a criterion only one polygon can be specified.

\item[g)] New criteria can easily be implemented into the Fortran code; to this aim the file {\em special.f} in directory {\em select} must be edited. The following example shows the implementation of an identification for warm conveyor belts (WCB). The calling sequence for the criterion is {\em SPECIAL:WCB:ascent,first,last}, the air stream must ascend at least {\em ascent} hPa between time {\em first} and time {\em last}. The corresponding Fortran looks as follows:

\begin{small}

\begin{verbatim}

> more select/special.f

      SUBROUTINE special (flag,cmd,tra,ntim,ncol,

     >                    vars,times,param,nparam)

c     ***************************************************************************

c     *                                                                         *

c     * OUTPUT:  flag           -> 1 if trajectory is selected, 0 if not        *

c     *                                                                         *

c     * INPUT:   cmd            <- command string (e.g. WCB)                    *

c     *          tra(ntim,ncol) <- single trajectory: indices time,column       *

c     *          ntim           <- number of times                              *

c     *          ncol           <- number of columns (including time,lon,lat,p) *

c     *          vars(ncol)     <- names of columns                             *

c     *          times(ntim)    <- List of times

c     *          param(nparam)  <- parameter values                             *

c     *          nparam         <- number of parameters                         *

c     *                                                                         *

c     ***************************************************************************

      implicit none

c     ---------------------------------------------------------------------------

c     Declaration of subroutine parameters

c     ---------------------------------------------------------------------------

      integer       flag           ! Boolean flag whether trajectory is selected

      character*80  cmd            ! Command string

      integer       ntim,ncol      ! Dimension of single trajectory

      real          tra(ntim,ncol) ! Single trajectory

      character*80  vars(ncol)     ! Name of columns

      real          times(ntim)    ! List of times

      integer       nparam         ! # parameters

      real          param(nparam)  ! List of parameters

c     ---------------------------------------------------------------------------

c     Declaration of local variables

c     ---------------------------------------------------------------------------

      integer       i

      integer       ip,i0,i1

c     --------------------------------------------------------------------------  %)

c     SPECIAL:WCB:ascent,first,last                                               %)

c         : Detect Warm Conveyor Belts (WCB); the air stream must ascend at least %)

c         : <ascent=param(1)> hPa between the two times <first=param(2)> and      %)

c         : <last=param(3)>. Note, the lowest pressure is allowed to occur at any %)

c         : time between <first> and <last>.                                      %)

c     --------------------------------------------------------------------------- %)

      if ( cmd.eq.'WCB' ) then

c        Reset the flag for selection

         flag = 0

c        Pressure is in the 4th column

         ip = 4

c        Get indices for times <first> and <last>

         i0 = 0

         i1 = 0

         do i=1,ntim

            if ( param(2).eq.times(i) ) i0 = i

            if ( param(3).eq.times(i) ) i1 = i

         enddo

         if ( (i0.eq.0).or.(i1.eq.0) ) then

            print*,' ERROR: invalid times in SPECIAL:WCB... Stop'

            stop

         endif

c        Check for ascent 

         do i=i0+1,i1

            if ( ( tra(1,ip)-tra(i,ip) ) .gt. param(1) ) flag = 1

         enddo

      endif

c     ---------------------------------------------------------------------------

      end

\end{verbatim}

\end{small}

\end{itemize}

\section{Trajectory Densities}

Single trajectories can be visualised e.g. with Matlab or wit NCL (see template scripts in the Lagranto folder). If many trajectories should be visualised instead, it is much more convenient to show trajectory densities. The easiiest way to get trajectory densities is:

\begin{verbatim}

> density wcb.1 densisty

> ncview density 

\end{verbatim}

This will project the trajectories in the trajectory file {\em wcb.1} onto a global longitude/latitude grid with 1 degree horizontal resolution. A filter radius of 100 km will be used

\includegraphics[width=0.9\textwidth]{screen1.ps}

\noindent

The CF-netCDF file contains several fields: a) the number of trajectory points associated to each grid point (COUNT); b) the residence time of the trajectories (in hours) associated to each grid point - the residence time being the time a trajectory stays at a certain grid cell (RESIDENCE); c) the area (im $km^2$) associated with each grid cell. The area allows to change the unit of the gridded trajectory from counts per grid point to counts per $km^2$.\\

\noindent

Often the trajectories do not spread over the whole globe; then the subdomain can be specified with

\begin{verbatim}

> density traj.1 density -latlon 300 150 -100 10 0.5 0.5 -create

\end{verbatim}

where the new grid has 300x150 grid points in zonal and meridional direction with south-eastern corner at 100\,W/10\,S and resolution of 0.5 degree in zonal and 0.5 degree in meridional direction. The flag {\em create} forces the netCDF file to be created anew, even if it already exists.

\includegraphics[width=0.9\textwidth]{screen2.ps}

\noindent

It is also possible to re-parameterise the trajectories before they are gridded. For instance, the following command interpolates the positions to a 1-h time interval and then performs the grissing. This results in a smoother density plot:

\begin{verbatim}

> density traj.1 density -latlon 300 150 -100 10 0.5 0.5 -create -interp 1 h

\end{verbatim}

\includegraphics[width=0.9\textwidth]{screen3.ps}

\noindent

In addition to a gridding of the complete trajectories, the single trajectory times can be gridded. 

\begin{verbatim}

> density traj.1 density -create -time  0.00 -create

> density traj.1 density -create -time  6.00 

> density traj.1 density -create -time 12.00 

> density traj.1 density -create -time 18.00 

> density traj.1 density -create -time 24.00 

\end{verbatim}

\noindent

In the previous figures, the density of the trajectories was determined - i.e. the number of trajectories associated with a grid point or the residence associated with a grid cell was determined. In addition to this most basic information, it is also possible to perform a gridding of any trajectory field. For instance, the trajectory file contains

\begin{verbatim}

Reference date 19891020_0000 / Time range    4320 min

  time      lon     lat     p

-----------------------------

   0.00   -79.61   40.45   862

   6.00   -80.57   43.23   791

  12.00   -82.23   45.89   782

  18.00   -84.94   47.07   744

\end{verbatim}

\noindent

and we would like to know at what height the trajectories typically (in the mean) are at a specific grid point. Then we would grid the pressure ``p'' instead of the position:

\begin{verbatim}

> density traj.1 density -create -latlon 300 150 -100 10 0.5 0.5 -field p -time 0.00

> density traj.1 density  -field p -time 24.00

> density traj.1 density  -field p -time 48.00

\end{verbatim}

\noindent

The following figures show the gridded pressure at time 0.00 and 24.00 h; note that the trajectories were initialised 100 hPa above ground.

\noindent

\begin{center}

\includegraphics[width=0.8\textwidth]{screen4.ps}

\\

\includegraphics[width=0.8\textwidth]{screen5.ps}

\end{center}

\section{Interface Script}

\subsection{Start from local directory}

In addition to the programs {\em create\_startf}, {\em caltra}, {\em trace}, {\em select}, Lagranto offers a ``master'' script which combines the call to the individual programs into one single call. For instance,

\begin{verbatim}

> lagranto local 19891020_00 19891024_18 startf nil -changet

\end{verbatim}

will start a Lagranto run starting from 00,UTC 20 Octiober 2010 to  18,UTC 24 Octiober 2010, based upon the starting positions in the file {\em startf}. No selection of trajectories is applied, as specified with the flag {\em nil}, and the times on the netCDF P and S files are set relative to the starting date prior  to the Lagranto run. Finally, {\em local} means that all input files are expected in the directory where Lagranto was called.\\

\noindent

The output for the above Lagranto call is saved in a newly created directory, which is located in the calling directory:

\begin{verbatim}

> ls -l ntr_19891020_00_f114_local_startf_nil/

-rw-r--r-- 1 michaesp wheel 5328945 2011-03-21 14:03 lsl_19891020_00

-rw-r--r-- 1 michaesp wheel   68195 2011-03-21 14:03 runscript.logfile

-rwxr--r-- 1 michaesp wheel    1025 2011-03-21 14:02 runscript.sh*

\end{verbatim}

The three different files are:

\begin{itemize}

\item[a)]  {\bf lsl\_19891020\_00}: the output trajectory file -the file name starts with {\em lsl} and contains the starting date of the Lagranto run:

\begin{verbatim}

> more lsl_19891020_00

Reference date 19891020_0000 / Time range    6840 min

  time      lon     lat     p        PS         Q        TH        RH

---------------------------------------------------------------------

   0.00   -79.61   40.45   862   961.659     6.434   290.838    97.722

   6.00   -80.57   43.23   791   980.984     5.334   293.824    98.773

   ...

\end{verbatim}

Note that additional fields have been traced along the trajectories, as specified in the tracing file {\em tracevars}:

\begin{verbatim}

> more tracevars

PS    1.  0 P

Q  1000.  0 P

TH    1.  0 S

RH    1.  1 *

\end{verbatim}

\item[b)] {\bf runscript.logfile}: a log file with all status and error information of the Lagranto run. If the flag {\em -log} is set in a Lagranto call, the log will be written to screen.

\item[c)] {\bf runscript.sh}: the calling script for the programs {\em create\_startf}, {\em caltra}, {\em trace} and {\em select}. The basic idea of {\em lagranto} is to create the output directory, to prepare all netCDF and other files in this output directory and to create a Shell script with name {\em runscript.sh}. If all these preparations were successfull, Lagranto will change into the output directory and launch {\em runscript.sh}. The {\em runscript.sh} for the previous Lagranto call looks as follows:

\begin{small}

\begin{verbatim}

#!/bin/csh

#

#----- Calling command

#

# lagranto local 19891020_00 19891024_18 startf nil -changet -log 

#

#----- Output file    

#

# lsl_19891020_00                

#

#------ Abort if no startf is available

#

if ( ! -f startf ) then

  echo " ERROR: no start file available .... Stop"

  exit 1

endif

#

#------ Remove existing trajectory files

#

if ( -f lsl_19891020_00.4 ) then

  \rm -f lsl_19891020_00.4

endif

if ( -f lsl_19891020_00 ) then

  \rm -f lsl_19891020_00

endif

#

#------ Run <caltra>

#

/home/sprenger/lagranto//bin/caltra.sh 19891020_00 19891024_18 startf lsl_19891020_00.4

#

#------ Abort if caltra was not successful

#

if ( ! -f lsl_19891020_00.4 ) then

  echo " ERROR: caltra failed .... Stop"

  exit 1

endif

#

#------ Run <trace>

#

/home/sprenger/lagranto//bin/trace.sh lsl_19891020_00.4 lsl_19891020_00 -v tracevars

#

#------ Abort if trace was not successful

#

if ( ! -f lsl_19891020_00 ) then

  echo " ERROR: trace failed .... Stop"

  exit 1

endif

\end{verbatim}

\end{small}

\end{itemize}

\noindent

Note that you are free to change to the output directory and manually launch {\em runscript.sh}, possibly after having modified it to your needs. This way of working is uspported by the optional flag {\em -prep} which will only prepare all files and then changes to the output directory:

\begin{verbatim}

> lagranto local 19891020_00 19891024_18 startf nil -changet -prep

\end{verbatim}

At the end of this call you will be asked to change to the output directory, which -after having agreed- will open a new {\em xterm} window. Note that you can always easily change to a output directory by calling

\begin{verbatim}

> lagranto -open local 

\end{verbatim}

If several trajectory runs are available in the local directory, you are asked to select one. Often, you would like to see the outcome of a run without changing to the output directory. This is most easily accomplished with the following call:

\begin{verbatim}

> lagranto -show local 

\end{verbatim}

\subsection{Start from case directory}

In this calling sequence, for instance

\begin{verbatim}

> lagranto tutorial 19891020_00 19891024_18 startf nil -changet

\end{verbatim}

the input files are not expected in the local directory, but are specified by means of a case identifier. For instance, a case has the identifier {\em tutorial}. Then Lagranto will expect the input netCDF P and S files to be located in

\begin{verbatim}

> ls -l ${HOME}/cdf/tutorial

/home/sprenger/cdf/tutorial/P19891020_00

/home/sprenger/cdf/tutorial/P19891020_06

/home/sprenger/cdf/tutorial/P19891020_12

/home/sprenger/cdf/tutorial/P19891020_18

/home/sprenger/cdf/tutorial/P19891021_00

/home/sprenger/cdf/tutorial/P19891021_06

/home/sprenger/cdf/tutorial/P19891021_12

/hom

\end{verbatim}

and all the other input files (starting positions, tracing file, region file, polygon specification) are expected in 

\begin{verbatim}

> ls -l ${HOME}/tra/tutorial

startf

tracevars

\end{verbatim}

The output of the trajectory calculation will be written to the following output directory, where now the case identifier {\em tutorial} is part of the directory name:

\begin{verbatim}

> cd /home/michaesp/tra/tutorial/ntr_19891020_00_f114_tutorial_startf_nil

> ls -1

-rw-r--r-- 1 michaesp wheel 5328945 2011-03-21 14:03 lsl_19891020_00

-rw-r--r-- 1 michaesp wheel   68195 2011-03-21 14:03 runscript.logfile

-rwxr--r-- 1 michaesp wheel    1025 2011-03-21 14:02 runscript.sh*

\end{verbatim}

All other aspects are identical to the ones described in the previous section. 

\section{Installation}

In this section you will find some hints how to install Lagranto on a Linux platform. Everthing is handled with the installation script {\em install.sh} which comes with the Lagranto distribution:

\begin{verbatim}

> install.csh 

install.sh [lib|core|goodies|links|all]

\end{verbatim}

The installation should proceed in several distinct steps:

\begin{itemize}

\item[a)] Find the place of the netCDF (http://www.unidata.ucar.edu/software/netcdf/) installation on your system - note that the netCDF comes as a pre-compiled package for many Linux distributions and most often can be installed with the Linux software management. Define an environmental variable {\em NETCDF} which directs to your installation, e.g.

\begin{verbatim}

> setenv NETCDF /usr/local/netcdf/

\end{verbatim}

\item[b)] Set the environmental variable {\em LAGRANTO} to the place where you have stored the Lagranto source code and include Lagranto in your search path. In {\em csh} this might look as follows:

\begin{verbatim}

> setenv LAGRANTO   /home/michaesp/lagranto/

> set isLAGRANTO=`echo $PATH | grep $LAGRANTO | wc -l`

> if ( $isLAGRANTO == 0 ) then

>   setenv PATH $LAGRANTO/bin:${PATH}

> endif

\end{verbatim}

You might include these statements also in your {\em .cshrc} file. If successful, you will then be able to open the lagranto help, e.g. with

\begin{verbatim}

> lagrantohelp

\end{verbatim} 

\item[c)] Install the different components of Lagranto and create links - proceed step by step to ensure that each one was successfully completed:

\begin{verbatim}

> install lib

> install core

> install goodies

> install links

\end{verbatim} 

\end{itemize}

Lagranto should now be ready to run! As a next step you might want to consider the tutorial, which can be invoked with the command:

\begin{verbatim}

> lagrantohelp tutorial

\end{verbatim}

If you are familiar with the most basic aspects of Lagranto, please refer to the reference guide which enlists all options of Lagranto:

\begin{verbatim}

> lagrantohelp refernce

\end{verbatim}

The contents of the reference guide can also be called from within the Linux shell, e.g. the documentation of {\em caltra} can be seen in man page format with:

\begin{verbatim}

> lagrantohelp caltra 

\end{verbatim}

\end{document}