source: LMDZ5/trunk/libf/cosp/dsd.F90 @ 2216

Last change on this file since 2216 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 8.9 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.