Meteorology 301 - Lab 6

Introduction

In this lab, you will build on work done in Lab 5. In the process, you will in effect write a computer program. You have already taken major steps in this direction in previous Matlab-based labs. You will continue here with some further computations of atmospheric water vapor.

First, you will add to your set of Matlab functions ("m" files) by copying two new ones to your Matlab space. Then you will assemble a program that calls these functions in sequence. Finally, you will write a function of your own to compute precipitable water.

This lab is in itself a step toward later work, in which we will use results of this lab as part of a simple model of the Earth's climate based on energy balance.

GOAL: Develop a program that will compute and plot precipitable water (PW) as a function of assumed temperature and relative humidity

Instructions

1. As with previous Matlab labs, follow the instructions for accessing meteorology computers to log into your Vincent account. Each student should open a window on the host machine, connect to Vincent in each window, and then take turns completing the rest of the lab.

2. I have written three new function files for you to use:

  file                  output
--------    ----------------------------------------
relhum.m    relative humidity as function of z
WSAT.m	    saturation mixing ratio as function of z
Main.m      a function that calls other function files

Your should copy them to your home directory or where you keep your Matlab work if it is not your home directory. (The commands below assume your Matlab work is in your home directory. Change them appropriately if it is elsewhere.)

Copy the file relhum.m to your home directory
cp /mthome/gutowski/MT301/relhum.m ~
(That last character in the line is the tilde, ~. Note that there is a space before the "~".)

Copy the file WSAT.m to your home directory
cp /mthome/gutowski/MT301/WSAT.m ~

Copy the file Main.m to your home directory
cp /mthome/gutowski/MT301/Main.m ~

3. Take a look at relhum.m

emacs relhum.m

Answer these questions

What input variables does the function require?
Why is negative relative humidity undesireable?
How do we prevent negative relative humidity in this function file?

4. Take a look at WSAT.m

emacs WSAT.m

Answer these questions

What input variables does the function require?
What items are computed in the loop over atmospheric layers?
What equation is used to compute saturation vapor pressure?

5. Now start Matlab, from a Vincent prompt. If you forget how, go to Lab 5 , step 4.

6. Let's see what these new functions give us. Enter these commands from within Matlab:

First, we need to set up the input variables. For our lapse rate, we will use the standard atmosphere; for our surface temperature, we will use the global average values.
>>N=51
>>Dz=200
>>Ts=288
>>gamma=6.5e-3
The relative-humidity computation needs some background information - its surface value and the rate it changes with height. Guided by the observed relative humidity vs. height, enter
>>RHs=0.80
>>dRHdz=-0.4/10000
(Note the minus sign and that the denominator in the last line is ten thousand.)
EXTRA CREDIT - Answer this question
1. Looking at the observed relative humidity vs. height , why are these choices of RHs and dRHdz reasonable? (The figure in the link shows relative humidity, called [u] in the figure, as a function of latitude and pressure. Pressure in the figure is measured in units of decibars, so the figure goes from 1000 mb to 200 mb.)
Compute relative humidity and specific humidity as functions of height, using the Main.m function:
>>[z,TEM,PR,RH,q]=Main(N,Dz,Ts,RHs,gamma,dRHdz)
Note that the output quantities are all contained within the square brackets, [ ].
Plot humidity quantities as functions of z:
>>plot(RH,z)
>>plot(q,z)
Answer this question
1. Why does q decrease so rapidly with height? (Think about T(z) - you may want to plot TEM vs. z to help you.)
Compute specific humidity, q(z), for other surface temperatures. What approximate values of Ts gives q(3km)=0.002 and q(3km)=0.005?
Answer this question
1. How does q(z) change the way it does when surface temperature changes? Why?

7. Now write your own function to compute precipitable water (PW) for a given q profile.

Recall that PW is the vertically integrated water vapor in the atmosphere.

          infinity            psfc
          /                   /
          |                   |
     PW = | q*(density) dz =  | q dp / g
          |                   |
         /                   /
         0                   0

where psfc is surface pressure.

We will assume that the computer can not compute integrals directly (generally true), but instead approximates them as sums. How will we do this? Consider the plot below showing q vs. atmospheric level:

pressure 
 level
   .
   .
   .
           |
   6      --  *
           |
   5      -    *
           |    
   4      --     *
           |
   3      --          *
           |
   2      --                *
           | 
   1      -----------------------*------------
                          q

It's easier to approximate the second integral form above, since g can be treated as a constant, but (density) can not. In this case, we can the write

  PW ~ - (1/g) *  SUM{ 0.5*( q(i)+q(i-1) ) * ( PR(i)-PR(i-1) )}
  
                                           , for i = 2, ..., N

where the SUM is over all level pairs (e.g., 1&2, 2&3, 3&4, ..., (N-1)&N ).

Why did we make this choice? We can make the approximation

 dp ~ (PR(i) - PR(i-1)).

Then we want the mean q in the layer between levels i and i-1:

 mean(q) ~ ( q(i) + q(i-1) )/2

We then sum over all the individual layers above. If the layer structure is fine enough, then these approximations are reasonably accurate.

Your task is to write such a function. You will probably find it useful to look at other functions you have used for help on how to write it, especially for writing the loop over atmospheric levels. Call the function prec_wat.m

I will help you out by giving you the function identification and comment lines for such a function. Comments of the form "%(loop)" are for line codes within the loop:

   function [PW] = prec_wat(N,PR,q)

   % Function to compute the vertical integral for precipitable water 
   % as a summation over atmospheric layers bounded by the 
   % pressure levels PR

   % First, intialize g to g = 9.8

   % Set a summation variable SUM equal to 0

   % Loop over all layers, starting from the layer bounded by pressure
   % levels 1 and 2 and continuing up to the combination N-1 and N

   %(loop)  Compute mean q for the layer

   %(loop)  Compute pressure difference across the layer.

   %(loop) Convert pressure difference from mb to Pascals, the MKS unit
   %(loop) (i.e., multiply by 100)

   %(loop)  Compute the product and sum:  SUM = SUM + qmean*dp/g

   %(loop)  End of loop

   %  Equate PW with the final value of SUM after looping

To test your function, go back to Matlab and re-issue the commands
 
>>N=51

>>Dz=200

>>Ts=288

>>gamma=6.5e-3

>>RHs=0.80

>>dRHdz=-0.4/10000

Then enter

>>[z,TEM,PR,RH,q]=Main(N,Dz,Ts,RHs,gamma,dRHdz)
 
>>PW=prec_wat(N,PR,q)



You should get an answer back of 20.8936, or something very close to
that.

8. You should turn in to me the answers to your questions along with a
printed copy of your function file prec_wat.m  

You should also email me a copy of prec_wat.m


Return to 
Meteorology 301 Lab List

Return to 
Meteorology 301 Homepage