The Kozu and Iguchi's DSD model includes parameters in 
not only R-Ze and K-Ze relations but also WC-Ze relations.
Using the coefficients of WC-Ze relations, we can calculate 
(a) water content (or mixing ratio of precipitation), and 
(b) mean terminal velocity of precipitation.  

R [mm hr^-1] is rainfall rate.
R*10e-3/3600 [m s^-1]
RhoL*R*10e-3/3600 [kg m^-2 s^-1]
  where RhoL=10e3 [kg m^-3] is water density.
R/(3600*Vt) [kg m^-3], ---- (1)
  where Vt [m s^-1] is terminal velocity

While, WC [g m^-3] is water content.
WC*10e-3 [kg m^-3]     ---- (2)

From (1) and (2), we can obtain the relations of R, WC, and Vt.
R=3.6*Vt*WC            ---- (3)
When the mixing ratio of precipitation qP [kg kg^-1] is used,
R=3.6*Vt*Rho*qP        ---- (4)
where Rho is the air-density in each height. The qP seems to 
include qR, qS, and qG (mixing ratio of rain, snow, grauple/hail).
Since the U.S. standard atmosphere has been used in 2A25 for vratio, 
the Rho in (4) is also U.S. standard value.
vratio=(Rho0/Rho)^0.4

The follows are extracts from my program source.  In this version, 
the coefficients of WC-Ze relations are used when epsilon=1. 
However, I am expecting adjusted WC-Ze parameters by epsilon 
in 5 nodes will be included in 2A25 version 6.  The new raintype 
flag, which is "Anvil", is also expected (see Fig. 5)

/* ------- WC-Ze relations ---------- */
static float ZWParmAlpha[3][5]={
  {0.00620868,0.00537114,0.00391752,0.00391752,0.00444452}, /* others */
  {0.00383613,0.00325046,0.00074301,0.00199806,0.00223787}, /* strat */
  {0.00620868,0.00391752,0.00391752,0.00391752,0.00444452}  /* conv */
};
static float ZWParmBeta[3][5]={
  {0.68902,0.68004,0.57855,0.57855,0.56232},                /* others */
  {0.71266,0.70472,0.66564,0.61342,0.59658},                /* strat */
  {0.68902,0.57855,0.57855,0.57855,0.56232}                 /* conv */
};
static float ParmRho[21]=     /* Rho in the US Standard Atmos. */
  {0.08891,0.10400,0.12165,0.14230,0.16647,0.19476,0.22786,
   0.26660,0.31194,0.36480,0.41351,0.46706,0.52579,0.59002,
   0.66011,0.73643,0.81935,0.90925,1.0066,1.1117,1.2250};

rho0=ParmRho[20];
vratio[kk]=pow(rho0/rho[kk],0.4);

fraction = (float)(k-read_data.ZRParmNode[j][l])/
           (float)(read_data.ZRParmNode[j][lp1]-read_data.ZRParmNode[j][l]);
zra[kk][n][m]=read_data.ZRParmA[j][l]+
           (read_data.ZRParmA[j][lp1]-read_data.ZRParmA[j][l])*fraction;
zrb[kk][n][m]=read_data.ZRParmB[j][l]+
           (read_data.ZRParmB[j][lp1]-read_data.ZRParmB[j][l])*fraction;
zwa[kk][n][m] = ZWParmAlpha[irt][l] +
           (ZWParmAlpha[irt][lp1]-ZWParmAlpha[irt][l])*fraction;
zwb[kk][n][m] = ZWParmBeta[irt][l] +
           (ZWParmBeta[irt][lp1]-ZWParmBeta[irt][l])*fraction;

xze  = (float)read_data.correctZFactor[j][k];
xz = pow(10.0,xze/10.0);
xrain2 = zra[kk][n][m] * pow(xz,zrb[kk][n][m]) / vratio[kk];
xwc = zwa[kk][n][m] * pow(xz,zwb[kk][n][m]);
xqp = xwc/rho[kk];
xvt = xrain2/(3.6*xwc);
xvt0 = xvt/vratio[kk];