| 1 | !! Diagnostics: Reflectivity |
|---|
| 2 | !! This routine - including comments were taken from the RIP4 program |
|---|
| 3 | !! The only change was the removing of the old "ivarint=0" option |
|---|
| 4 | |
|---|
| 5 | MODULE module_calc_dbz |
|---|
| 6 | |
|---|
| 7 | CONTAINS |
|---|
| 8 | SUBROUTINE calc_dbz (SCR, cname, cdesc, cunits) |
|---|
| 9 | |
|---|
| 10 | ! |
|---|
| 11 | ! This routine computes equivalent reflectivity factor (in dBZ) at |
|---|
| 12 | ! each model grid point. In calculating Ze, the RIP algorithm makes |
|---|
| 13 | ! assumptions consistent with those made in an early version |
|---|
| 14 | ! (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5 |
|---|
| 15 | ! model (i.e., the scheme known as "Resiner-2"). For each species: |
|---|
| 16 | ! |
|---|
| 17 | ! 1. Particles are assumed to be spheres of constant density. The |
|---|
| 18 | ! densities of rain drops, snow particles, and graupel particles are |
|---|
| 19 | ! taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and |
|---|
| 20 | ! rho_g = 400 kg m^-3, respectively. (l refers to the density of |
|---|
| 21 | ! liquid water.) |
|---|
| 22 | ! |
|---|
| 23 | ! 2. The size distribution (in terms of the actual diameter of the |
|---|
| 24 | ! particles, rather than the melted diameter or the equivalent solid |
|---|
| 25 | ! ice sphere diameter) is assumed to follow an exponential |
|---|
| 26 | ! distribution of the form N(D) = N_0 * exp( lambda*D ). |
|---|
| 27 | ! |
|---|
| 28 | ! 3. If ivarint=0, the intercept parameters are assumed constant (as |
|---|
| 29 | ! in early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4, |
|---|
| 30 | ! for rain, snow, and graupel, respectively. If ivarint=1, variable |
|---|
| 31 | ! intercept parameters are used, as calculated in Thompson, Rasmussen, |
|---|
| 32 | ! and Manning (2004, Monthly Weather Review, Vol. 132, No. 2, pp. 519-542.) |
|---|
| 33 | ! |
|---|
| 34 | ! More information on the derivation of simulated reflectivity in RIP |
|---|
| 35 | ! can be found in Stoelinga (2005, unpublished write-up). Contact |
|---|
| 36 | ! Mark Stoelinga (stoeling@atmos.washington.edu) for a copy. |
|---|
| 37 | ! |
|---|
| 38 | |
|---|
| 39 | USE constants_module |
|---|
| 40 | USE module_model_basics |
|---|
| 41 | |
|---|
| 42 | IMPLICIT NONE |
|---|
| 43 | |
|---|
| 44 | !Arguments |
|---|
| 45 | real, pointer, dimension(:,:,:) :: SCR |
|---|
| 46 | character (len=128) :: cname, cdesc, cunits |
|---|
| 47 | |
|---|
| 48 | !Local Variables |
|---|
| 49 | integer :: i, j, k |
|---|
| 50 | real :: temp_c |
|---|
| 51 | real :: gonv, ronv, sonv |
|---|
| 52 | real :: factor_g, factor_r, factor_s |
|---|
| 53 | real :: factorb_g, factorb_r, factorb_s |
|---|
| 54 | real :: rhoair, z_e |
|---|
| 55 | real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: qvp, qra, qsn, qgr |
|---|
| 56 | real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: tmk, prs |
|---|
| 57 | |
|---|
| 58 | ! Constants used to calculate variable intercepts |
|---|
| 59 | real, parameter :: r1 = 1.e-15 |
|---|
| 60 | real, parameter :: ron = 8.e6 |
|---|
| 61 | real, parameter :: ron2 = 1.e10 |
|---|
| 62 | real, parameter :: son = 2.e7 |
|---|
| 63 | real, parameter :: gon = 5.e7 |
|---|
| 64 | real, parameter :: ron_min = 8.e6 |
|---|
| 65 | real, parameter :: ron_qr0 = 0.00010 |
|---|
| 66 | real, parameter :: ron_delqr0 = 0.25*ron_qr0 |
|---|
| 67 | real, parameter :: ron_const1r = (ron2-ron_min)*0.5 |
|---|
| 68 | real, parameter :: ron_const2r = (ron2+ron_min)*0.5 |
|---|
| 69 | ! Other constants |
|---|
| 70 | real, parameter :: gamma_seven = 720. |
|---|
| 71 | real, parameter :: rho_r = RHOWAT ! 1000. kg m^-3 |
|---|
| 72 | real, parameter :: rho_s = 100. ! kg m^-3 |
|---|
| 73 | real, parameter :: rho_g = 400. ! kg m^-3 |
|---|
| 74 | real, parameter :: alpha = 0.224 |
|---|
| 75 | |
|---|
| 76 | ALLOCATE(SCR(west_east_dim,south_north_dim,bottom_top_dim)) |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | prs = PRES * 0.01 ! pressure in hPa |
|---|
| 80 | tmk = TK ! temperature in K |
|---|
| 81 | |
|---|
| 82 | qvp = QV |
|---|
| 83 | qra = QR |
|---|
| 84 | qsn = 0.0 |
|---|
| 85 | qgr = 0.0 |
|---|
| 86 | IF ( have_QS ) qsn = QS |
|---|
| 87 | IF ( have_QG ) qgr = QG |
|---|
| 88 | IF ( maxval(qsn) == 0.0 .AND. maxval(qgr) == 0.0 ) THEN |
|---|
| 89 | !! No ice in this run |
|---|
| 90 | WHERE ( tmk .lt. CELKEL ) |
|---|
| 91 | qra = 0.0 |
|---|
| 92 | qsn = qr |
|---|
| 93 | END WHERE |
|---|
| 94 | END IF |
|---|
| 95 | |
|---|
| 96 | factor_r = gamma_seven * 1.e18 * (1./(PI*rho_r))**1.75 |
|---|
| 97 | factor_s = gamma_seven * 1.e18 * (1./(PI*rho_s))**1.75 & |
|---|
| 98 | * (rho_s/RHOWAT)**2 * alpha |
|---|
| 99 | factor_g = gamma_seven * 1.e18 * (1./(PI*rho_g))**1.75 & |
|---|
| 100 | * (rho_g/RHOWAT)**2 * alpha |
|---|
| 101 | |
|---|
| 102 | DO k = 1,bottom_top_dim |
|---|
| 103 | DO j = 1,south_north_dim |
|---|
| 104 | DO i = 1,west_east_dim |
|---|
| 105 | |
|---|
| 106 | rhoair=prs(i,j,k)*100./ ( Rd*virtual(tmk(i,j,k),qvp(i,j,k)) ) ! air density |
|---|
| 107 | |
|---|
| 108 | ! Adjust factor for brightband, where snow or graupel particle |
|---|
| 109 | ! scatters like liquid water (alpha=1.0) because it is assumed to |
|---|
| 110 | ! have a liquid skin. |
|---|
| 111 | |
|---|
| 112 | IF (tmk(i,j,k) .gt. CELKEL) THEN |
|---|
| 113 | factorb_s=factor_s/alpha |
|---|
| 114 | factorb_g=factor_g/alpha |
|---|
| 115 | ELSE |
|---|
| 116 | factorb_s=factor_s |
|---|
| 117 | factorb_g=factor_g |
|---|
| 118 | ENDIF |
|---|
| 119 | |
|---|
| 120 | ! Calculate variable intercept parameters |
|---|
| 121 | |
|---|
| 122 | temp_c = amin1(-0.001, tmk(i,j,k)-CELKEL) |
|---|
| 123 | sonv = amin1(2.0e8, 2.0e6*exp(-0.12*temp_c)) |
|---|
| 124 | |
|---|
| 125 | gonv = gon |
|---|
| 126 | IF (qgr(i,j,k).gt.r1) THEN |
|---|
| 127 | gonv = 2.38*(PI*rho_g/(rhoair*qgr(i,j,k)))**0.92 |
|---|
| 128 | gonv = max(1.e4, min(gonv,gon)) |
|---|
| 129 | ENDIF |
|---|
| 130 | |
|---|
| 131 | ronv = ron2 |
|---|
| 132 | IF (qra(i,j,k).gt. r1) THEN |
|---|
| 133 | ronv = ron_const1r*tanh((ron_qr0-qra(i,j,k)) & |
|---|
| 134 | /ron_delqr0) + ron_const2r |
|---|
| 135 | ENDIF |
|---|
| 136 | |
|---|
| 137 | ! Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is |
|---|
| 138 | ! the sum of z_e for each hydrometeor species: |
|---|
| 139 | |
|---|
| 140 | z_e = factor_r * (rhoair*qra(i,j,k))**1.75 /ronv**.75 & |
|---|
| 141 | + factorb_s * (rhoair*qsn(i,j,k))**1.75 /sonv**.75 & |
|---|
| 142 | + factorb_g * (rhoair*qgr(i,j,k))**1.75 /gonv**.75 |
|---|
| 143 | |
|---|
| 144 | ! Adjust small values of Z_e so that dBZ is no lower than -30 |
|---|
| 145 | z_e = max(z_e,.001) |
|---|
| 146 | |
|---|
| 147 | ! Convert to dBZ |
|---|
| 148 | SCR(i,j,k) = 10. * log10(z_e) |
|---|
| 149 | |
|---|
| 150 | ENDDO |
|---|
| 151 | ENDDO |
|---|
| 152 | ENDDO |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | cname = "dbz" |
|---|
| 156 | cdesc = "Reflectivity" |
|---|
| 157 | cunits = "-" |
|---|
| 158 | |
|---|
| 159 | END SUBROUTINE calc_dbz |
|---|
| 160 | |
|---|
| 161 | END MODULE module_calc_dbz |
|---|
| 162 | |
|---|
