Re: Question of theta-e

Paul:

> Does anyone know if it is possible to determine the value of a wet adiabat
> on a thermodynamic diagram given the pressure and temperature. I would lke
> to interrogate a sounding, determine the value of the wet adiabat that
> passes through each point, thereby determining what theta-e ought to
> approach in order to break the observed cap.

Many years ago, we started doing this by an iterative process that produces
excellent results.  Below are a bunch of static methods taken from a Java class
that deals with sounding and some thermodynamic computations.  The last method
will show you the technique.  Hope this helps.

tom



  /** Virtual Temperature
  *
  */
  public static double virtualTemperature(double t, double td, double p) {
    return ( t * (1.0 + .000609*mixingRatio(td, p)) );
  }

  /** saturation vapor pressure over water.  t in kelvin.
  *
  */
  public static double satVapPres(double t) {
    double coef[]={6.1104546,0.4442351,1.4302099e-2, 2.6454708e-4,
              3.0357098e-6, 2.0972268e-8, 6.0487594e-11,-1.469687e-13};

    // sat vap pressures every 5C from -50 to -200
    double escold[] = {
      0.648554685769663908E-01, 0.378319512256073479E-01,
      0.222444934288790197E-01, 0.131828928424683120E-01,
      0.787402077141244848E-02, 0.473973049488473318E-02,
      0.287512035504357928E-02, 0.175743037675810294E-02,
      0.108241739518850975E-02, 0.671708939185605941E-03,
      0.419964702632039404E-03, 0.264524363863469876E-03,
      0.167847963736813220E-03, 0.107285397631620379E-03,
      0.690742634496135612E-04, 0.447940489768084267E-04,
      0.292570419563937303E-04, 0.192452912634994161E-04,
      0.127491372410747951E-04, 0.850507010275505138E-05,
      0.571340025334971129E-05, 0.386465029673876238E-05,
      0.263210971965005286E-05, 0.180491072930570428E-05,
      0.124607850555816049E-05, 0.866070571346870824E-06,
      0.605982217668895538E-06, 0.426821197943242768E-06,
      0.302616508514379476E-06, 0.215963854234913987E-06,
      0.155128954578336869E-06};

    double temp = t - 273.16;
    double retval;

    if (temp > -50.) {

      retval = ( coef[0] + temp*(coef[1] + temp*(coef[2] + temp*(coef[3] +
      temp*(coef[4] + temp*(coef[5] + temp*(coef[6] + temp*coef[7])))))) );

    } else {
       double tt = (-temp - 50.)/5.;
       int inx = (int) tt;
       if (inx < escold.length) { 
         retval = escold[inx] + (tt % 1.)*(escold[inx+1]-escold[inx]);
       } else {
         retval = 1e-7;
       }

    }
    
    return retval;
  }


  /** mixing ratio
  *
  */
  public static double mixingRatio(double t, double p) {
    double e = satVapPres(t);
    return ( 621.97*e/(p - e) );
  }

  /** temperature at LCL (old SELS method)
  *
  */
  public static double tempAtLCL(double t, double td) {
    return (td - (.001296*td - .15772)*(t-td) );
  }

  /** define the parcel's theta-E.  Pressure should be Station Pressure!
  *
  */
  public static double parcelThetaE(double t, double td, double p) {
    double theta = t * Math.pow ( 1000.0/p, .286);
    return (theta * Math.exp(2481.9e-3*mixingRatio(td,p)/tempAtLCL(t, td)) );
  }

  /** pressure at the LCL
  *
  */
  public static double pressureAtLCL(double t, double td, double p) {
    double theta = t * Math.pow( 1000.0/p, .286);
    return (1000./ Math.pow( theta/tempAtLCL(t, td), 3.5) );
  }

  /** compute temperature along Pseudo Adiabats.  This iterates
  * until it converges....
  *
  */
  public static double tempAlongSatAdiabat(double thetaE, double p) {
    double s, th;

    double pcon = Math.pow( 1000.0/p, .286);
    double t = 273.0;
    double delta = 20.;
    do {
      s = mixingRatio(t, p);
      th = t * pcon * Math.exp(2.5*s/t);
      if ( (th-thetaE)*delta > 0.0) delta = -.5 * delta;
      t = t + delta;

    } while (Math.abs(delta) > .1);
    
    return (t);
  }

}




-- 
Tom Whittaker (tomw@xxxxxxxxxxxxx)
University of Wisconsin-Madison
Space Science and Engineering Center
Phone/VoiceMail: 608/262-2759
Fax: 608/262-5974