Meteorology 301 - Lab 5

Introduction

GOAL: Compute and plot p(z) for an atmosphere where T = T(z), using more advanced Matlab function techniques.

Instructions

1. As with Lab 3, follow the instructions for accessing meteorology computers to log into your Meteorology account. If you are sharing a computer, the second student should open a second window, as per the instructions, and then login as her/himself, entering in the second window,

ssh LOGNAME@0

where the student's login name should be substituted for LOGNAME. Then take turns completing the rest of the lab.

2. First, we want to create a new function file, height.m, that will give us a vector of heights z. The function will be given the number of vector elements, N, and the spacing between height levels, Dz. Our formula will be

      z(i) = i * Dz ,					(5.2.1)

• Open an emacs window:
  emacs height.m
• Then enter the following lines (being careful with upper and lower case lettering!):

  
  function [z]=height(N,Dz)
  
  % Assume z starts at z = 0
  % Then loop N times to get an array of
  % z values extending upward from the
  % surface with spacing Dz
  
  for i = 1 : N
  
  % Note semi-colon (;) at end of line.  This tells
  % Matlab NOT to echo its computation back to the
  % screen.
  
     z(i)=(i-1)*Dz ;
  
  % End loop
  
  end
  

Note a new operation here: the code will repeat the line z(i)=(i-1)*Dz for N times. This is known as a "loop", which extends from the line "for i = 1 : N" to the line "end". The iteration variable i starts with an initial value i = 1. It is incremented by i each pass through the loop. For the final pass, i = N.
• Finally exit emacs
  C-x C-c

3. Next, we want to create a new function file, temperature.m, that will give us a vector of temperatures vs. height. We will assume constant lapse rate, gamma. Our formula will be

      T(i) = Ts + z(i)*gamma ,					(5.3.1)

• Once again, open an emacs window:
  emacs temperature.m
• Then enter the following lines (again being careful with upper and lower case lettering!):

  
  function [T]=temperature(Ts,gamma,z)
  
  % Compute T(z) given the surface temperature Ts
  % and the lapse rate gamma.
  
  T = Ts - gamma*z ;
  

• Then exit emacs
  C-x C-c

4. Now let's go into Matlab to see how these functions work for us.
Start matlab:
matlab
If you are logged into Project Vincent, but not on the meteorology computers, then you have to take some extra steps to start Matlab. Enter
add matlab
add math
matlab
Do these "adds" ONLY if you are logged into Project Vincent using non-meteorology computers!

A window with a fancy 3-d graph may appear and then disappear. Then you may see a symbol like an upside-down "L" appear at your mouse pointer on the screen. Click your mouse button, and a window with 3 panels will open.

You may need to move the window to center it. This can be done by holding down the mouse button when the mouse pointer is on the upper bar of the window and dragging the window to where you want it.

5. First make sure that Matlab is working in your home directory (where presumably all your Matlab files are located):
>>cd ~
Assign an array of z(i), using a relatively coarse spacing of 1000 m (1 km):
>>N=11
>>Dz=1000
>>z=height(N,Dz)

Answer these questions:

1. What values of z are returned to your screen?
2. Why did we assign N=11 in order to get our highest z at 10 km?

6. Now assign z(i) using a much finer spacing:
>>N=51
>>Dz=200
>>z=height(N,Dz)
We will use this finer resolution for integrating numerically the hydrostatic equation.

7. Now that z is assigned values, test our function temperature.m

• Assign some surface temperature values with units [K]:
  >>Ts1=250
  >>Ts2=288
  
  and some lapse rate values with units [deg/m]:
  >>gamma1=6.5e-3
  >>gamma2=9.8e-3
  Note that these gamma's are given values in exponential notation . Thus,

                        -3
       3.0e-3 = 3.0 x 10
  

  Note also that Ts2 and gamma1 are approximately the global-average values of surface temperature and lapse rate, respectively.
• Then create some temperature vectors:
  >>TEM1=temperature(Ts1,gamma1,z)
  >>TEM2=temperature(Ts1,gamma2,z)
  >>TEM3=temperature(Ts2,gamma1,z)
  
• Plot your T(z) profiles:
  >>plot(TEM1,z,'y-',TEM2,z,'c--',TEM3,z,'y-.')
• Answer these questions:
  1. What is the effect of changing Ts on the temperature profile?
  2. What is the effect of changing lapse rate on the temperature profile?

8. Now that we have z and T(z), we are ready to compute p(z).
We will compute p(z) by integrating the hydrostatic equation up from the surface. If we have

	dp
	-- = -g*rho
	dz
then
	dp = -g*rho*dz = -g*p*dz/(Rd*Tv)
or
	d(lnp) = -g*dz/(Rd*Tv)

This formula can be used to tell us how log(p) changes over an interval dz. We will assume here that for the virtual temperature, Tv ~ T, and that g = g0 = 9.8 m/s**2.

I have written a function file to do this integration. You should copy it to your home directory:

• Exit Matlab
  >>quit
• Copy the file pressure_2.m from my MT031 directory on the Meteorology computer system to your home directory
  cp /mthome/gutowski/MT301/pressure_2.m ~

  That last character in the line is the tilde, ~. Note that there is a blank space before it. The "~" refers to your home directory. No matter where you are in the file space, you can always copy files back to your home directory by telling the computer to copy to "~".

9. Take a look at pressure_2.m
emacs pressure_2.m

Answer these questions

1. Why do we convert p(1) to lnp(1)?
2. Why do we use "{T(i-1)+T(i)}/2" as the temperature in the formula for the change in log(p) between z(i-1) and z(i)?
3. Why do we have the line "p(i)=exp(lnp(i))" in the code?

10. Now restart Matlab.
Now enter (note the finer grid spacing)
>>N=51
>>Dz=200
>>Ts1=250
>>Ts2=288
>>gamma1=6.5e-3
>>gamma2=9.8e-3
>>z=height(N,Dz)
>>TEM1=temperature(Ts1,gamma1,z)
>>TEM2=temperature(Ts1,gamma2,z)
>>TEM3=temperature(Ts2,gamma1,z)

Then compute p(z) profiles using the function pressure_2.m:
>>PR1=pressure_2(N,Dz,TEM1)
>>PR2=pressure_2(N,Dz,TEM2)
>>PR3=pressure_2(N,Dz,TEM3)

Plot the p(z) profiles you just created:
>>plot(PR1,z,'y-',PR2,z,'c--',PR3,z,'y-.')

Answer these questions

1. What is the effect on p(z) of changing surface temperature (Ts_)?
2. What is the effect on p(z) of changing lapse rate?
3. Using gamma1, what approximate range of Ts values give PR at 5 km ranging between 400 and 600 mb?

(Hint: you may find this easiest to do using plots of p(z). I only need approximate Ts values.)

EXTRA CREDIT: As in Lab 3, use "help" to learn how to put on your plot axis labels, a title, and text labels for each curve. Then add any or all of these to the plot.

Finally, whether or not you do the extra credit, save the current plot by "printing" it to a file called pressure_2.ps:
>>print -dps pressure_2.ps
As in Lab 3, this will create a file called pressure.ps in the directory you are working in. (The file is in so-called postscript format.) Following instructions given in Lab 3, step 14., print a copy of your plot..

11. Finally, compare your p(z) when T=T(z) with the p(z) computed for an isothermal temperature field, using the pressure.m file you created in Lab 3:
>>PR4=pressure(z,Ts2)
>>plot(PR3,z,'y-',PR4,z,'c--')
Here, we are assuming that geopotential height Z = distance above the surface z, a good approximation for the troposphere.

Answer these questions

1. How do the PR3 and PR4 profiles differ from each other?
2. What causes this difference? (Tell me how the atmospheric fields used to compute p(z) in each case differ and how this difference affects the computation of p(z).)

As in step 10., print a copy of the plot, this time to a file called "pressure_3.ps". You will again get extra credit for labeling it.

Return to Meteorology 301 Lab List

Return to Meteorology 301 Homepage