wrf-python/fortran/wrf_user_dbz.f90

!     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

    !   Force all Q arrays to be 0.0 or greater.
    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

    !   Input pressure is Pa, but we need hPa in calculations

    IF (sn0 .EQ. 0) THEN
        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
    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,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 = 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

                ! 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

    RETURN

END SUBROUTINE CALCDBZ