Home > atmlab > arts > scenegen > asg_dbz2pnd.m

asg_dbz2pnd

PURPOSE ^

asg_dbz2pnd calculates iwc and pnd fields

SYNOPSIS ^

function [G]= asg_dbz2pnd( G, Q ,P);

DESCRIPTION ^

asg_dbz2pnd calculates iwc and pnd fields

 FORMAT   G=asg_dbz2pnd( G, Q, P)
 
 OUT      G        Modified gformat data,now including
                   IWC and PND data
 IN     
          G        original gformat data 
          Q        Qarts setting structure.            
          P        conversion structure of length 2
                   P(1) corresponds to the radar frequency                    
                   P(2) corresponds to the instrument frequencys 
                   with fields                  
            SSP    single scattering properties    
            PSD    particle size dist. choice
                   'MH97' is the only valid option
            method method for calculating particle
                   number density fields
                   'gauss-laguerre' is the only valid option
            x      particle diameter (see function gauss_laguerre
                   for help)                     
            w      weights
            x_norm normalisation factor 
            shape  particle shape
                   'sphere' is the only valid option 

  2007-11-13 created by Bengt Rydberg

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

asg_dbz2pnd.m

SOURCE CODE ^

0001 %asg_dbz2pnd calculates iwc and pnd fields
0002 %
0003 % FORMAT   G=asg_dbz2pnd( G, Q, P)
0004 %
0005 % OUT      G        Modified gformat data,now including
0006 %                   IWC and PND data
0007 % IN
0008 %          G        original gformat data
0009 %          Q        Qarts setting structure.
0010 %          P        conversion structure of length 2
0011 %                   P(1) corresponds to the radar frequency
0012 %                   P(2) corresponds to the instrument frequencys
0013 %                   with fields
0014 %            SSP    single scattering properties
0015 %            PSD    particle size dist. choice
0016 %                   'MH97' is the only valid option
0017 %            method method for calculating particle
0018 %                   number density fields
0019 %                   'gauss-laguerre' is the only valid option
0020 %            x      particle diameter (see function gauss_laguerre
0021 %                   for help)
0022 %            w      weights
0023 %            x_norm normalisation factor
0024 %            shape  particle shape
0025 %                   'sphere' is the only valid option
0026 %
0027 %  2007-11-13 created by Bengt Rydberg
0028 
0029 function [G]= asg_dbz2pnd(  G, Q ,P);
0030 
0031 r_ind=find(strncmp(lower({G(:).NAME}),'radar',5));
0032 t_ind=find(strncmp(lower({G(:).NAME}),'temperature',11));
0033 z_ind=find(strncmp(lower({G(:).NAME}),'altitude',8));
0034 
0035 dBZ=G(r_ind).DATA;
0036 T=G(t_ind).DATA;
0037 ALT=G(z_ind).DATA;
0038 
0039 warning off
0040 
0041 if length(P)~=2
0042    error('P structure must be of length 2')
0043 end
0044 
0045 if ~exist('T','var') & ~exist('dBZ','var')
0046    error('G must hold Temperature and Radar_Reflectivity data')
0047 end
0048 
0049 
0050 if strcmp(P(1).method,'gauss-laguerre')
0051    len=size(P(1).radar_back,1);
0052    if length(P(1).x)~=len & length(P(1).w)~=len
0053       error('mismatch in size between P(1).radar_back,P(1).x, and P(1).w')
0054    end
0055    Q_b1=P(1).radar_back';
0056    %create the look-up-table for iwc
0057    c=constants('SPEED_OF_LIGHT');
0058    f=P(1).F_GRID;
0059    lambda=c/f;
0060    nwater=sqrt(eps_water_liebe93(f,273.15));
0061    Kwater=( abs( (nwater.^2-1)./ (nwater.^2+2) ) ).^2;
0062    IWCv=logspace(-3.65,0.5,100);
0063    Tvi=180:2:273;
0064    T_grid=P(1).T_grid;
0065    %loop over Tv and IWCv to find dBZ(T,IWC)
0066    for j=1:length(Tvi)
0067        sigma_b=interp1(T_grid,Q_b1,Tvi(j));
0068        %sigma_c=interp1(T_grid,Q_b2,Tvi(j));
0069        for i=1:length(IWCv)
0070            if strcmp(P(1).PSD,'MH97')
0071               y(i,:) = ice_psd_Mcfar_97(Tvi(j),IWCv(i),P(1).x,1);
0072               %rho=0.91e6;
0073               %a=1;
0074               %N=6*IWCv(i).*((3.67+a)./dm).^(4+a)/(rho*pi*gamma(4+a));
0075               %y(i,:)=N*P(1).x.^a.*exp(-(3.67+a)*P(1).x/dm);
0076               Y(i,:) = gauss_laguerre_apply(y(i,:)',P(1).x,P(1).w,P(1).x_norm);
0077            else
0078              error('MH97 is the only valid option for P.PSD')
0079            end
0080        end
0081        Z(j,:)=lambda^4 /(pi^5)/Kwater*...
0082                            sum([Y.*[ones(length(IWCv),1)*sigma_b]]');
0083        %ext(j,:)=sum([Y.*[ones(length(IWCv),1)*sigma_c]]');
0084    end
0085    dBZm=10*log10(Z*1e18);
0086    
0087    len=size(P(2).radar_back,1);
0088    if length(P(2).x)~=len & length(P(2).w)~=len
0089       error('mismatch in size between P(2).SSP,P(2).x, and P(2).w')
0090    end
0091    clear Q_b1 Q_b2 y Y
0092    Q_b1=P(2).radar_back';
0093   
0094    %create the look-up-table for lwc
0095    LWCv=logspace(-3.65,2,100);
0096    %IWCv=logspace(-4,6,100);
0097    Tvw=273:2:313;
0098    T_grid=P(2).T_grid;
0099    %loop over Tv and IWCv to find dBZ(T,IWC)
0100    for j=1:length(Tvw)
0101        sigma_b=interp1(T_grid,Q_b1,Tvw(j));
0102        %sigma_c=interp1(T_grid,Q_b2,Tvw(j));
0103        for i=1:length(IWCv)
0104            if strcmp(P(2).PSD,'Water')
0105               c1=6;c2=1;rc=10;
0106               y(i,:)=water_psd(LWCv(i),P(2).x,rc,c1,c2)';
0107               %rho=0.91e6;
0108               %a=1;
0109               %N=6*IWCv(i).*((3.67+a)./dm).^(4+a)/(rho*pi*gamma(4+a));
0110               %y(i,:)=N*P(1).x.^a.*exp(-(3.67+a)*P(1).x/dm);
0111               Y(i,:) = gauss_laguerre_apply(y(i,:)',P(2).x,P(2).w,P(2).x_norm);
0112            else
0113              error('not a valid option for P.PSD')
0114            end
0115        end
0116        Zw(j,:)=lambda^4 /(pi^5)/Kwater*...
0117                            sum([Y.*[ones(length(IWCv),1)*sigma_b]]');
0118        %extw(j,:)=sum([Y.*[ones(length(IWCv),1)*sigma_c]]');
0119    end
0120    dBZmw=10*log10(Zw*1e18);
0121 
0122 
0123    %estimate the unattenuated dBZ field
0124    if 0
0125 
0126    dBZp=zeros(size(dBZ));
0127    for i=1:size(dBZ,2) 
0128     for j=1:size(dBZ,3)
0129      for k=size(dBZ,1):-1:1
0130        %ind of the closest temperature
0131        if T(k,i,j)<273
0132         ind=find(abs(T(k,i,j)-Tvi)==min(abs(T(k,i,j)-Tvi)));
0133         ice=1;
0134        else
0135         ind=find(abs(T(k,i,j)-Tvw)==min(abs(T(k,i,j)-Tvw)));
0136         ice=0;
0137        end
0138 
0139        if k==size(dBZ,1)
0140         dz=ALT(k,i,j)-ALT(k-1,i,j);
0141         if ice
0142          qext(k)=interp1(10.^(dBZm(ind,:)/10),ext(ind,:),...
0143                         10^(dBZ(k,i,j)/10),'linear','extrap')*dz;
0144         else
0145          qext(k)=interp1(10.^(dBZmw(ind,:)/10),extw(ind,:),...
0146                         10^(dBZ(k,i,j)/10),'linear','extrap')*dz;  
0147         end
0148        elseif k<size(dBZ,1) & k>1
0149 
0150         dz=( [ALT(k,i,j)-ALT(k-1,i,j)] + [ALT(k+1,i,j)-ALT(k,i,j)] )/2;
0151         if ice
0152          qext(k)= interp1(10.^(dBZm(ind,:)/10),ext(ind,:),...
0153                          10^(dBZp(k,i,j)/10),'linear','extrap')*dz;
0154         else
0155          qext(k)= interp1(10.^(dBZmw(ind,:)/10),extw(ind,:),...
0156                          10^(dBZp(k,i,j)/10),'linear','extrap')*dz;
0157         end
0158        else
0159 
0160         dz=ALT(k+1,i,j)-ALT(k,i,j);
0161         if ice
0162          qext(k)=interp1(10.^(dBZm(ind,:)/10),ext(ind,:),...
0163                         10^(dBZp(k,i,j)/10),'linear','extrap')*dz;
0164         else
0165          qext(k)=interp1(10.^(dBZmw(ind,:)/10),extw(ind,:),...
0166                         10^(dBZp(k,i,j)/10),'linear','extrap')*dz;
0167         end 
0168        end
0169        %if qext(k)>10
0170        %   qext(k)=10;
0171        %end
0172        %two way transmission
0173        sqext2(k)=2*sum(qext(k:size(ALT,1)));
0174        trans2(k)=exp(-sqext2(k));
0175        %self layer transmission
0176        trans1(k)=exp(-qext(k));
0177 
0178        if k~=1 & k~=size(ALT,1)
0179         dBZp(k,i,j)=10*log10(10^(dBZp(k,i,j)/10)/trans1(k)); 
0180         dBZp(k-1,i,j)=10*log10(10^(dBZ(k-1,i,j)/10)/trans2(k)); 
0181        elseif k==1
0182         dBZp(k,i,j)=10*log10(10^(dBZp(k,i,j)/10)/trans1(k));
0183        else
0184         dBZp(k,i,j)=10*log10(10^(dBZ(k,i,j)/10)/trans1(k));
0185         dBZp(k-1,i,j)=10*log10(10^(dBZ(k-1,i,j)/10)/trans2(k));
0186        end
0187      end
0188     end
0189    end
0190       
0191    
0192 
0193    %set all modified dBZp to min(dBZ) for these original dBZ==min(dBZ)
0194    %since otherwise we may strecth out the clouds
0195    mindBZ=min(min(min(dBZ)));
0196    dBZp(find(dBZ==mindBZ))=mindBZ;
0197 
0198    end
0199 
0200    %re-arrange the look-up-table [IWC(T,dBZ)]
0201    dBZv=[-49:1:30];
0202    %dBZv=[-49:1:50];
0203    for i=1:size(dBZm,1)
0204        IWCC(i,:)=interp1(dBZm(i,:),IWCv,dBZv);
0205    end
0206    %IWC=interp2(dBZv,Tvi,IWCC,dBZp,T);
0207    IWC=interp2(dBZv,Tvi,IWCC,dBZ,T);
0208    IWC(isnan(IWC))=0;
0209    IWC(find(T>273))=0;
0210    
0211    %re-arrange the look-up-table [LWC(T,dBZ)]
0212    dBZv=[-49:1:30];
0213    %dBZv=[-49:1:50];
0214    for i=1:size(dBZmw,1)
0215        LWCC(i,:)=interp1(dBZmw(i,:),LWCv,dBZv);
0216    end
0217    %LWC=interp2(dBZv,Tvw,LWCC,dBZp,T);
0218    LWC=interp2(dBZv,Tvw,LWCC,dBZ,T);
0219    LWC(find(T<=273))=0;
0220    LWC(isnan(LWC))=0;
0221   
0222  
0223 
0224 
0225    %put in IWC and LWC in G
0226    glen=length(G);   
0227    G(glen+1)=G(glen);
0228    G(glen+1).NAME='IWC field';
0229    G(glen+1).DATA=IWC;
0230    G(glen+1).DATA_NAME='IWC';
0231    G(glen+1).DATA_UNIT='gm-3';
0232    G(glen+1).SOURCE='';
0233    
0234    %setting in PND fields in G
0235    Ni=length(P(1).x);
0236    if ~strcmp(P(1).shape,'sphere')
0237       error('sphere is the only valid option for P.shape')
0238    end
0239    if strcmp(P(1).PSD,'MH97')
0240       y=zeros(length(IWC(:)),Ni);
0241       ivec=find(IWC);
0242       for i=ivec'
0243           y(i,:) = ice_psd_Mcfar_97(T(i),IWC(i),P(1).x,1);
0244       end
0245       Y = gauss_laguerre_apply(y',P(1).x,P(1).w,P(1).x_norm);
0246       for i=1:Ni
0247           G(glen+1+i)=G(glen+1);
0248           G(glen+1+i).DATA=reshape(Y(i,:),size(IWC));
0249           G(glen+1+i).PROPS=[];
0250           G(glen+1+i).NAME=['Particles ',num2str(P(1).x(i)*1e6),' um'];
0251           G(glen+1+i).DATA_NAME=['Particle number density'];
0252           G(glen+1+i).DATA_UNIT=['m-3'];
0253       end
0254    end
0255  
0256    glen=length(G);   
0257    G(glen+1)=G(glen);
0258    G(glen+1).NAME='LWC field';
0259    G(glen+1).DATA=LWC;
0260    G(glen+1).DATA_NAME='IWC';
0261    G(glen+1).DATA_UNIT='gm-3';
0262    G(glen+1).SOURCE='';
0263    
0264    %setting in PND fields in G
0265    Ni=length(P(2).x);
0266    if ~strcmp(P(2).shape,'sphere')
0267       error('sphere is the only valid option for P.shape')
0268    end
0269    if 0 %strcmp(P(4).PSD,'Water')
0270       y=zeros(length(LWC(:)),Ni);
0271       ivec=find(LWC);
0272       for i=ivec'
0273           %N=6*IWC(i).*((3.67+a)./dm).^(4+a)/(rho*pi*gamma(4+a));
0274           %y(i,:)=N*P(4).x.^a.*exp(-(3.67+a)*P(2).x/dm);
0275           %y(i,:) = ice_psd_Mcfar_97(T(i),IWC(i),P(2).x,1);
0276           y(i,:)=water_psd(LWC(i),P(4).x,rc,c1,c2)';
0277       end
0278       Y = gauss_laguerre_apply(y',P(4).x,P(4).w,P(4).x_norm);
0279       for i=1:Ni
0280           G(glen+1+i)=G(glen+1);
0281           G(glen+1+i).DATA=reshape(Y(i,:),size(IWC));
0282           G(glen+1+i).PROPS=P(4).SSP(i);
0283           G(glen+1+i).NAME=['Particles ',num2str(P(4).x(i)*1e6),' um'];
0284           G(glen+1+i).DATA_NAME=['Particle number density'];
0285           G(glen+1+i).DATA_UNIT=['m-3'];
0286       end
0287    end
0288    
0289 
0290 else 
0291   error('gauss-laguerre is the only valid method')
0292 end
0293 
0294 
0295 warning on

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