| 1 | subroutine dsd(Q,Re,D,N,nsizes,dtype,rho_a,tc, & |
|---|
| 2 | dmin,dmax,apm,bpm,rho_c,p1,p2,p3,fc,scaled) |
|---|
| 3 | use array_lib |
|---|
| 4 | use math_lib |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | ! Purpose: |
|---|
| 8 | ! Create a discrete drop size distribution |
|---|
| 9 | ! Part of QuickBeam v1.03 by John Haynes |
|---|
| 10 | ! http://reef.atmos.colostate.edu/haynes/radarsim |
|---|
| 11 | ! |
|---|
| 12 | ! Inputs: |
|---|
| 13 | ! [Q] hydrometeor mixing ratio (g/kg) |
|---|
| 14 | ! [Re] Optional Effective Radius (microns). 0 = use default. |
|---|
| 15 | ! [D] discrete drop sizes (um) |
|---|
| 16 | ! [nsizes] number of elements of [D] |
|---|
| 17 | ! [dtype] distribution type |
|---|
| 18 | ! [rho_a] ambient air density (kg m^-3) |
|---|
| 19 | ! [tc] temperature (C) |
|---|
| 20 | ! [dmin] minimum size cutoff (um) |
|---|
| 21 | ! [dmax] maximum size cutoff (um) |
|---|
| 22 | ! [rho_c] alternate constant density (kg m^-3) |
|---|
| 23 | ! [p1],[p2],[p3] distribution parameters |
|---|
| 24 | ! |
|---|
| 25 | ! Input/Output: |
|---|
| 26 | ! [fc] scaling factor for the distribution |
|---|
| 27 | ! [scaled] has this hydrometeor type been scaled? |
|---|
| 28 | ! [apm] a parameter for mass (kg m^[-bpm]) |
|---|
| 29 | ! [bmp] b params for mass |
|---|
| 30 | ! |
|---|
| 31 | ! Outputs: |
|---|
| 32 | ! [N] discrete concentrations (cm^-3 um^-1) |
|---|
| 33 | ! or, for monodisperse, a constant (1/cm^3) |
|---|
| 34 | ! |
|---|
| 35 | ! Requires: |
|---|
| 36 | ! function infind |
|---|
| 37 | ! |
|---|
| 38 | ! Created: |
|---|
| 39 | ! 11/28/05 John Haynes (haynes@atmos.colostate.edu) |
|---|
| 40 | ! Modified: |
|---|
| 41 | ! 01/31/06 Port from IDL to Fortran 90 |
|---|
| 42 | ! 07/07/06 Rewritten for variable DSD's |
|---|
| 43 | ! 10/02/06 Rewritten using scaling factors (Roger Marchand and JMH) |
|---|
| 44 | |
|---|
| 45 | ! ----- INPUTS ----- |
|---|
| 46 | |
|---|
| 47 | integer*4, intent(in) :: nsizes |
|---|
| 48 | integer, intent(in) :: dtype |
|---|
| 49 | real*8, intent(in) :: Q,D(nsizes),rho_a,tc,dmin,dmax, & |
|---|
| 50 | rho_c,p1,p2,p3 |
|---|
| 51 | |
|---|
| 52 | ! ----- INPUT/OUTPUT ----- |
|---|
| 53 | |
|---|
| 54 | real*8, intent(inout) :: fc(nsizes),apm,bpm,Re |
|---|
| 55 | logical, intent(inout) :: scaled |
|---|
| 56 | |
|---|
| 57 | ! ----- OUTPUTS ----- |
|---|
| 58 | |
|---|
| 59 | real*8, intent(out) :: N(nsizes) |
|---|
| 60 | |
|---|
| 61 | ! ----- INTERNAL ----- |
|---|
| 62 | |
|---|
| 63 | real*8 :: & |
|---|
| 64 | N0,D0,vu,np,dm,ld, & ! gamma, exponential variables |
|---|
| 65 | dmin_mm,dmax_mm,ahp,bhp, & ! power law variables |
|---|
| 66 | rg,log_sigma_g, & ! lognormal variables |
|---|
| 67 | rho_e ! particle density (kg m^-3) |
|---|
| 68 | |
|---|
| 69 | real*8 :: tmp1, tmp2 |
|---|
| 70 | real*8 :: pi,rc |
|---|
| 71 | |
|---|
| 72 | integer k,lidx,uidx |
|---|
| 73 | |
|---|
| 74 | pi = acos(-1.0) |
|---|
| 75 | |
|---|
| 76 | ! // if density is constant, store equivalent values for apm and bpm |
|---|
| 77 | if ((rho_c > 0) .and. (apm < 0)) then |
|---|
| 78 | apm = (pi/6)*rho_c |
|---|
| 79 | bpm = 3. |
|---|
| 80 | endif |
|---|
| 81 | |
|---|
| 82 | select case(dtype) |
|---|
| 83 | |
|---|
| 84 | ! ---------------------------------------------------------! |
|---|
| 85 | ! // modified gamma ! |
|---|
| 86 | ! ---------------------------------------------------------! |
|---|
| 87 | ! :: N0 = total number concentration (m^-3) |
|---|
| 88 | ! :: np = fixed number concentration (kg^-1) |
|---|
| 89 | ! :: D0 = characteristic diameter (um) |
|---|
| 90 | ! :: dm = mean diameter (um) |
|---|
| 91 | ! :: vu = distribution width parameter |
|---|
| 92 | |
|---|
| 93 | case(1) |
|---|
| 94 | if (abs(p1+1) < 1E-8) then |
|---|
| 95 | |
|---|
| 96 | ! // D0, vu are given |
|---|
| 97 | vu = p3 |
|---|
| 98 | |
|---|
| 99 | if(Re.le.0) then |
|---|
| 100 | dm = p2 |
|---|
| 101 | D0 = gamma(vu)/gamma(vu+1)*dm |
|---|
| 102 | else |
|---|
| 103 | D0 = 2.0*Re*gamma(vu+2)/gamma(vu+3) |
|---|
| 104 | endif |
|---|
| 105 | |
|---|
| 106 | if (scaled .eqv. .false.) then |
|---|
| 107 | |
|---|
| 108 | fc = ( & |
|---|
| 109 | ((D*1E-6)**(vu-1)*exp(-1*D/D0)) / & |
|---|
| 110 | (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm)) & |
|---|
| 111 | ) * 1E-12 |
|---|
| 112 | scaled = .true. |
|---|
| 113 | |
|---|
| 114 | endif |
|---|
| 115 | |
|---|
| 116 | N = fc*rho_a*(Q*1E-3) |
|---|
| 117 | |
|---|
| 118 | elseif (abs(p2+1) < 1E-8) then |
|---|
| 119 | |
|---|
| 120 | ! // N0, vu are given |
|---|
| 121 | np = p1 |
|---|
| 122 | vu = p3 |
|---|
| 123 | tmp1 = (Q*1E-3)**(1./bpm) |
|---|
| 124 | |
|---|
| 125 | if (scaled .eqv. .false.) then |
|---|
| 126 | |
|---|
| 127 | fc = (D*1E-6 / (gamma(vu)/(apm*np*gamma(vu+bpm)))** & |
|---|
| 128 | (1./bpm))**vu |
|---|
| 129 | |
|---|
| 130 | scaled = .true. |
|---|
| 131 | |
|---|
| 132 | endif |
|---|
| 133 | |
|---|
| 134 | N = ( & |
|---|
| 135 | (rho_a*np*fc*(D*1E-6)**(-1.))/(gamma(vu)*tmp1**vu) * & |
|---|
| 136 | exp(-1.*fc**(1./vu)/tmp1) & |
|---|
| 137 | ) * 1E-12 |
|---|
| 138 | |
|---|
| 139 | else |
|---|
| 140 | |
|---|
| 141 | ! // vu isn't given |
|---|
| 142 | print *, 'Error: Must specify a value for vu' |
|---|
| 143 | stop |
|---|
| 144 | |
|---|
| 145 | endif |
|---|
| 146 | |
|---|
| 147 | ! ---------------------------------------------------------! |
|---|
| 148 | ! // exponential ! |
|---|
| 149 | ! ---------------------------------------------------------! |
|---|
| 150 | ! :: N0 = intercept parameter (m^-4) |
|---|
| 151 | ! :: ld = slope parameter (um) |
|---|
| 152 | |
|---|
| 153 | case(2) |
|---|
| 154 | if (abs(p1+1) > 1E-8) then |
|---|
| 155 | |
|---|
| 156 | ! // N0 has been specified, determine ld |
|---|
| 157 | N0 = p1 |
|---|
| 158 | |
|---|
| 159 | if(Re>0) then |
|---|
| 160 | |
|---|
| 161 | ! if Re is set and No is set than the distribution is fully defined. |
|---|
| 162 | ! so we assume Re and No have already been chosen consistant with |
|---|
| 163 | ! the water content, Q. |
|---|
| 164 | |
|---|
| 165 | ! print *,'using Re pass ...' |
|---|
| 166 | |
|---|
| 167 | ld = 1.5/Re ! units 1/um |
|---|
| 168 | |
|---|
| 169 | N = ( & |
|---|
| 170 | N0*exp(-1*ld*D) & |
|---|
| 171 | ) * 1E-12 |
|---|
| 172 | |
|---|
| 173 | else |
|---|
| 174 | |
|---|
| 175 | tmp1 = 1./(1.+bpm) |
|---|
| 176 | |
|---|
| 177 | if (scaled .eqv. .false.) then |
|---|
| 178 | fc = ((apm*gamma(1.+bpm)*N0)**tmp1)*(D*1E-6) |
|---|
| 179 | scaled = .true. |
|---|
| 180 | |
|---|
| 181 | endif |
|---|
| 182 | |
|---|
| 183 | N = ( & |
|---|
| 184 | N0*exp(-1.*fc*(1./(rho_a*Q*1E-3))**tmp1) & |
|---|
| 185 | ) * 1E-12 |
|---|
| 186 | |
|---|
| 187 | endif |
|---|
| 188 | |
|---|
| 189 | elseif (abs(p2+1) > 1E-8) then |
|---|
| 190 | |
|---|
| 191 | ! // ld has been specified, determine N0 |
|---|
| 192 | ld = p2 |
|---|
| 193 | |
|---|
| 194 | if (scaled .eqv. .false.) then |
|---|
| 195 | |
|---|
| 196 | fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))* & |
|---|
| 197 | exp(-1.*(ld*1E6)*(D*1E-6))*1E-12 |
|---|
| 198 | scaled = .true. |
|---|
| 199 | |
|---|
| 200 | endif |
|---|
| 201 | |
|---|
| 202 | N = fc*rho_a*(Q*1E-3) |
|---|
| 203 | |
|---|
| 204 | else |
|---|
| 205 | |
|---|
| 206 | ! // ld will be determined from temperature, then N0 follows |
|---|
| 207 | ld = 1220*10.**(-0.0245*tc)*1E-6 |
|---|
| 208 | N0 = ((ld*1E6)**(1+bpm)*Q*1E-3*rho_a)/(apm*gamma(1+bpm)) |
|---|
| 209 | |
|---|
| 210 | N = ( & |
|---|
| 211 | N0*exp(-1*ld*D) & |
|---|
| 212 | ) * 1E-12 |
|---|
| 213 | |
|---|
| 214 | endif |
|---|
| 215 | |
|---|
| 216 | ! ---------------------------------------------------------! |
|---|
| 217 | ! // power law ! |
|---|
| 218 | ! ---------------------------------------------------------! |
|---|
| 219 | ! :: ahp = Ar parameter (m^-4 mm^-bhp) |
|---|
| 220 | ! :: bhp = br parameter |
|---|
| 221 | ! :: dmin_mm = lower bound (mm) |
|---|
| 222 | ! :: dmax_mm = upper bound (mm) |
|---|
| 223 | |
|---|
| 224 | case(3) |
|---|
| 225 | |
|---|
| 226 | ! :: br parameter |
|---|
| 227 | if (abs(p1+2) < 1E-8) then |
|---|
| 228 | ! :: if p1=-2, bhp is parameterized according to Ryan (2000), |
|---|
| 229 | ! :: applicatable to cirrus clouds |
|---|
| 230 | if (tc < -30) then |
|---|
| 231 | bhp = -1.75+0.09*((tc+273)-243.16) |
|---|
| 232 | elseif ((tc >= -30) .and. (tc < -9)) then |
|---|
| 233 | bhp = -3.25-0.06*((tc+273)-265.66) |
|---|
| 234 | else |
|---|
| 235 | bhp = -2.15 |
|---|
| 236 | endif |
|---|
| 237 | elseif (abs(p1+3) < 1E-8) then |
|---|
| 238 | ! :: if p1=-3, bhp is parameterized according to Ryan (2000), |
|---|
| 239 | ! :: applicable to frontal clouds |
|---|
| 240 | if (tc < -35) then |
|---|
| 241 | bhp = -1.75+0.09*((tc+273)-243.16) |
|---|
| 242 | elseif ((tc >= -35) .and. (tc < -17.5)) then |
|---|
| 243 | bhp = -2.65+0.09*((tc+273)-255.66) |
|---|
| 244 | elseif ((tc >= -17.5) .and. (tc < -9)) then |
|---|
| 245 | bhp = -3.25-0.06*((tc+273)-265.66) |
|---|
| 246 | else |
|---|
| 247 | bhp = -2.15 |
|---|
| 248 | endif |
|---|
| 249 | else |
|---|
| 250 | ! :: otherwise the specified value is used |
|---|
| 251 | bhp = p1 |
|---|
| 252 | endif |
|---|
| 253 | |
|---|
| 254 | ! :: Ar parameter |
|---|
| 255 | dmin_mm = dmin*1E-3 |
|---|
| 256 | dmax_mm = dmax*1E-3 |
|---|
| 257 | |
|---|
| 258 | ! :: commented lines are original method with constant density |
|---|
| 259 | ! rc = 500. ! (kg/m^3) |
|---|
| 260 | ! tmp1 = 6*rho_a*(bhp+4) |
|---|
| 261 | ! tmp2 = pi*rc*(dmax_mm**(bhp+4))*(1-(dmin_mm/dmax_mm)**(bhp+4)) |
|---|
| 262 | ! ahp = (Q*1E-3)*1E12*tmp1/tmp2 |
|---|
| 263 | |
|---|
| 264 | ! :: new method is more consistent with the rest of the distributions |
|---|
| 265 | ! :: and allows density to vary with particle size |
|---|
| 266 | tmp1 = rho_a*(Q*1E-3)*(bhp+bpm+1) |
|---|
| 267 | tmp2 = apm*(dmax_mm**bhp*dmax**(bpm+1)-dmin_mm**bhp*dmin**(bpm+1)) |
|---|
| 268 | ahp = tmp1/tmp2 * 1E24 |
|---|
| 269 | ! ahp = tmp1/tmp2 |
|---|
| 270 | |
|---|
| 271 | lidx = infind(D,dmin) |
|---|
| 272 | uidx = infind(D,dmax) |
|---|
| 273 | do k=lidx,uidx |
|---|
| 274 | |
|---|
| 275 | N(k) = ( & |
|---|
| 276 | ahp*(D(k)*1E-3)**bhp & |
|---|
| 277 | ) * 1E-12 |
|---|
| 278 | |
|---|
| 279 | enddo |
|---|
| 280 | |
|---|
| 281 | ! print *,'test=',ahp,bhp,ahp/(rho_a*Q),D(100),N(100),bpm,dmin_mm,dmax_mm |
|---|
| 282 | |
|---|
| 283 | ! ---------------------------------------------------------! |
|---|
| 284 | ! // monodisperse ! |
|---|
| 285 | ! ---------------------------------------------------------! |
|---|
| 286 | ! :: D0 = particle diameter (um) |
|---|
| 287 | |
|---|
| 288 | case(4) |
|---|
| 289 | |
|---|
| 290 | if (scaled .eqv. .false.) then |
|---|
| 291 | |
|---|
| 292 | D0 = p1 |
|---|
| 293 | rho_e = (6/pi)*apm*(D0*1E-6)**(bpm-3) |
|---|
| 294 | fc(1) = (6./(pi*D0**3*rho_e))*1E12 |
|---|
| 295 | scaled = .true. |
|---|
| 296 | |
|---|
| 297 | endif |
|---|
| 298 | |
|---|
| 299 | N(1) = fc(1)*rho_a*(Q*1E-3) |
|---|
| 300 | |
|---|
| 301 | ! ---------------------------------------------------------! |
|---|
| 302 | ! // lognormal ! |
|---|
| 303 | ! ---------------------------------------------------------! |
|---|
| 304 | ! :: N0 = total number concentration (m^-3) |
|---|
| 305 | ! :: np = fixed number concentration (kg^-1) |
|---|
| 306 | ! :: rg = mean radius (um) |
|---|
| 307 | ! :: log_sigma_g = ln(geometric standard deviation) |
|---|
| 308 | |
|---|
| 309 | case(5) |
|---|
| 310 | if (abs(p1+1) < 1E-8) then |
|---|
| 311 | |
|---|
| 312 | ! // rg, log_sigma_g are given |
|---|
| 313 | log_sigma_g = p3 |
|---|
| 314 | tmp2 = (bpm*log_sigma_g)**2. |
|---|
| 315 | if(Re.le.0) then |
|---|
| 316 | rg = p2 |
|---|
| 317 | else |
|---|
| 318 | rg =Re*exp(-2.5*(log_sigma_g**2)) |
|---|
| 319 | endif |
|---|
| 320 | |
|---|
| 321 | if (scaled .eqv. .false.) then |
|---|
| 322 | |
|---|
| 323 | fc = 0.5 * ( & |
|---|
| 324 | (1./((2.*rg*1E-6)**(bpm)*apm*(2.*pi)**(0.5) * & |
|---|
| 325 | log_sigma_g*D*0.5*1E-6)) * & |
|---|
| 326 | exp(-0.5*((log(0.5*D/rg)/log_sigma_g)**2.+tmp2)) & |
|---|
| 327 | ) * 1E-12 |
|---|
| 328 | scaled = .true. |
|---|
| 329 | |
|---|
| 330 | endif |
|---|
| 331 | |
|---|
| 332 | N = fc*rho_a*(Q*1E-3) |
|---|
| 333 | |
|---|
| 334 | elseif (abs(p2+1) < 1E-8) then |
|---|
| 335 | |
|---|
| 336 | ! // N0, log_sigma_g are given |
|---|
| 337 | Np = p1 |
|---|
| 338 | log_sigma_g = p3 |
|---|
| 339 | N0 = np*rho_a |
|---|
| 340 | tmp1 = (rho_a*(Q*1E-3))/(2.**bpm*apm*N0) |
|---|
| 341 | tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2. |
|---|
| 342 | rg = ((tmp1/tmp2)**(1/bpm))*1E6 |
|---|
| 343 | |
|---|
| 344 | N = 0.5*( & |
|---|
| 345 | N0 / ((2.*pi)**(0.5)*log_sigma_g*D*0.5*1E-6) * & |
|---|
| 346 | exp((-0.5*(log(0.5*D/rg)/log_sigma_g)**2.)) & |
|---|
| 347 | ) * 1E-12 |
|---|
| 348 | |
|---|
| 349 | else |
|---|
| 350 | |
|---|
| 351 | ! // vu isn't given |
|---|
| 352 | print *, 'Error: Must specify a value for sigma_g' |
|---|
| 353 | stop |
|---|
| 354 | |
|---|
| 355 | endif |
|---|
| 356 | |
|---|
| 357 | end select |
|---|
| 358 | |
|---|
| 359 | end subroutine dsd |
|---|
