forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
209 lines
7.1 KiB
209 lines
7.1 KiB
C NCLFORTSTART |
|
SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, |
|
+ SN0,IVARINT,ILIQSKIN) |
|
c |
|
c This routine computes equivalent reflectivity factor (in dBZ) at |
|
c each model grid point. In calculating Ze, the RIP algorithm makes |
|
c assumptions consistent with those made in an early version |
|
c (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5 |
|
c model (i.e., the scheme known as "Resiner-2"). For each species: |
|
c |
|
c 1. Particles are assumed to be spheres of constant density. The |
|
c densities of rain drops, snow particles, and graupel particles are |
|
c taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and |
|
c rho_g = 400 kg m^-3, respectively. (l refers to the density of |
|
c liquid water.) |
|
c |
|
c 2. The size distribution (in terms of the actual diameter of the |
|
c particles, rather than the melted diameter or the equivalent solid |
|
c ice sphere diameter) is assumed to follow an exponential |
|
c distribution of the form N(D) = N_0 * exp( lambda*D ). |
|
c |
|
c 3. If ivarint=0, the intercept parameters are assumed constant |
|
c (as in early Reisner-2), with values of 8x10^6, 2x10^7, |
|
c and 4x10^6 m^-4, for rain, snow, and graupel, respectively. |
|
c If ivarint=1, variable intercept parameters are used, as |
|
c calculated in Thompson, Rasmussen, and Manning (2004, Monthly |
|
c Weather Review, Vol. 132, No. 2, pp. 519-542.) |
|
c |
|
c 4. If iliqskin=1, frozen particles that are at a temperature above |
|
c freezing are assumed to scatter as a liquid particle. |
|
c |
|
c More information on the derivation of simulated reflectivity in |
|
c RIP can be found in Stoelinga (2005, unpublished write-up). |
|
c Contact Mark Stoelinga (stoeling@atmos.washington.edu) for a copy. |
|
c |
|
|
|
c Arguments |
|
INTEGER WEDIM,SNDIM,BTDIM |
|
INTEGER SN0,IVARINT,ILIQSKIN |
|
DOUBLE PRECISION DBZ(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION PRS(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION TMK(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION QVP(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION QRA(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION QSN(WEDIM,SNDIM,BTDIM) |
|
DOUBLE PRECISION QGR(WEDIM,SNDIM,BTDIM) |
|
|
|
C NCLEND |
|
|
|
c Local Variables |
|
INTEGER I,J,K |
|
DOUBLE PRECISION TEMP_C,VIRTUAL_T |
|
DOUBLE PRECISION GONV,RONV,SONV |
|
DOUBLE PRECISION FACTOR_G,FACTOR_R,FACTOR_S |
|
DOUBLE PRECISION FACTORB_G,FACTORB_R,FACTORB_S |
|
DOUBLE PRECISION RHOAIR,Z_E |
|
|
|
c Constants used to calculate variable intercepts |
|
DOUBLE PRECISION R1,RON,RON2,SON,GON |
|
DOUBLE PRECISION RON_MIN,RON_QR0,RON_DELQR0 |
|
DOUBLE PRECISION RON_CONST1R,RON_CONST2R |
|
c Constant intercepts |
|
DOUBLE PRECISION RN0_R,RN0_S,RN0_G |
|
c Other constants |
|
DOUBLE PRECISION RHO_R,RHO_S,RHO_G |
|
DOUBLE PRECISION GAMMA_SEVEN,ALPHA |
|
DOUBLE PRECISION RHOWAT,CELKEL,PI,RD |
|
|
|
|
|
c Constants used to calculate variable intercepts |
|
R1 = 1.D-15 |
|
RON = 8.D6 |
|
RON2 = 1.D10 |
|
SON = 2.D7 |
|
GON = 5.D7 |
|
RON_MIN = 8.D6 |
|
RON_QR0 = 0.00010D0 |
|
RON_DELQR0 = 0.25D0*RON_QR0 |
|
RON_CONST1R = (RON2-RON_MIN)*0.5D0 |
|
RON_CONST2R = (RON2+RON_MIN)*0.5D0 |
|
|
|
c Constant intercepts |
|
RN0_R = 8.D6 |
|
RN0_S = 2.D7 |
|
RN0_G = 4.D6 |
|
|
|
c Other constants |
|
GAMMA_SEVEN = 720.D0 |
|
RHOWAT = 1000.D0 |
|
RHO_R = RHOWAT |
|
RHO_S = 100.D0 |
|
RHO_G = 400.D0 |
|
ALPHA = 0.224D0 |
|
CELKEL = 273.15D0 |
|
PI = 3.141592653589793D0 |
|
RD = 287.04D0 |
|
|
|
c Force all Q arrays to be 0.0 or greater. |
|
DO K = 1,BTDIM |
|
DO J = 1,SNDIM |
|
DO I = 1,WEDIM |
|
IF (QVP(I,J,K).LT.0.0) THEN |
|
QVP(I,J,K) = 0.0 |
|
END IF |
|
IF (QRA(I,J,K).LT.0.0) THEN |
|
QRA(I,J,K) = 0.0 |
|
END IF |
|
IF (QSN(I,J,K).LT.0.0) THEN |
|
QSN(I,J,K) = 0.0 |
|
END IF |
|
IF (QGR(I,J,K).LT.0.0) THEN |
|
QGR(I,J,K) = 0.0 |
|
END IF |
|
END DO |
|
END DO |
|
END DO |
|
|
|
c Input pressure is Pa, but we need hPa in calculations |
|
|
|
IF (SN0.EQ.0) THEN |
|
DO K = 1,BTDIM |
|
DO J = 1,SNDIM |
|
DO I = 1,WEDIM |
|
IF (TMK(I,J,K).LT.CELKEL) THEN |
|
QSN(I,J,K) = QRA(I,J,K) |
|
QRA(I,J,K) = 0.D0 |
|
END IF |
|
END DO |
|
END DO |
|
END DO |
|
END IF |
|
|
|
|
|
FACTOR_R = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_R))**1.75D0 |
|
FACTOR_S = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_S))**1.75D0* |
|
+ (RHO_S/RHOWAT)**2*ALPHA |
|
FACTOR_G = GAMMA_SEVEN*1.D18* (1.D0/ (PI*RHO_G))**1.75D0* |
|
+ (RHO_G/RHOWAT)**2*ALPHA |
|
|
|
|
|
DO K = 1,BTDIM |
|
DO J = 1,SNDIM |
|
DO I = 1,WEDIM |
|
|
|
VIRTUAL_T = TMK(I,J,K)* (0.622D0+QVP(I,J,K))/ |
|
+ (0.622D0* (1.D0+QVP(I,J,K))) |
|
RHOAIR = PRS(I,J,K) / (RD*VIRTUAL_T) |
|
|
|
c Adjust factor for brightband, where snow or graupel particle |
|
c scatters like liquid water (alpha=1.0) because it is assumed to |
|
c have a liquid skin. |
|
|
|
IF (ILIQSKIN.EQ.1 .AND. TMK(I,J,K).GT.CELKEL) THEN |
|
FACTORB_S = FACTOR_S/ALPHA |
|
FACTORB_G = FACTOR_G/ALPHA |
|
ELSE |
|
FACTORB_S = FACTOR_S |
|
FACTORB_G = FACTOR_G |
|
END IF |
|
|
|
c Calculate variable intercept parameters |
|
|
|
IF (IVARINT.EQ.1) THEN |
|
|
|
TEMP_C = DMIN1(-0.001D0,TMK(I,J,K)-CELKEL) |
|
SONV = DMIN1(2.0D8,2.0D6*EXP(-0.12D0*TEMP_C)) |
|
|
|
GONV = GON |
|
IF (QGR(I,J,K).GT.R1) THEN |
|
GONV = 2.38D0* (PI*RHO_G/ |
|
+ (RHOAIR*QGR(I,J,K)))**0.92D0 |
|
GONV = MAX(1.D4,MIN(GONV,GON)) |
|
END IF |
|
|
|
RONV = RON2 |
|
IF (QRA(I,J,K).GT.R1) THEN |
|
RONV = RON_CONST1R*TANH((RON_QR0-QRA(I,J,K))/ |
|
+ RON_DELQR0) + RON_CONST2R |
|
END IF |
|
|
|
ELSE |
|
|
|
RONV = RN0_R |
|
SONV = RN0_S |
|
GONV = RN0_G |
|
|
|
END IF |
|
|
|
c Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is |
|
c the sum of z_e for each hydrometeor species: |
|
|
|
Z_E = FACTOR_R* (RHOAIR*QRA(I,J,K))**1.75D0/ |
|
+ RONV**.75D0 + FACTORB_S* |
|
+ (RHOAIR*QSN(I,J,K))**1.75D0/SONV**.75D0 + |
|
+ FACTORB_G* (RHOAIR*QGR(I,J,K))**1.75D0/ |
|
+ GONV**.75D0 |
|
|
|
c Adjust small values of Z_e so that dBZ is no lower than -30 |
|
Z_E = MAX(Z_E,.001D0) |
|
|
|
c Convert to dBZ |
|
DBZ(I,J,K) = 10.D0*LOG10(Z_E) |
|
|
|
END DO |
|
END DO |
|
END DO |
|
|
|
RETURN |
|
END
|
|
|