Home > atmlab > arts > scenegen > asg_pathiwc.m

asg_pathiwc

PURPOSE ^

ASG_PATHIWC calulates iwc and iwp and rhi along geometric propagation path

SYNOPSIS ^

function P=asg_pathiwc(G,Q,z_max)

DESCRIPTION ^

ASG_PATHIWC calulates iwc and iwp and rhi along geometric propagation path

     This function calculates the geometric propagation
     path from a given height to the tangent point in steps 
     of 100 m. The IWC field in G is simply linearly
     3-D interpolated (in altitude, latitude,longitude)
     on this path. Reguired data in G is IWC, Altitude
     Temperature, and Water vapour field. 

OUT  P.z   altitude along propagation path
     P.lat latitude along propagation path
     P.lon longitude along propagation path
     P.iwc ice water content along propagation path
     P.iwp vertical ice water path above the tangent point
     P.rhi relative humidity w.r.t. ice
    
IN   
     G     ASG data.
     Q     Qarts setting structure
           Only Q.SENSOR_LOS and Q.SENSOR_POS are considered
     z_max maximum altitude of the propagation path
           taking into consideration

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

asg_pathiwc.m

SOURCE CODE ^

0001 %ASG_PATHIWC calulates iwc and iwp and rhi along geometric propagation path
0002 %
0003 %     This function calculates the geometric propagation
0004 %     path from a given height to the tangent point in steps
0005 %     of 100 m. The IWC field in G is simply linearly
0006 %     3-D interpolated (in altitude, latitude,longitude)
0007 %     on this path. Reguired data in G is IWC, Altitude
0008 %     Temperature, and Water vapour field.
0009 %
0010 %OUT  P.z   altitude along propagation path
0011 %     P.lat latitude along propagation path
0012 %     P.lon longitude along propagation path
0013 %     P.iwc ice water content along propagation path
0014 %     P.iwp vertical ice water path above the tangent point
0015 %     P.rhi relative humidity w.r.t. ice
0016 %
0017 %IN
0018 %     G     ASG data.
0019 %     Q     Qarts setting structure
0020 %           Only Q.SENSOR_LOS and Q.SENSOR_POS are considered
0021 %     z_max maximum altitude of the propagation path
0022 %           taking into consideration
0023 
0024 % 2007-12-06   Created by Bengt Rydberg
0025 
0026 function P=asg_pathiwc(G,Q,z_max)
0027 
0028 
0029 %find the data of interest
0030 iwc_ind=find(strncmp(lower({G.NAME}),'iwc',3));
0031 alt_ind=find(strncmp(lower({G.NAME}),'altitude',8));
0032 hum_ind=min(find(strcmp(lower({G.NAME}),'h2o')));
0033 tem_ind=find(strncmp(lower({G.NAME}),'temperature',11));
0034 lwc_ind=find(strncmp(lower({G.NAME}),'iwc',3));
0035 
0036 
0037 if isempty( iwc_ind )
0038    error('G must contain IWC field')
0039 end
0040 
0041 if isempty( alt_ind )
0042    error('G must contain Altitude field')
0043 end
0044 
0045 if isempty( hum_ind )
0046    error('G must contain Water vapour')
0047 end
0048 
0049 if isempty( tem_ind )
0050    error('G must contain Temperature field')
0051 end
0052 
0053 if isempty( lwc_ind )
0054    error('G must contain LWC field')
0055 end
0056 
0057 deg2rad = constants( 'DEG2RAD' );
0058 
0059 %regrid the fields on the same grids
0060 
0061 Q2 = asg_atmgrids( D, G, Q );
0062 
0063 %remove endpoints since these can cause problems
0064 %in interpolation, we are not anyway interested
0065 %in these
0066     
0067 Q2.P_GRID=Q2.P_GRID(2:end-1);
0068 
0069 if ~isempty(iwc_ind)
0070     %H=asg_regrid(G([alt_ind,hum_ind,tem_ind,iwc_ind,lwc_ind]),Q2);
0071     H=G([alt_ind,hum_ind,tem_ind,iwc_ind,lwc_ind]);
0072 else
0073     %H=asg_regrid(G([alt_ind,hum_ind,tem_ind]),Q2);
0074     G([alt_ind,hum_ind,tem_ind]);
0075 end
0076 
0077 %create a pressure matrix
0078 
0079 Pr=zeros(size(H(1).DATA));
0080 for i=1:length(H(1).GRID3)
0081     Pr(:,:,i)=vec2col(H(1).GRID1)*ones(1,length(H(1).GRID2));
0082 end
0083 
0084 %create meshgrids to be used in interpolation
0085 [X,Y,Z]=meshgrid(H(1).GRID2,H(1).DATA(:,round(end/2),round(end/2)),H(1).GRID3);
0086 
0087 %loop over sensor positions and line of sights
0088 %
0089 
0090 for pp=1:size(Q.SENSOR_POS,1)
0091 
0092     pos=Q.SENSOR_POS(pp,:);
0093     los=Q.SENSOR_LOS(pp,:);
0094     
0095     [x,y,z,dx,dy,dz] = arts_poslos2cart( pos(1),pos(2),pos(3),los(1),los(2) );
0096 
0097     %find the propagation path from the sensor to the tangent point
0098     %
0099 
0100     %the length from the sensor to the tangent point
0101     l_sen=pos(1)*cos( (180-los(1)) *deg2rad);
0102 
0103     %length of each step
0104     l_step=1e2;
0105 
0106     %make a vector of propagation steps
0107     l_sen=unique([ [0:l_step:l_sen] l_sen]);
0108 
0109     %the propagation path in cartesian cord.
0110     [x,y,z]=deal(x+dx*l_sen,y+dy*l_sen,z+dz*l_sen);
0111 
0112     %the propagation path in spherical cord.
0113     [r,lat,lon] = arts_cart2sph(x,y,z);
0114 
0115     %the geoid radius at the tangent point
0116     %re = interp1( [-90:1:90], wgs84( 2, [-90:1:90]), lat(end) );
0117     re = Q.R_GEOID;  
0118 
0119     %interpolate the data field on the propagation path
0120     %
0121 
0122     %remove high altitude points
0123     z_grid=r-re;
0124     ind=find(z_grid<z_max & z_grid>-1e3);
0125     
0126     z_grid=z_grid(ind);
0127     lat=lat(ind);
0128     lon=lon(ind);
0129    
0130     if ~isempty(iwc_ind)
0131        IWC=interp3(X,Y,Z,H(4).DATA,lat,z_grid,lon);
0132     end
0133     HUM=interp3(X,Y,Z,H(2).DATA,lat,z_grid,lon);
0134     TEM=interp3(X,Y,Z,H(3).DATA,lat,z_grid,lon);
0135     PRE=exp(interp3(X,Y,Z,log(Pr),lat,z_grid,lon));
0136     LWC=interp3(X,Y,Z,H(5).DATA,lat,z_grid,lon);
0137 
0138     %Equilibrium water vapor pressure
0139     ei = e_eq_ice(TEM);
0140 
0141     %Relative humidity w.r.t. ice
0142     RH=PRE.*HUM./ei*100;
0143   
0144     P(pp).z=z_grid;
0145     P(pp).lat=lat;
0146     P(pp).lon=lon;
0147     if ~isempty(iwc_ind)
0148        %set all nan values to 0
0149        IWC(find(isnan(IWC)))=0;
0150        LWC(find(isnan(LWC)))=0;
0151        P(pp).iwc=IWC;
0152        P(pp).lwc=LWC;
0153        P(pp).iwp=sum( (P(pp).z(1:end-1)-P(pp).z(2:end))...
0154                     .* (IWC(1:end-1)+IWC(2:end))/2 );
0155     else
0156       P(pp).iwc=zeros(size(z_grid));
0157       P(pp).iwp=0;
0158     end
0159     P(pp).rhi=RH;
0160     P(pp).t=TEM;
0161 
0162 end

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