Home > atmlab > geophysics > twvcalc.m

twvcalc

PURPOSE ^

TWVCALC Total column water vapor

SYNOPSIS ^

function twv = twvcalc(p, t, z, h2o, plogspacing, verbose)

DESCRIPTION ^

 TWVCALC Total column water vapor

 Calculate the column water vapor by vertically integrating the
 water vapor density. 
 
 For better accuarcy, first interpolate onto a fine vertical
 grid. The interpolation is linear in log(p). For humidity, the
 interpolation is done for RH, not VMR.
 
 FORMAT twv = twvcalc(p, t, z, h2o, plogspacing, verbose)

 OUT    y     Total column water vapor in kg/m^2.
 IN     p     Vector of pressures [Pa]. Must be strictly decreasing.
        t     Vector of temperatures [K].
        h2o   Vater vapor [VMR]. 
        plogspacing Max spacing (in log units) of the pressure grid on
                    which we interpolate. Must not be too coarse,
                    recommended value: 0.01
        verbose     Output some information if this is present and
                    not zero.

 2007-08-06 Created by Stefan Buehler.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

twvcalc.m

SOURCE CODE ^

0001 % TWVCALC Total column water vapor
0002 %
0003 % Calculate the column water vapor by vertically integrating the
0004 % water vapor density.
0005 %
0006 % For better accuarcy, first interpolate onto a fine vertical
0007 % grid. The interpolation is linear in log(p). For humidity, the
0008 % interpolation is done for RH, not VMR.
0009 %
0010 % FORMAT twv = twvcalc(p, t, z, h2o, plogspacing, verbose)
0011 %
0012 % OUT    y     Total column water vapor in kg/m^2.
0013 % IN     p     Vector of pressures [Pa]. Must be strictly decreasing.
0014 %        t     Vector of temperatures [K].
0015 %        h2o   Vater vapor [VMR].
0016 %        plogspacing Max spacing (in log units) of the pressure grid on
0017 %                    which we interpolate. Must not be too coarse,
0018 %                    recommended value: 0.01
0019 %        verbose     Output some information if this is present and
0020 %                    not zero.
0021 %
0022 % 2007-08-06 Created by Stefan Buehler.
0023 
0024 function twv = twvcalc(p, t, z, h2o, plogspacing, verbose)
0025 
0026 %= Check input
0027 
0028 rqre_nargin( 5, nargin );
0029 
0030 rqre_datatype( p, {@istensor1} );                                           %&%
0031 rqre_datatype( t, {@istensor1} );                                           %&%
0032 rqre_datatype( z, {@istensor1} );                                           %&%
0033 rqre_datatype( h2o, {@istensor1} );                                         %&%
0034 rqre_datatype( plogspacing, {@istensor0} );                                 %&%
0035 
0036 np = length( p );
0037 
0038 log_p = log(p);
0039 
0040 mindiff = min(abs(diff(log_p)));
0041 
0042 if max(diff(p)) >= 0
0043   error('The pressure grid *p* must be strictly decreasing.');
0044 end 
0045 
0046 if length(t) ~= np
0047   error('The length of *p* and *t* must be identical.');
0048 end 
0049 
0050 if length(z) ~= np
0051   error('The length of *p* and *z* must be identical.');
0052 end 
0053 
0054 if min(diff(z)) <= 0
0055   error('The vector *z* must be strictly increasing.');
0056 end 
0057 
0058 if length(h2o) ~= np
0059   error('The length of *p* and *h2o* must be identical.');
0060 end 
0061 
0062 if plogspacing <= 0
0063   error('The value of *plogspacing* must be greater than zero.');
0064 end 
0065 
0066 if plogspacing > mindiff
0067   error(['The value of *plogspacing* is larger than the original log\n' ...
0068          'grid spacing, which is %g'], mindiff);
0069 end 
0070 
0071 
0072 %= Start to do some business now.
0073 
0074 % Create fine log pressure grid:
0075 
0076 % Construct a fine pressure grid. We go from max pressure to min pressure,
0077 % in steps given by the desired spacing.
0078 log_p_fine = log_p(1):-plogspacing:log_p(end);
0079 
0080 % Interpolate T, z, RH:
0081   
0082 t_int = interp1( log_p, t, log_p_fine ); 
0083 z_int = interp1( log_p, z, log_p_fine ); 
0084 
0085 % RH (RH = p*VMR_H2O/e_eq_water(T)):
0086 rh = p .* h2o ./ e_eq_water(t);
0087 
0088 rh_int = interp1( log_p, rh, log_p_fine ); 
0089 
0090 % H2O In VMR units:
0091 h2o_int = rh_int .* e_eq_water(t_int) ./ exp(log_p_fine); 
0092 
0093 
0094 % The gas constant of water vapor in J/(K kg), from Wallace&Hobbs,
0095 % 2nd edition:
0096 Rv = 461;
0097 
0098 % The density is
0099 % rho = p*h2o / (Rv*T)
0100 % (from W&H)
0101 
0102 rho = exp(log_p_fine) .* h2o_int ./ (Rv * t_int);
0103 
0104 
0105 % Now vertically integrate rho:
0106 twv = sum(layermean(rho) .* diff(z_int));
0107 
0108 % Output some information, if verbose is set:
0109 if exist('verbose','var') 
0110   if verbose ~= 0
0111     disp(sprintf('Original log_p spacing: %g',mindiff));
0112     disp(sprintf('No of old/new pressure grid points: %d/%d', ...
0113                  length(log_p), ...
0114                  length(log_p_fine)));
0115     disp(sprintf('Min and max RH: %g, %g', min(rh), max(rh)));
0116   end
0117 end
0118 
0119 
0120 % That's it, we are done.

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