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