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.
191 lines
7.0 KiB
191 lines
7.0 KiB
! This routine computes equivalent reflectivity factor (in dBZ) at |
|
! each model grid point. 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.) |
|
! |
|
! 4. If iliqskin=1, frozen particles that are at a temperature above |
|
! freezing are assumed to scatter as a liquid particle. |
|
! |
|
! More information on the derivation of simulated reflectivity in |
|
! RIP can be found in Stoelinga (2005, unpublished write-up). |
|
! Contact Mark Stoelinga (stoeling@atmos.washington.edu) for a copy. |
|
|
|
!NCLFORTSTART |
|
SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx, ny, nz) |
|
USE wrf_constants, ONLY : GAMMA_SEVEN, RHOWAT, RHO_R, RHO_S, RHO_G, ALPHA, & |
|
CELKEL, PI, RD |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: dbz |
|
|
|
! Arguments |
|
INTEGER, INTENT(IN) :: nx, ny, nz |
|
INTEGER, INTENT(IN) :: sn0, ivarint, iliqskin |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: dbz |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: prs |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: tmk |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qvp |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qra |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qsn |
|
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qgr |
|
|
|
!NCLEND |
|
|
|
! Local Variables |
|
INTEGER :: i, j, k |
|
REAL(KIND=8) :: temp_c, virtual_t |
|
REAL(KIND=8) :: gonv, ronv, sonv |
|
REAL(KIND=8) :: factor_g, factor_r, factor_s |
|
REAL(KIND=8) :: factorb_g, factorb_s |
|
REAL(KIND=8) :: rhoair, z_e |
|
|
|
! Constants used to calculate variable intercepts |
|
REAL(KIND=8), PARAMETER :: R1 = 1.D-15 |
|
REAL(KIND=8), PARAMETER :: RON = 8.D6 |
|
REAL(KIND=8), PARAMETER :: RON2 = 1.D10 |
|
REAL(KIND=8), PARAMETER :: SON = 2.D7 |
|
REAL(KIND=8), PARAMETER :: GON = 5.D7 |
|
REAL(KIND=8), PARAMETER :: RON_MIN = 8.D6 |
|
REAL(KIND=8), PARAMETER :: RON_QR0 = 0.00010D0 |
|
REAL(KIND=8), PARAMETER :: RON_DELQR0 = 0.25D0*RON_QR0 |
|
REAL(KIND=8), PARAMETER :: RON_CONST1R = (RON2-RON_MIN)*0.5D0 |
|
REAL(KIND=8), PARAMETER :: RON_CONST2R = (RON2+RON_MIN)*0.5D0 |
|
|
|
! Constant intercepts |
|
REAL(KIND=8), PARAMETER :: RN0_R = 8.D6 |
|
REAL(KIND=8), PARAMETER :: RN0_S = 2.D7 |
|
REAL(KIND=8), PARAMETER :: RN0_G = 4.D6 |
|
|
|
!$OMP PARALLEL |
|
|
|
! Force all Q arrays to be 0.0 or greater. |
|
!$OMP DO COLLAPSE(3) SCHEDULE(runtime) |
|
DO k = 1,nz |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
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 |
|
!$OMP END DO |
|
|
|
! Input pressure is Pa, but we need hPa in calculations |
|
|
|
IF (sn0 .EQ. 0) THEN |
|
!$OMP DO COLLAPSE(3) SCHEDULE(runtime) |
|
DO k = 1,nz |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
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 |
|
!$OMP 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 |
|
|
|
|
|
!$OMP DO COLLAPSE(3) PRIVATE(i, j, k, temp_c, virtual_t, gonv, ronv, sonv, & |
|
!$OMP factorb_g, factorb_s, rhoair, z_e) & |
|
!$OMP FIRSTPRIVATE(factor_r, factor_s, factor_g) SCHEDULE(runtime) |
|
DO k = 1,nz |
|
DO j = 1,ny |
|
DO i = 1,nx |
|
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) |
|
|
|
! 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. |
|
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 |
|
|
|
! Calculate variable intercept parameters |
|
IF (ivarint .EQ. 1) THEN |
|
|
|
temp_c = MIN(-0.001D0, tmk(i,j,k)-CELKEL) |
|
sonv = MIN(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 |
|
|
|
! 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(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 |
|
|
|
! Adjust small values of Z_e so that dBZ is no lower than -30 |
|
z_e = MAX(z_e, .001D0) |
|
|
|
! Convert to dBZ |
|
dbz(i,j,k) = 10.D0*LOG10(z_e) |
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END DO |
|
|
|
!$OMP END PARALLEL |
|
|
|
RETURN |
|
|
|
END SUBROUTINE CALCDBZ |
|
|
|
|