source: LMDZ6/trunk/libf/phylmd/cosp/dsd.f90 @ 5306

Last change on this file since 5306 was 5268, checked in by abarral, 8 days ago

.f90 <-> .F90 depending on cpp key use

  • 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
  • Property svn:keywords set to Id
File size: 10.5 KB
Line 
1!
2! $Id: dsd.f90 5268 2024-10-23 17:02:39Z fairhead $
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
Note: See TracBrowser for help on using the repository browser.