function iwc_expo_ZV_bvbzbm, Z, Vdq, bm, bv, bz, am, av, az, ba, rho_air

; assume inputs are in cgs units

; This function calculates the iwc.  The formulation is based
; on pg6 of notes.    This iwc is consistent with the assumed
; characteristics of the ice crystal as represented by the mass and fall speed power laws that have been optimized.

;parameter inputs are bv, bz, and bm.

; begin by choosing which value of the best number is appropriate.  Use a generic area relation.  This should only cause some
; error around the boundaries between coeffecients

dyn_visc=(1.7d-5)*(10.d)     ; the 1000/100 converts from kg/m*s to g/cm*s
kin_visc=dyn_visc/rho_air
g=981.d


;ab=(0.04394+0.06049)/2.d & bb=(0.970+0.831)/2.   ; this make the most sense from a physical perspective
ab=(0.04394) & bb=(0.970)
;ab=0.06049 &
;bb=0.8


; cv

term1=(gamma(6.+bv+bz))/(gamma(6.+bz))
term2=1.+(bv/(6.+bz))
cv=av*term1*term2



exp1=(bm)/(((bb*bm)+(2.*bb)-(bb*ba))-1.d)

const=2.*(0.19)*((6./(3.14159*0.92))^2)*((am))

; the exponent to vdq is term1

iwc=((gamma(bm))/(gamma(2.*bm)))*(Z/const)*((cv/vdq)^exp1)


return, iwc
end