1 | ! |
---|
2 | ! $Id: dsd.F90 4446 2023-03-01 17:22:39Z evignon $ |
---|
3 | ! |
---|
4 | subroutine dsd(Q,Re_,Np,D,N,nsizes,dtype,rho_a,tk, & |
---|
5 | dmin,dmax,apm,bpm,rho_c,p1,p2,p3) |
---|
6 | use array_lib |
---|
7 | use math_lib |
---|
8 | implicit none |
---|
9 | |
---|
10 | ! Purpose: |
---|
11 | ! Create a discrete drop size distribution |
---|
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 |
---|
21 | ! http://reef.atmos.colostate.edu/haynes/radarsim |
---|
22 | ! |
---|
23 | ! Inputs: |
---|
24 | ! |
---|
25 | ! [Q] hydrometeor mixing ratio (g/kg) |
---|
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). |
---|
29 | ! [nsizes] number of elements of [D] |
---|
30 | ! |
---|
31 | ! [dtype] distribution type |
---|
32 | ! [rho_a] ambient air density (kg m^-3) |
---|
33 | ! [tk] temperature (K) |
---|
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 |
---|
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). |
---|
57 | |
---|
58 | ! ----- INPUTS ----- |
---|
59 | |
---|
60 | integer, intent(in) :: nsizes |
---|
61 | integer, intent(in) :: dtype |
---|
62 | real*8, intent(in) :: Q,Re_,Np,D(nsizes) |
---|
63 | real*8, intent(in) :: rho_a,tk,dmin,dmax,rho_c,p1,p2,p3 |
---|
64 | |
---|
65 | real*8, intent(inout) :: apm,bpm |
---|
66 | |
---|
67 | ! ----- OUTPUTS ----- |
---|
68 | |
---|
69 | real*8, intent(out) :: N(nsizes) |
---|
70 | |
---|
71 | ! ----- INTERNAL ----- |
---|
72 | |
---|
73 | real*8 :: fc(nsizes) |
---|
74 | |
---|
75 | real*8 :: & |
---|
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) |
---|
80 | |
---|
81 | real*8 :: tmp1, tmp2 |
---|
82 | real*8 :: pi,rc,tc |
---|
83 | real*8 :: Re |
---|
84 | |
---|
85 | integer k,lidx,uidx |
---|
86 | |
---|
87 | Re = Re_ |
---|
88 | |
---|
89 | tc = tk - 273.15 |
---|
90 | pi = acos(-1.0) |
---|
91 | |
---|
92 | ! // if density is constant, store equivalent values for apm and bpm |
---|
93 | if ((rho_c > 0) .and. (apm < 0)) then |
---|
94 | apm = (pi/6)*rho_c |
---|
95 | bpm = 3. |
---|
96 | endif |
---|
97 | |
---|
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 | |
---|
109 | select case(dtype) |
---|
110 | |
---|
111 | ! ---------------------------------------------------------! |
---|
112 | ! // modified gamma ! |
---|
113 | ! ---------------------------------------------------------! |
---|
114 | ! :: np = total number concentration |
---|
115 | ! :: D0 = characteristic diameter (um) |
---|
116 | ! :: dm = mean diameter (um) - first moment over zeroth moment |
---|
117 | ! :: vu = distribution width parameter |
---|
118 | |
---|
119 | case(1) |
---|
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 |
---|
135 | |
---|
136 | ! vu is fixed in hp structure |
---|
137 | vu = p3 |
---|
138 | |
---|
139 | else |
---|
140 | |
---|
141 | ! vu isn't specified |
---|
142 | |
---|
143 | print *, 'Error: Must specify a value for vu for Modified Gamma distribution' |
---|
144 | stop |
---|
145 | |
---|
146 | endif |
---|
147 | |
---|
148 | if(Re>0) then |
---|
149 | |
---|
150 | D0 = 2.0*Re*gamma(vu+2)/gamma(vu+3) |
---|
151 | |
---|
152 | fc = ( & |
---|
153 | ((D*1E-6)**(vu-1)*exp(-1*D/D0)) / & |
---|
154 | (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm)) & |
---|
155 | ) * 1E-12 |
---|
156 | |
---|
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 |
---|
163 | |
---|
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 |
---|
172 | |
---|
173 | local_np = p1 ! total number concentration / pa check |
---|
174 | |
---|
175 | tmp1 = (Q*1E-3)**(1./bpm) |
---|
176 | |
---|
177 | fc = (D*1E-6 / (gamma(vu)/(apm*local_np*gamma(vu+bpm)))** & |
---|
178 | (1./bpm))**vu |
---|
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 |
---|
184 | |
---|
185 | else |
---|
186 | |
---|
187 | print *, 'Error: No default value for Dm or Np provided! ' |
---|
188 | stop |
---|
189 | |
---|
190 | endif |
---|
191 | |
---|
192 | |
---|
193 | ! ---------------------------------------------------------! |
---|
194 | ! // exponential ! |
---|
195 | ! ---------------------------------------------------------! |
---|
196 | ! :: N0 = intercept parameter (m^-4) |
---|
197 | ! :: ld = slope parameter (um) |
---|
198 | |
---|
199 | case(2) |
---|
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 |
---|
211 | |
---|
212 | ! use N0 default value |
---|
213 | |
---|
214 | N0 = p1 |
---|
215 | |
---|
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 |
---|
223 | |
---|
224 | elseif (abs(p2+1) > 1E-8) then |
---|
225 | |
---|
226 | ! used default value for lambda |
---|
227 | ld = p2 |
---|
228 | |
---|
229 | fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))* & |
---|
230 | exp(-1.*(ld*1E6)*(D*1E-6))*1E-12 |
---|
231 | |
---|
232 | N = fc*rho_a*(Q*1E-3) |
---|
233 | |
---|
234 | else |
---|
235 | |
---|
236 | ! ld "parameterized" from temperature (carry over from original Quickbeam). |
---|
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) |
---|
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 |
---|
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 |
---|
299 | ! rc = 500. ! (kg/m^3) |
---|
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 | |
---|
315 | N(k) = ( & |
---|
316 | ahp*(D(k)*1E-3)**bhp & |
---|
317 | ) * 1E-12 |
---|
318 | |
---|
319 | enddo |
---|
320 | |
---|
321 | ! print *,'test=',ahp,bhp,ahp/(rho_a*Q),D(100),N(100),bpm,dmin_mm,dmax_mm |
---|
322 | |
---|
323 | ! ---------------------------------------------------------! |
---|
324 | ! // monodisperse ! |
---|
325 | ! ---------------------------------------------------------! |
---|
326 | ! :: D0 = particle diameter (um) |
---|
327 | |
---|
328 | case(4) |
---|
329 | |
---|
330 | if (Re>0) then |
---|
331 | D0 = Re |
---|
332 | else |
---|
333 | D0 = p1 |
---|
334 | endif |
---|
335 | |
---|
336 | rho_e = (6/pi)*apm*(D0*1E-6)**(bpm-3) |
---|
337 | fc(1) = (6./(pi*D0**3*rho_e))*1E12 |
---|
338 | N(1) = fc(1)*rho_a*(Q*1E-3) |
---|
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) |
---|
349 | if (abs(p1+1) < 1E-8 .or. Re>0 ) then |
---|
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 |
---|
355 | rg = p2 |
---|
356 | else |
---|
357 | rg =Re*exp(-2.5*(log_sigma_g**2)) |
---|
358 | endif |
---|
359 | |
---|
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 | |
---|
366 | N = fc*rho_a*(Q*1E-3) |
---|
367 | |
---|
368 | elseif (abs(p2+1) < 1E-8 .or. Np>0) then |
---|
369 | |
---|
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 | |
---|
377 | log_sigma_g = p3 |
---|
378 | N0 = local_np*rho_a |
---|
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) * & |
---|
385 | exp((-0.5*(log(0.5*D/rg)/log_sigma_g)**2.)) & |
---|
386 | ) * 1E-12 |
---|
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 |
---|