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