ATOC 5060: Atmospheric Dynamics

This tutorial is intended to get you started on the computing assignments, and thus is very much non comprehensive, and should be considered a survival guide.

Introduction to atoc.colorado.edu

atoc is a Linux computer. You can log onto it from any machine connected to the internet using ssh.:

    ssh username@atoc.coloroado.edu

Enter your password. If you don't know what this is, or if you don't know your username, contact the system administration at trouble@atoc.colorado.edu and ask. They will likely ask you to call so that that can set (reset) your passwork over the phone.

Once logged on, you'll find yourself in you home directory and be presented with the command line prompt. Now, you can start typing. Some basic UNIX commends include cp (copy), ls (list contents of directory), pwd (tell you which directory you're in), and cd (change directory). If you are not familiar with these, see some additional hints on this help page. In all cases, when you need assistance, the fastest source of information is google. Try a google search for "basic unix help".

Before you begin the assignments, I suggest you make sure you can access your output via the web. From your home directory check that you have a subdirectory called public_html. (type ls, and if it is there, you will see it). If it does not exist, make it! type:

    mkdir public_html

Now any files you put in this directory can be seen via a web browser at:

    http://atoc.colorado.edu/~username

For the class assignments, it maybe convenient to output graphic to this directory, such that you can see them easily and copy them to word documents as needed.

Getting started with IDL

To start IDL, type, "idl". You will be presented with the IDL prompt ("IDL>"). Notice that this is the command line interface.

IDL is a rather detailed analysis package, that contains a large array of graphing and data processing functions. It also has a syntax that is somewhat similar to other languages (Fortran, C, etc). In many respects IDL is similar to Matlab.

IDL variables, assignments and operators

Variables in IDL (like other languages) can be scalars or arrays. Arrays can be though of as vectors, where each element of the vector can be referenced by an index. Arrays may also be multi dimensional (like fields that depend on, say, space and time, or x and y, etc.). Also, variable can be different types. 3 types of particular note are floating point variables ("floats"), integers, and complex numbers. Consider the difference between 1.0, 1, and (1 + i 0) (where i is the square root of 1). Mathematically, about the same! But the way these are represented by the computer differ, which we need to be careful about.

You can assign values to quantities like this

     IDL> a = 1.234

Then print its value on the screen

     IDL> print, a

Now try some math!

     IDL> print, a*a

     IDL> print, 1/a

     IDL> print, a + a

     IDL> print, a^2

you net the idea. You can also assign values:

     IDL> b = a + 2

and so on.

Rather than typing at the prompt you can put these commands in a file (called somtheing.pro, and yes, the .pro is important). The last line of the file must be "end". Make a new file with a text editor (see about text editors below), such that it contains the following:

     atoc> pico myfile.pro

     a = 1.234
     b = a+2
     print, b
     end

Now, back in idl, run this "script":

     IDL> .run myfile

(You should see IDL print the answer on the screen)

You can always check what type of variable something is by typing "help, variablename" (e.g., "help,a"). Is should say "FLOAT". This is particularly helpful if you have arrays and want to know what their dimensions are.

IDL arrays and loops

Arrays work just the same. First you must create an array, of some length. Consider the following, and add this to you myfile.pro file

     nlon = 144
     lons = fltarr(nlon)

This sets up an (empty) array of length "nlon".

Now you can fill up your array.

     lons[0] = 0.
     lons[1] = 2.5

etc etc etc... actually this gets to be very tedious... perhaps IDL can automate this... SURE CAN. With a "for loop". Consider the following.

     for i = 0,nlon-1 do begin
        lons[i] = 2.5*i
     endfor
     print, lons

Now .run myfile again, you should see a list of longitudes printed out. Notice arrays must have square brackets when referring to individual elements. Notice that the first element of the array is element number zero (i.e., lons[0]). This means the last element will be number nlon-1. So lons[nlon-1]. The line that starts the for loop says: loop over the chunk of code between the 'begin' and the 'endfor', and each time though make i larger by 1, starting from i = 0, and finish looping when i = nlon-1. This is a fundamental building block of any complicated and computationally demanding analysis.

Computing various averages and deviations

Now lets try a more useful example. Imagine you need to compute the mean of some quantity, let's imagine, say, you have a variable called "tg" which is an array with dimensions nlon and ntime. (i.e., temperature T = T(x,t)).

     tg = fltarr(nlon,ntime)
     ; Some code here to read in time temperature data from file to
     ; fill in tg with actual values

In IDL the semicolon acts as a "comment" line. Anything that follows the semicolon will be ignored when you run you script.
Now you can make a new array to save the time mean.

     tbar = fltarr(nlon)
     do i = 0, nlon-1 do begin
       tbar[i] = mean(tg[i,*])
     endfor

Here, tbar is the new time mean. You can use the IDL function "mean" to get the mean of the tg values at at particular longitude (i), over all times. Time is the second dimension of the array, and * means "use all values in this dimension". One could of course also compute the zonal mean. Perhaps,

     tzon = fltarr(ntime)
     do n = 0, ntime-1 do begin
       tzon[n] = mean(tg[*,n])
     endfor

(Check the difference in dimensions). If you want the time and zonal mean,

     tmean = mean(tg)

Now that we have some averages, lets make some anomalies. One can take the definition that we can decompose any quantity into a mean and a varying component. In space, T = [T] + T*, and in time T = Tbar + T'. And a combination of the two would give, T = [Tbar] + Tbar* + T'. Above, we have already calculated [Tbar], the time and space mean, as "tmean". Also, Tbar, is, well, "tbar". So consider the following.

     tstar = fltarr(nlon)
     for i = 0, nlon-1 do begin
       tstar[i] = tbar[i] - tmean
     endfor

(In fact IDL "knows" tbar is an array, so you could just do "tstar = tbar - tmean" as a single line an IDL will know to do the loop part. However, I recommend writing the loop explicitly, since it helps avoid errors when getting started)

Now since we know that T is the sum of the three parts, we can get the T' as a residual! That is, the transient component is the total minus the stationary part. Notice that T' = T'(x,y), so need to work in 2 dimensions.

     tprime = fltarr(nlon,ntime)
     for n = 0, ntime-1 do begin
       for i = 0, nlon-1 do begin
          tprime[i,n] = tg[i,n]  - tbar[i]
       endfor
     endfor

Where, as above, tbar = tmean + tstar. If we've done this all correctly, the time mean of tprime, should be zero! So too, the space (zonal) mean of tstar should be zero.

Fourier transforms

You may be interested in converting you "grid point" data into a representation as a series of (complex) Fourier coefficients. IDL has a built in function called FFT, which can do both the forward (grid to Fourier) and inverse (Fourier to grid) transformations. E.g.,

     tstarfour = fft(tstar)

Let's see what the type and dimensions are.

     help,tstar
     TSTAR FLOAT = Array[144]

     help,tstarfour
     TSTARFOUR COMPLEX = Array[144]

That is, FFT returns the series of complex coefficients. Notice that element zero is the "zeroth" coefficient, and thus, is just the mean value. Since tstar is a deviation, the mean is (SHOULD BE!) zero. Of the remaining coefficients, only half have useful value - recall the Fourier transform is symmetric about zero with the negative wavenumbers (frequencies, if in time) just symmetric, or, the complex conjugate. As such, to look at the amplitude the harmonics, we only need to look at half of the spectrum. The number of useful waves is nlon/2-1.

The amplitude is the square root of the real part squared plus the complex part squared. This is done for us with the IDL function "abs" (i.e., absolute value). So we can look at the amplitude as

     print, abs(tstarfour[1:nlon/2-1])

After reading about plotting (below) you will be ready to plot the amplitude as a function of wave number.

     wavenumber = findgen(nlon)
     plot,wavenumber[1:nlon/2-1], abs(tstarfour[1:nlon/2-1])

Remember that amplitude is not the same as power. Power is proportional to variance, and is the real part of the coefficients multiplied by their complex conjugate. It is useful to keep in mind you can look at the real and imaginary parts separately:

     tstarreal = real_part(tstarfour)
     tstarimag = imaginary(tstarfour)

You can get the complex conjugate as

     tstarconj = conj(tstarfour)

As such, one obtains the power as

     tstarpower = 2*real_part(tstarfour*conj(tstarfour))

Notice the factor of two is in there since we will only look at half of the spectrum.

Which you can now plot (again, look at plotting below):

     plot,wavenumber[1:nlon/2-1], tstarpower[1:nlon/2-1], /xlog

Notice, it is often convenient to make power spectrum plots with a log axis.

What are the units of power? Same as variance. So the total variance (integral over wave number) should equal the total variance in grid point form (i.e., before doing the FFT). Let's check! First the total variance in the initial (grid) data:

     print, variance(tstar)

Now let's get the variance from the power spectrum:

     print, total(tstarpower[1:nlon/2-1])

(Should be close - the difference is due to numerical precision, and should be less than about 0.1%). Notice we only do the total over the useful half of the spectrum, else we would have twice the variance because of the factor 2 above (remember the full spectrum is symmetric, so we only need to look at half).

It is a good habit to always check this to make sure you haven't lost variance! (And can give a reminder that there should be a factor of 2...). I strongly advice you convince yourself of this power/variance spectrum by making a plot of tstar as in this example, with the plotting example below before moving onto the next step as in the assignment.

For more details on the use of the FFT routine, see http://idlastro.gsfc.nasa.gov/idl_html_help/FFT.html
If your interested in inverse transforms, look at the keyword "direction". If you don't give this, it implies the forward transform. i.e., direction = -1 by default.

Plotting and seeing graphs

Ultimately you will want to look at your results graphically. IDL has a large range of plotting routines, which is one of its strengths. The graphics can be output to the screen (if you have graphics set up), or to a file (ultimately what you need, to have a graph in a paper). Let's work with sending the graphics to a postscript file. Essentially all you need to do is tell IDL to send the graphics to a postscript file, then make the graph, then close the file.

     set_plot,'ps'
     plot,lons,tstar
     device,/close

Easy! You've just made a graph of tstar, the stationary wave variation in T. By default, idl will output this to a file called "idl.ps" You can add a line to you script to automatically convert this to a gif or png file that you can view on the web. This is actually done on the UNIX command line, but you can get IDL to do it for you.

     spawn,'convert idl.ps ~/public_html/plot.png'

You can now go to http://atoc.colorado.edu/~usename/plot.png to see your figure. (Yes, you have to change username to YOUR username!).

If you want to add additional lines to the same axes, you can use oplot, to "over plot", Consider the following very messy graph that results if you have:

     set_plot,'ps'
     plot,lons,tstar, thick=5
     for n = 0, nime-1 do begin
        plot,lons,tstar[*,n]
     endfor
     device,/close
     spawn,'convert idl.ps ~/public_html/plot.png'

IDL has lots of options for the plot function that you can use to make you graph aesthetically pleasing. See the plot documentation at: http://idlastro.gsfc.nasa.gov/idl_html_help/PLOT.html. You can also delve into the world of color, but you will need to set up you PS device correctly first.

     device,/color,bits=8

You get fancy, you can also have contour plot, maps with various projections, vector field, even 3d volumes, and animations. But this is not needed for this class.

Notice here that I suggest you write the graphics to a file, and use a web browser to see it. In general you can also have the graphic output to the screen (and "X" device), but this requires you are using X windows. If you are, you can have set_plot,'X'. (Linux and Mac users, you already have this!) For out class assignments, the time taken to set up X windows environment may not be worth it.

More?

Cora Randall also has a brief tutorial for ATOC 5235. There are other good example of beginner tutorials available lots of places around the web - as always, try a google search for "beginner idl tutorial" and find one you like.

Simple text exiting

You will need to edit/create your idl scripts using a text editor. There are many of these, and they all have various quirks. For the die-hard, I highly recommend learning the editor "vi" (or its rival emacs). Notice that vi is very powerful, but there is a steep learning curve.

A more user-friendly editor is called "pico". This is the editor that is build into the good-old email program pine. You open up files from the command line as:

pico mynewfile.pro

(here the "pro" suggest that this will be an idl procedure). You can then edit as you might expect - move around wit the arrow keys, and type things, delete things. Finally, you can hit control-O to "save" or control-X to exit. It will ask if you if you want to "save changes to the buffer" if you haven't yet saved the changes. Hit Y if yes, or N if you want to not save your change.

One could now start up idl, and run your new script with ".run mynewfile". Usually I find it most convenient to have two windows open: One with IDL running, and a second with the editor running. This way, you can make changes (i.e., fix bugs) and save in one window. Then jump back to the other window to run the script in IDL.

Last updated:

ATOC 5060

Atmospheric Dynamics

Computing tutorial