Home > atmlab > h2o > thermodynamics > pseudo_adiabatic_liwc.m

pseudo_adiabatic_liwc

PURPOSE ^

PSEUDO_ADIABATIC_LIWC Calculate pseudo-adabatic liquid

SYNOPSIS ^

function [T1,xwc]=pseudo_adiabatic_liwc(T0,P0,P1,flag)

DESCRIPTION ^

PSEUDO_ADIABATIC_LIWC   Calculate pseudo-adabatic liquid
                       water content and temperature

          Calculates pseudo-adiabatic lwc or iwc
          content and pseudo-adiabatic 
          temperature following eq. 2.34
          in R.R. Rodgers, "A short course
          in cloud physics"  
    
 FORMAT   [T1,xwc]=pseudo_adiabatic_liwc(T0,P0,P1,flag)

 IN    T0  temperature at   (scalar) [K]
           cloudbase
       P0  pressure at      (scalar) [Pa]
           cloudbase
       P1  pressure levels  (vector) [Pa]
           above cloud-base
       fLag tell whether to calculate the icw (1) or lwc (0)

 OUT   T1  temperature at   (vector) [K] 
           P1 levels
       lwc adiabatic liquid (vector) [kg/m3]
           water content at
           P1 levels
       iwc adiabatic ice (vector) [kg/m3]
           water content at
           P1 levels

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

pseudo_adiabatic_liwc.m

SOURCE CODE ^

0001 %PSEUDO_ADIABATIC_LIWC   Calculate pseudo-adabatic liquid
0002 %                       water content and temperature
0003 %
0004 %          Calculates pseudo-adiabatic lwc or iwc
0005 %          content and pseudo-adiabatic
0006 %          temperature following eq. 2.34
0007 %          in R.R. Rodgers, "A short course
0008 %          in cloud physics"
0009 %
0010 % FORMAT   [T1,xwc]=pseudo_adiabatic_liwc(T0,P0,P1,flag)
0011 %
0012 % IN    T0  temperature at   (scalar) [K]
0013 %           cloudbase
0014 %       P0  pressure at      (scalar) [Pa]
0015 %           cloudbase
0016 %       P1  pressure levels  (vector) [Pa]
0017 %           above cloud-base
0018 %       fLag tell whether to calculate the icw (1) or lwc (0)
0019 %
0020 % OUT   T1  temperature at   (vector) [K]
0021 %           P1 levels
0022 %       lwc adiabatic liquid (vector) [kg/m3]
0023 %           water content at
0024 %           P1 levels
0025 %       iwc adiabatic ice (vector) [kg/m3]
0026 %           water content at
0027 %           P1 levels
0028 
0029 % example usage:
0030 %               T0=273;P0=87e3;
0031 %               P1=logspace(log10(870e2),log10(300e2),20);
0032 %               [T1,xwc]=pseudo_adiabatic_liwc(T0,P0,P1,flag)
0033 
0034 % 2009-03-30 Created by Bengt Rydberg, Marston Johnston
0035 
0036 % Assumptions:
0037 %   1. The release of latent heat when supercooled water freezes onto
0038 %       ice is ignored.
0039 %   2. The process where by the liquid water evaporates is not considered
0040 %   3. The exact point at where ice is formed is assumed to exist below 0
0041 %       degrees. In the atmosphere this may or may not be true
0042 
0043 
0044 function [T1,xwc]=pseudo_adiabatic_liwc(T0,P0,P1,flag)
0045 
0046 %some checks
0047 if length(T0)~=1
0048  error('T0 must be a scalar')
0049 end
0050 if length(P0)~=1
0051  error('P0 must be a scalar')
0052 end
0053 if ~any(size(P1)==1)
0054  error('P1 must be a vector')
0055 end
0056 
0057 Rp  = 287;               % individual gas constant for dry air
0058 cp  = 1005;             % specific heat capacity at constant pressure
0059 k   = Rp/cp;
0060 ep  = 0.622; 
0061 Rv  = 461.5;            % individual gas constant for water vapor
0062 Lv  = 2.5*1e6;     % Latent heat of condensation/evaporation @ 0 deg
0063 Ls  = 2.85*1e6;       % Latent heat of sublimation/deposition
0064 Lf  = 3.34*1e5;       % Latent heat of fusion/melting
0065 dp  = -10;                % small pressure change
0066 P   = P0:dp:min(P1)+dp; % pressure grid
0067 T   = zeros(size(P));   % temperature vector
0068 T(1)=T0;
0069 
0070 for i=1:length(P)
0071 
0072  % saturation water vapor pressure
0073  if flag == 0 
0074      es    =  e_eq_water(T(i));
0075  elseif flag == 1
0076      es    = e_eq_ice(T(i));
0077  else
0078      error('Variable for LWC/IWC flag not set');
0079  end
0080      
0081  % saturation mixing ratio
0082  ws(i) = ep*es/(P(i)-es);
0083  % ambient temperature change
0084  if T(i) > 273 
0085      L = Lv; 
0086  else
0087      L = Ls; 
0088  end
0089 
0090  dT    = (k/P(i)+2*L*ws/T(i)/cp/(P(i)-es))/...
0091          (1/T(i)+L^2*ws*P(i)/(Rv*T(i)^3*cp*(P(i)-es)))*dp;
0092  T(i+1) = T(i) + dT;
0093 end
0094 
0095 % weight of humid air
0096 rho_air = (P - e_eq_water(T(1:end-1)))./Rp./T(1:end-1)+...
0097            e_eq_water(T(1:end-1))/Rv./T(1:end-1);
0098 
0099 % adiabatic liquid/ice water content
0100 
0101 awc = rho_air(1)*ws(1)-rho_air.*ws;
0102 if flag == 1 % too warm for ice development
0103     awc(T>273) = 0;
0104 end
0105 
0106 % interpolate temperature and adaiabatic liquid water content
0107 % on desired output levels
0108 xwc     = interp1(log(P),awc,log(P1));
0109 T1      = interp1(log(P),T(1:end-1),log(P1));

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