Home > atmlab > geophysics > pt2z.m

pt2z

PURPOSE ^

PT2Z Hydrostatic altitudes

SYNOPSIS ^

function z = pt2z(p,t,h2o,p0,z0,varargin)

DESCRIPTION ^

 PT2Z   Hydrostatic altitudes

    Calculates altitudes fulfilling hydrostatic equilibrium, based on
    vertical profiles of pressure, temperature and water vapour. Pressure
    and altitude of a reference point must be specified.

    Molecular weights and gravitational constants are hard coded and 
    function is only valid for the Earth.

    As the gravitation changes with altitude, an iterative process is
    needed. The accuracy can be controlled by *z_acc*. The calculations
    are repeated until the max change of the altitudes is below *z_acc*. If
    z_acc<0, the calculations are run twice, which should give an accuracy
    better than 1 m.

 FORMAT   z = pt2z( p, t, h2o, p0, z0 [,lat,z_acc,refell] )
        
 OUT   z         Altitudes [m].
 IN    p         Column vector of pressures [Pa].
       t         Column vector of temperatures [K].
       h2o       Water vapour [VMR]. Vector or a scalar, e.g. 0. 
       p0        Pressure of reference point [Pa].
       z0        Altitude of reference point [m].
       lat       Latitude. Default is 45.
       z_acc     Accuracy for z. Default is -1.
       ellipsoid Reference ellipsoid data, see *ellipsoidmodels*. 
                 Default is data matching WGS84.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

pt2z.m

SOURCE CODE ^

0001 % PT2Z   Hydrostatic altitudes
0002 %
0003 %    Calculates altitudes fulfilling hydrostatic equilibrium, based on
0004 %    vertical profiles of pressure, temperature and water vapour. Pressure
0005 %    and altitude of a reference point must be specified.
0006 %
0007 %    Molecular weights and gravitational constants are hard coded and
0008 %    function is only valid for the Earth.
0009 %
0010 %    As the gravitation changes with altitude, an iterative process is
0011 %    needed. The accuracy can be controlled by *z_acc*. The calculations
0012 %    are repeated until the max change of the altitudes is below *z_acc*. If
0013 %    z_acc<0, the calculations are run twice, which should give an accuracy
0014 %    better than 1 m.
0015 %
0016 % FORMAT   z = pt2z( p, t, h2o, p0, z0 [,lat,z_acc,refell] )
0017 %
0018 % OUT   z         Altitudes [m].
0019 % IN    p         Column vector of pressures [Pa].
0020 %       t         Column vector of temperatures [K].
0021 %       h2o       Water vapour [VMR]. Vector or a scalar, e.g. 0.
0022 %       p0        Pressure of reference point [Pa].
0023 %       z0        Altitude of reference point [m].
0024 %       lat       Latitude. Default is 45.
0025 %       z_acc     Accuracy for z. Default is -1.
0026 %       ellipsoid Reference ellipsoid data, see *ellipsoidmodels*.
0027 %                 Default is data matching WGS84.
0028 
0029 % 2005-05-11   Created by Patrick Eriksson.
0030 
0031 
0032 function z = pt2z(p,t,h2o,p0,z0,varargin)
0033 %
0034 [lat,z_acc,ellipsoid] = optargs( varargin, { 45, -1, NaN } );
0035 %
0036 if isnan(ellipsoid)
0037   ellipsoid = ellipsoidmodels('wgs84');
0038 end
0039                                                                             %&%
0040 rqre_nargin( 5, nargin );                                                   %&%
0041 rqre_datatype( p, @istensor1 );                                             %&%
0042 rqre_datatype( t, @istensor1 );                                             %&%
0043 rqre_datatype( h2o, @istensor1 );                                           %&%
0044 rqre_datatype( p0, @istensor0 );                                            %&%
0045 rqre_datatype( z0, @istensor0 );                                            %&%
0046 rqre_datatype( lat, @istensor0 );                                           %&%
0047 np = length( p );
0048 if length(t) ~= np                                                          %&%
0049   error('The length of *p* and *t* must be identical.');                    %&%
0050 end                                                                         %&%
0051 if ~( length(h2o) == np  |  length(h2o) == 1 )                              %&%
0052   error('The length of *h2o* must be 1 or match *p*.');                     %&%
0053 end                                                                         %&%
0054 if p0 > p(1)  |  p0 < p(np)                                                 %&%
0055   error('Reference point (p0) can not be outside range of *p*.');           %&%
0056 end                                                                         %&%
0057 
0058 
0059 %= Expand *h2o* if necessary
0060 %
0061 if  length(h2o) == 1
0062   h2o = repmat( h2o, np, 1 );
0063 end
0064 
0065 
0066 %= Make rough estimate of *z*
0067 %
0068 z = p2z_simple( p );
0069 z = shift2refpoint( p, z, p0, z0 );
0070 
0071 
0072 %= Set Earth radius and g at z=0
0073 %
0074 re = ellipsoidradii( ellipsoid, lat );
0075 g0 = lat2g0( lat );
0076 
0077 
0078 %= Gas constant and molecular weight of dry air and water vapour
0079 %
0080 r  = constants( 'GAS_CONST' );
0081 md = 28.966;
0082 mw = 18.016;
0083 %
0084 k  = 1-mw/md;        % 1 - eps
0085 rd = 1e3 * r / md;   % Gas constant for 1 kg dry air
0086 
0087 
0088 %= How to end iterations
0089 %
0090 if z_acc < 0
0091   niter = 2;
0092 else
0093   niter = 99;
0094 end  
0095 
0096 for iter = 1:niter
0097 
0098   zold = z;
0099   
0100   g = z2g( re, g0, z );
0101 
0102   for i = 1 : (np-1)
0103       
0104     gp  = ( g(i) + g(i+1) ) / 2;
0105   
0106     %-- Calculate average water VMR (= average e/p)
0107     hm  = (h2o(i)+h2o(i+1)) / 2;
0108   
0109     %--  The virtual temperature (no liquid water)
0110     tv = (t(i)+t(i+1)) / ( 2 * (1-hm*k) );   % E.g. 3.16 in Wallace&Hobbs
0111   
0112     %-- The change in vertical altitude from i to i+1
0113     dz = rd * (tv/gp) * log( p(i)/p(i+1) );
0114     z(i+1) = z(i) + dz;
0115       
0116   end
0117   
0118   %-- Match the altitude of the reference point
0119   z = shift2refpoint( p, z, p0, z0 );
0120 
0121   if z_acc >= 0 & max(abs(z-zold)) < z_acc
0122     break;
0123   end
0124   
0125 end
0126 
0127 return
0128 %----------------------------------------------------------------------------
0129 
0130 function z = shift2refpoint( p, z, p0, z0 )
0131   %
0132   z = z - ( interpp( p, z, p0 ) - z0 );
0133   %
0134 return
0135 
0136 
0137 function g = z2g(r_geoid,g0,z)
0138   %
0139   g = g0 * (r_geoid./(r_geoid+z)).^2;
0140   %
0141 return
0142 
0143 % Expression below taken from Wikipedia page "Gravity of Earth", that is stated
0144 % to be: International Gravity Formula 1967, the 1967 Geodetic Reference System
0145 % Formula, Helmert's equation or Clairault's formula.
0146 function g0 = lat2g0(lat)
0147   %
0148   x  = abs( lat );
0149   g0 = 9.780327 * ( 1 + 5.3024e-3*sind(x).^2 + 5.8e-6*sind(2*x).^2 );
0150   %
0151 return
0152

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005