<p><b>laura@ucar.edu</b> 2012-10-22 16:22:09 -0600 (Mon, 22 Oct 2012)</p><p>updated the wsm6 cloud microphysics scheme to include the calculation of radar reflectivity<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_mp_wsm6.F
===================================================================
--- branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_mp_wsm6.F 2012-10-22 22:05:17 UTC (rev 2252)
+++ branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_mp_wsm6.F 2012-10-22 22:22:09 UTC (rev 2253)
@@ -8,6 +8,13 @@
MODULE module_mp_wsm6
!
+!#if defined(non_hydrostatic_core) || defined(hydrostatic_code)
+! USE mpas_atmphys_utilities
+!#else
+! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
+! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
+!#endif
+ USE module_mp_radar
!
REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
@@ -64,6 +71,7 @@
,rain, rainncv &
,snow, snowncv &
,sr &
+ ,refl_10cm, diagflag, do_radar_ref &
,graupel, graupelncv &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
@@ -112,6 +120,16 @@
INTENT(INOUT) :: rain, &
rainncv, &
sr
+
+!+---+-----------------------------------------------------------------+
+#if defined(non_hydrostatic_core) | defined(hydrostatic_core)
+ REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT),optional:: & ! GT
+ refl_10cm
+#else
+ REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
+#endif
+!+---+-----------------------------------------------------------------+
+
REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
INTENT(INOUT) :: snow, &
snowncv
@@ -123,7 +141,13 @@
REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
INTEGER :: i,j,k
-!-------------------------------------------------------------------
+
+!+---+-----------------------------------------------------------------+
+ REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
+ LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
+ INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
+!+---+-----------------------------------------------------------------+
+
DO j=jts,jte
DO k=kts,kte
DO i=its,ite
@@ -163,6 +187,29 @@
qg(i,k,j) = qrs(i,k,3)
ENDDO
ENDDO
+
+!+---+-----------------------------------------------------------------+
+ IF ( PRESENT (diagflag) ) THEN
+ if (diagflag .and. do_radar_ref == 1) then
+ DO I=its,ite
+ DO K=kts,kte
+ t1d(k)=th(i,k,j)*pii(i,k,j)
+ p1d(k)=p(i,k,j)
+ qv1d(k)=q(i,k,j)
+ qr1d(k)=qr(i,k,j)
+ qs1d(k)=qs(i,k,j)
+ qg1d(k)=qg(i,k,j)
+ ENDDO
+ call refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, &
+ t1d, p1d, dBZ, kts, kte, i, j)
+ do k = kts, kte
+ refl_10cm(i,k,j) = MAX(-35., dBZ(k))
+ enddo
+ ENDDO
+ endif
+ ENDIF
+!+---+-----------------------------------------------------------------+
+
ENDDO
END SUBROUTINE wsm6
!===================================================================
@@ -1489,6 +1536,24 @@
rsloper3max = rsloper2max * rslopermax
rslopes3max = rslopes2max * rslopesmax
rslopeg3max = rslopeg2max * rslopegmax
+
+!+---+-----------------------------------------------------------------+
+!..Set these variables needed for computing radar reflectivity. These
+!.. get used within radar_init to create other variables used in the
+!.. radar module.
+ xam_r = PI*denr/6.
+ xbm_r = 3.
+ xmu_r = 0.
+ xam_s = PI*dens/6.
+ xbm_s = 3.
+ xmu_s = 0.
+ xam_g = PI*deng/6.
+ xbm_g = 3.
+ xmu_g = 0.
+
+ call radar_init
+!+---+-----------------------------------------------------------------+
+
!
END SUBROUTINE wsm6init
!------------------------------------------------------------------------------
@@ -2215,4 +2280,182 @@
enddo i_loop
!
END SUBROUTINE nislfv_rain_plm6
+
+!+---+-----------------------------------------------------------------+
+
+ subroutine refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, &
+ t1d, p1d, dBZ, kts, kte, ii, jj)
+
+ IMPLICIT NONE
+
+!..Sub arguments
+ INTEGER, INTENT(IN):: kts, kte, ii, jj
+ REAL, DIMENSION(kts:kte), INTENT(IN):: &
+ qv1d, qr1d, qs1d, qg1d, t1d, p1d
+ REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
+
+!..Local variables
+ REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
+ REAL, DIMENSION(kts:kte):: rr, rs, rg
+ REAL:: temp_C
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
+ DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
+ DOUBLE PRECISION:: lamr, lams, lamg
+ LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
+
+ REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
+ DOUBLE PRECISION:: fmelt_s, fmelt_g
+
+ INTEGER:: i, k, k_0, kbot, n
+ LOGICAL:: melti
+
+ DOUBLE PRECISION:: cback, x, eta, f_d
+ REAL, PARAMETER:: R=287.
+
+!+---+
+
+ do k = kts, kte
+ dBZ(k) = -35.0
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Put column of data into local arrays.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ temp(k) = t1d(k)
+ temp_C = min(-0.001, temp(K)-273.15)
+ qv(k) = MAX(1.E-10, qv1d(k))
+ pres(k) = p1d(k)
+ rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
+
+ if (qr1d(k) .gt. 1.E-9) then
+ rr(k) = qr1d(k)*rho(k)
+ N0_r(k) = n0r
+ lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
+ ilamr(k) = 1./lamr
+ L_qr(k) = .true.
+ else
+ rr(k) = 1.E-12
+ L_qr(k) = .false.
+ endif
+
+ if (qs1d(k) .gt. 1.E-9) then
+ rs(k) = qs1d(k)*rho(k)
+ N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
+ lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
+ ilams(k) = 1./lams
+ L_qs(k) = .true.
+ else
+ rs(k) = 1.E-12
+ L_qs(k) = .false.
+ endif
+
+ if (qg1d(k) .gt. 1.E-9) then
+ rg(k) = qg1d(k)*rho(k)
+ N0_g(k) = n0g
+ lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
+ ilamg(k) = 1./lamg
+ L_qg(k) = .true.
+ else
+ rg(k) = 1.E-12
+ L_qg(k) = .false.
+ endif
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Locate K-level of start of melting (k_0 is level above).
+!+---+-----------------------------------------------------------------+
+ melti = .false.
+ k_0 = kts
+ do k = kte-1, kts, -1
+ if ( (temp(k).gt.273.15) .and. L_qr(k) &
+ .and. (L_qs(k+1).or.L_qg(k+1)) ) then
+ k_0 = MAX(k+1, k_0)
+ melti=.true.
+ goto 195
+ endif
+ enddo
+ 195 continue
+
+!+---+-----------------------------------------------------------------+
+!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
+!.. and non-water-coated snow and graupel when below freezing are
+!.. simple. Integrations of m(D)*m(D)*N(D)*dD.
+!+---+-----------------------------------------------------------------+
+
+ do k = kts, kte
+ ze_rain(k) = 1.e-22
+ ze_snow(k) = 1.e-22
+ ze_graupel(k) = 1.e-22
+ if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
+ if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
+ * (xam_s/900.0)*(xam_s/900.0) &
+ * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
+ if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
+ * (xam_g/900.0)*(xam_g/900.0) &
+ * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
+ enddo
+
+
+!+---+-----------------------------------------------------------------+
+!..Special case of melting ice (snow/graupel) particles. Assume the
+!.. ice is surrounded by the liquid water. Fraction of meltwater is
+!.. extremely simple based on amount found above the melting level.
+!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
+!.. routines).
+!+---+-----------------------------------------------------------------+
+
+ if (melti .and. k_0.ge.kts+1) then
+ do k = k_0-1, kts, -1
+
+!..Reflectivity contributed by melting snow
+ if (L_qs(k) .and. L_qs(k_0) ) then
+ fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
+ eta = 0.d0
+ lams = 1./ilams(k)
+ do n = 1, nrbins
+ x = xam_s * xxDs(n)**xbm_s
+ call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
+ fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
+ CBACK, mixingrulestring_s, matrixstring_s, &
+ inclusionstring_s, hoststring_s, &
+ hostmatrixstring_s, hostinclusionstring_s)
+ f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
+ eta = eta + f_d * CBACK * simpson(n) * xdts(n)
+ enddo
+ ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
+ endif
+
+
+!..Reflectivity contributed by melting graupel
+
+ if (L_qg(k) .and. L_qg(k_0) ) then
+ fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
+ eta = 0.d0
+ lamg = 1./ilamg(k)
+ do n = 1, nrbins
+ x = xam_g * xxDg(n)**xbm_g
+ call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
+ fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
+ CBACK, mixingrulestring_g, matrixstring_g, &
+ inclusionstring_g, hoststring_g, &
+ hostmatrixstring_g, hostinclusionstring_g)
+ f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
+ eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
+ enddo
+ ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
+ endif
+
+ enddo
+ endif
+
+ do k = kte, kts, -1
+ dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
+ enddo
+
+
+ end subroutine refl10cm_wsm6
+!+---+-----------------------------------------------------------------+
+
END MODULE module_mp_wsm6
</font>
</pre>