source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/cosp/dsd.F90 @ 4003

Last change on this file since 4003 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 10.4 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.