Home > atmlab > h2o > parametrisations > dBZ2iwcMH97.m

dBZ2iwcMH97

PURPOSE ^

dBZ2iwcMH97 returns iwc from radar reflectivity and temperature

SYNOPSIS ^

function [IWC] = dBZ2iwcMH97(dBZ,T,f,clo)

DESCRIPTION ^

 dBZ2iwcMH97 returns iwc from radar reflectivity and temperature
             while assuming McFarquhar and Heymsfield (1997)
             particle size dsitribution
                          

 OUT      IWC is a matrix with iwc [g/m^3] profiles

 IN       dBZ [mm^6/m^-3] matrix 
          outside this range nan is reported

          T [K] matrix [limits: 200 - 272]
          outside this range nan is reported

          f frequency (scalar)

          clo.do create look-up-table (0 or 1) 
          clo.file filename where the look-up-table
                   is/will be stored          

 
          comments: IWC values are linearly interpolated in 
                    a look-up-table of IWC(dBZ,T).
                    By advantage this look-up-table
                    is precalculated and stored
                    as a matlab file. The file is not large
                    but take some time to calculate.
          Example usage:
              tmpfolder='/home/bengt/tmp';
              clo.do=0;
              clo.file=[tmpfolder,'/dBZ-data.mat']
              dBZ=[-30 -40;-20 16];
              T=[230 240;220 220];
              [IWC] = dBZ2iwcMH97(dBZ,T,clo)


 FORMAT   [IWC] = dBZ2iwcMH97(dBZ,T,clo)


 History: 2007-03-27  Created by Bengt Rydberg

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

dBZ2iwcMH97.m

SOURCE CODE ^

0001 % dBZ2iwcMH97 returns iwc from radar reflectivity and temperature
0002 %             while assuming McFarquhar and Heymsfield (1997)
0003 %             particle size dsitribution
0004 %
0005 %
0006 % OUT      IWC is a matrix with iwc [g/m^3] profiles
0007 %
0008 % IN       dBZ [mm^6/m^-3] matrix
0009 %          outside this range nan is reported
0010 %
0011 %          T [K] matrix [limits: 200 - 272]
0012 %          outside this range nan is reported
0013 %
0014 %          f frequency (scalar)
0015 %
0016 %          clo.do create look-up-table (0 or 1)
0017 %          clo.file filename where the look-up-table
0018 %                   is/will be stored
0019 %
0020 %
0021 %          comments: IWC values are linearly interpolated in
0022 %                    a look-up-table of IWC(dBZ,T).
0023 %                    By advantage this look-up-table
0024 %                    is precalculated and stored
0025 %                    as a matlab file. The file is not large
0026 %                    but take some time to calculate.
0027 %          Example usage:
0028 %              tmpfolder='/home/bengt/tmp';
0029 %              clo.do=0;
0030 %              clo.file=[tmpfolder,'/dBZ-data.mat']
0031 %              dBZ=[-30 -40;-20 16];
0032 %              T=[230 240;220 220];
0033 %              [IWC] = dBZ2iwcMH97(dBZ,T,clo)
0034 %
0035 %
0036 % FORMAT   [IWC] = dBZ2iwcMH97(dBZ,T,clo)
0037 %
0038 %
0039 % History: 2007-03-27  Created by Bengt Rydberg
0040 %
0041 
0042 function [IWC] = dBZ2iwcMH97(dBZ,T,f,clo)
0043 
0044 if  sum(size(T)==size(dBZ))~=2
0045    error('size of T must be equal size of dBZ!!!')
0046 end
0047 if clo.do~=0 & clo.do~=1
0048    error('clo.do must be 0 or 1!!!')
0049 end
0050 if ~isstr(clo.file)
0051     error('clo.file must be a string specifying a filename!!!')
0052 end    
0053 if  clo.do==0
0054     if exist(clo.file,'file')~=2 
0055         error(['the file ',clo.file,' does not exist!!!'])
0056     end
0057 end
0058 
0059 if clo.do
0060     
0061   %calculate the radar backscattering cross-section for spheres
0062   %particle diameter
0063   D=5e-6:10e-6:1e-2;
0064   x_i=D/2;
0065   F=create_ssp('sphere',f,x_i/2);
0066 
0067   %loop over particle diameters
0068   for i=1:length(x_i)
0069       %radar backscattering per particle
0070       Q_b1(:,i)=4*pi*squeeze(F(i).pha_mat_data(1,:,end,1,1,1,1))';
0071       %extinction per particle
0072       Q_b2(:,i)=F(i).ext_mat_data;
0073   end
0074   
0075   %create the look-up-table
0076   c=3e8;
0077   lambda=c/f;
0078   nwater=sqrt(eps_water_liebe93(f,273.15));
0079   Kwater=( abs( (nwater.^2-1)./ (nwater.^2+2) ) ).^2;
0080   %
0081   IWCv = logspace(-3.65,0.58,100);
0082   Tv   = 190:2:272;
0083   mode = 1;
0084   %
0085   dx=D(2)-D(1);
0086   T_grid=F(1).T_grid;
0087   %loop over Tv and IWCv to find dBZ(T,IWC)
0088   clear y Z ext
0089   for j=1:length(Tv)
0090       sigma_b=interp1(T_grid,Q_b1,Tv(j));
0091       sigma_c=interp1(T_grid,Q_b2,Tv(j));
0092       for i=1:length(IWCv)
0093           y(i,:) = ice_psd_Mcfar_97(Tv(j),IWCv(i),D,mode);
0094       end
0095       Z(j,:)=lambda^4 /(pi^5)/Kwater*...
0096              sum([y.*[ones(length(IWCv),1)*sigma_b]*dx]');
0097       ext(j,:)=sum([y.*[ones(length(IWCv),1)*sigma_c]*dx]');
0098   end
0099   dBZm=10*log10(Z*1e18);
0100 
0101   %plot a figure showing the relationships
0102   if 0
0103      semilogy(dBZm,IWCv)
0104      legend(num2str(Tv'))
0105   end
0106 
0107   %re-arrange the look-up-table [IWC(T,dBZ)]
0108   dBZv=[-40:1:15];
0109   for i=1:size(dBZm,1)
0110       IWCC(i,:)=interp1(dBZm(i,:),IWCv,dBZv);
0111   end
0112   P.IWC=IWCC;
0113   P.T=Tv;
0114   P.dBZ=dBZv; 
0115   
0116   save(clo.file,'P')
0117 end
0118 
0119 load(clo.file)
0120 
0121 ind1=find(isnan(dBZ(:)));
0122 ind2=find(isnan(T(:)));
0123 ind=union(ind1,ind2);
0124 ind=setdiff(1:length(dBZ(:)),ind); 
0125 
0126 
0127 IWC=zeros(size(dBZ));
0128 IWC(ind)=interp2(P.dBZ,P.T,P.IWC,dBZ(ind),T(ind));
0129 
0130 
0131 %- Handle data outisde covered dBZ-T area
0132 %
0133 % Low dBZ and high/low T -> IWC = 0
0134 ind      = find( dBZ<min(P.dBZ) |  T>max(P.T) |  T<min(P.T) );
0135 IWC(ind) = 0;
0136 %
0137 % High dBZ -> max(IWC)
0138 ind      = find( dBZ>max(P.dBZ) );
0139 IWC(ind) = mean( max( P.IWC' ) );
0140 
0141 
0142 if any(isnan(IWC))
0143   fprintf( 'NaN found. Check cause!\n' );
0144   keyboard
0145 end

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