wrf_post_wrf后处理包_matlab_wrf_WRF后处理_wrf_post.m_

function [ dbz,maxdbz ] = calc_dbz( prs,tmk,QV,QR,QS,QG ) %CALC_DBZ To compute equivalent reflectivity factor Ze(in dBZ) at each % model grid point(nx*ny*nz*nt), and maximum reflectivity Ze-max of all % levels at each horizontal model grid point(nx*ny). Where nx,ny are % number of horizontal grid points and nz is vertical. % %OUTPUT % dbz: equivalent reflectivity factor Ze(unit dBZ, size nx*ny*nz); % maxdbz: maximum reflectivity Ze-max(unit dBZ, size nx*ny); %INPUT % prs: pressure on each level (unit hPa); % tmk: temperature (unit K); % QV: Water vapor mixing ratio (unit kg/kg); % QR: Rain water mixing ratio (unit kg/kg); % QS: Snow mixing ratio (unit kg/kg), missing to set as -1; % QG: Graupel mixing ratio (unit kg/kg), missing to set as -1. % %USAGE % DBZ=clac_dbz(prs,tmk,QV,QR,QS,QG) calculate Ze in dBZ with QS and QG % data available. % [DBZ,MAXDBZ]=calc_dbz(prs,tmk,QV,QR,QS,QG) calculate Ze and Ze-max in % dBZ with QS and QG data available. % DBZ=clac_dbz(prs,tmk,QV,QR,-1,-1) or % [DBZ,MAXDBZ]=calc_dbz(prs,tmk,QV,QR,-1,QG) when QS or QG data missing. % %ASSUMPTIONS % In calculating Ze, the RIP algorithm makes assumptions consistent with % those made in an early version (ca. 1996) of the bulk mixed-phase % microphysical scheme in the MM5 model (i.e., the scheme known as % "Resiner-2"). For each species: % % 1. Particles are assumed to be spheres of constant density. The % densities of rain drops, snow particles, and graupel particles are % taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and % rho_g = 400 kg m^-3, respectively. (l refers to the density of % liquid water.) % % 2. The size distribution (in terms of the actual diameter of the % particles, rather than the melted diameter or the equivalent solid % ice sphere diameter) is assumed to follow an exponential % distribution of the form N(D) = N_0 * exp( lambda*D ). % % 3. If ivarint=0, the intercept parameters are assumed constant (as in % early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4, for % rain, snow, and graupel, respectively. If ivarint=1, variable intercept % parameters are used, as calculated in Thompson, Rasmussen, and Manning % (2004, Monthly Weather Review, Vol. 132, No. 2, pp. 519-542.) % %STATEMENT % This program is adapted from Fortran program module_calc_dbz.f90. % % AUTHOR: sfhstcn2 % CONTACT: sfh_st_cn2@163.com % Checking size of arrays if any([ndims(prs) ndims(tmk) ndims(QV) ndims(QR) ndims(QS) ndims(QG)]>4) error('Input arrays has exceeded dimensions.') elseif ~isequal(size(prs),size(tmk),size(QV),size(QR)) error('Inconsistent dimensions for tmk, QV and QR.') end if QS~=-1 if ~isequal(size(tmk),size(QV),size(QR),size(QS)) error('Inconsistent dimensions for QS and others.') end end if QG~=-1 if ~isequal(size(tmk),size(QV),size(QR),size(QG)) error('Inconsistent dimensions for QG and others.') end end [nx,ny,nz,nt]=size(tmk); % Parameters CELKEL = 273.15; gamma7 = 720.; RHOWAT = 1000.; % kg*m^3 rho_r = RHOWAT; rho_s = 100.; % kg*m^3 rho_g = 400.; % kg*m^3 alpha = .224; Rd = 287.04; gon = 5.e7; r1 = 1.e-15; ron2 = 1.e10; ron_min = 8.e6; ron_qr0 = .00010; ron_delqr0 = .25*ron_qr0; ron_const1r = (ron2-ron_min)*.5; ron_const2r = (ron2+ron_min)*.5; qvp = QV; qra = QR; qsn = zeros(nx,ny,nz,nt); qgr = zeros(nx,ny,nz,nt); if QS~=-1 qsn = QS; end if QG~=-1 qgr = QG; end if max(max(max(max(qsn))))==0.0 && max(max(max(max(qgr))))==0.0 % No ice in this run qra(tmk<CELKEL) = 0.; qsn(tmk<CELKEL) = QR(tmk<CELKEL); end qvp(~isnan(qvp)) = max(qvp(~isnan(qvp)),0.); qra(~isnan(qra)) = max(qra(~isnan(qra)),0.); qsn(~isnan(qsn)) = max(qsn(~isnan(qsn)),0.); qgr(~isnan(qgr)) = max(qgr(~isnan(qgr)),0.); factor_r = gamma7*1.e18*(1/pi/rho_r)^1.75; factor_s = gamma7*1.e18*(1/pi/rho_s)^1.75*(rho_s/RHOWAT)^2*alpha; factor_g = gamma7*1.e18*(1/pi/rho_g)^1.75*(rho_g/RHOWAT)^2*alpha; virtual = tmk.*(.622+qvp)./(.622*(1.+qvp)); rhoair = prs*100./(Rd*virtual); % Adjust factor for brightband, where snow or graupel particle scatters % like liquid water (alpha=1.0) because it is assumed to have a liquid % skin. factorb_s(1:nx,1:ny,1:nz,1:nt) = factor_s; factorb_g(1:nx,1:ny,1:nz,1:nt) = factor_g; factorb_s(tmk>CELKEL) = factor_s/double(alpha); factorb_g(tmk>CELKEL) = factor_g/double(alpha); % Calculate variable intercept parameters temp_c = min(-.001,tmk-CELKEL); sonv = min(2.e8,2.e6*exp(-.12*temp_c)); gonv(1:nx,1:ny,1:nz,1:nt) = gon; gonv(qgr>r1) = 2.38*(pi*rho_g./(rhoair(qgr>r1).*qgr(qgr>r1))).^.92; gonv(qgr>r1) = max(1.e4,min(gonv(qgr>r1),gon)); ronv(1:nx,1:ny,1:nz,1:nt) = ron2; ronv(qra>r1) = ron_const1r*tanh((ron_qr0-qra(qra>r1))/ron_delqr0)+ron_const2r; % Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is the sum of % z_e for each hydrometeor species: z_e = factor_r.*(rhoair.*qra).^1.75 ./ronv.^.75 + ... factorb_s.*(rhoair.*qsn).^1.75 ./sonv.^.75 + ... factorb_g.*(rhoair.*qgr).^1.75 ./gonv.^.75; % Adjust small values of Z_e so that dBZ is no lower than -30 z_e = max(z_e,.001); % Convert to dBZ dbz = 10.*log10(z_e); % Calculate the maximum reflectivity maxdbz(:,:,:) = max(dbz,[],3); end