source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/cosp/calc_Re.F90 @ 3373

Last change on this file since 3373 was 2432, checked in by idelkadi, 9 years ago

Correction : rajout de fichiers manquants pour les simulateurs Modis, Cloudsat, Misr et Isccp

File size: 7.1 KB
Line 
1  subroutine calc_Re(Q,Np,rho_a, &
2             dtype,dmin,dmax,apm,bpm,rho_c,p1,p2,p3, &
3             Re)
4  use math_lib
5  implicit none
6
7! Purpose:
8!   Calculates Effective Radius (1/2 distribution 3rd moment / 2nd moment).
9!
10!   For some distribution types, the total number concentration (per kg), Np
11!   may be optionally specified.   Should be set to zero, otherwise.
12!
13!   Roj Marchand July 2010
14
15
16! Inputs:
17!
18!   [Q]        hydrometeor mixing ratio (g/kg)  ! not needed for some distribution types
19!   [Np]       Optional Total number concentration (per kg).  0 = use defaults (p1, p2, p3)
20!
21!   [rho_a]    ambient air density (kg m^-3)   
22!
23!   Distribution parameters as per quickbeam documentation.
24!   [dtype]    distribution type
25!   [dmin]     minimum size cutoff (um)
26!   [dmax]     maximum size cutoff (um)
27!   [apm]      a parameter for mass (kg m^[-bpm])
28!   [bmp]      b params for mass
29!   [p1],[p2],[p3]  distribution parameters
30!
31!
32! Outputs:
33!   [Re]       Effective radius, 1/2 the 3rd moment/2nd moment (um)
34!
35! Created:
36!   July 2010  Roj Marchand
37!
38
39 
40! ----- INPUTS ----- 
41 
42  real*8, intent(in) :: Q,Np,rho_a
43 
44  integer, intent(in):: dtype
45  real*8, intent(in) :: dmin,dmax,rho_c,p1,p2,p3
46   
47  real*8, intent(inout) :: apm,bpm 
48   
49! ----- OUTPUTS -----
50
51  real*8, intent(out) :: Re
52 
53! ----- INTERNAL -----
54 
55  integer :: local_dtype
56  real*8  :: local_p3,local_Np
57
58  real*8 :: pi, &
59  N0,D0,vu,dm,ld, &         ! gamma, exponential variables
60  rg,log_sigma_g
61 
62  real*8 :: tmp1,tmp2
63
64  pi = acos(-1.0)
65
66  ! // if density is constant, set equivalent values for apm and bpm
67  if ((rho_c > 0) .and. (apm < 0)) then
68    apm = (pi/6)*rho_c
69    bpm = 3.
70  endif
71
72  ! Exponential is same as modified gamma with vu =1
73  ! if Np is specified then we will just treat as modified gamma
74  if(dtype.eq.2 .and. Np>0) then
75    local_dtype=1;
76    local_p3=1;
77  else
78    local_dtype=dtype;
79    local_p3=p3;
80  endif
81 
82  select case(local_dtype)
83 
84! ---------------------------------------------------------!
85! // modified gamma                                        !
86! ---------------------------------------------------------!
87! :: Np = total number concentration (1/kg) = Nt / rho_a
88! :: D0 = characteristic diameter (um)
89! :: dm = mean diameter (um) - first moment over zeroth moment
90! :: vu = distribution width parameter
91
92  case(1) 
93 
94    if( abs(local_p3+2) < 1E-8) then
95 
96    if(Np>1E-30) then
97        ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1)
98        ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane
99        vu = (1/(0.2714 + 0.00057145*Np*rho_a*1E-6))**2 ! units of Nt = Np*rhoa = #/cm^3
100    else
101        print *, 'Error: Must specify a value for Np in each volume', &
102             ' with Morrison/Martin Scheme.'
103            stop   
104    endif
105   
106    elseif (abs(local_p3+1) > 1E-8) then
107
108      ! vu is fixed in hp structure 
109      vu = local_p3
110
111    else
112
113      ! vu isn't specified
114     
115      print *, 'Error: Must specify a value for vu for Modified Gamma distribution'
116      stop   
117     
118    endif
119   
120
121    if( Np.eq.0 .and. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default
122     
123        dm = p2             ! by definition, should have units of microns
124    D0 = gamma(vu)/gamma(vu+1)*dm
125       
126    else   ! use value of Np
127       
128        if(Np.eq.0) then
129       
130            if( abs(p1+1) > 1E-8 ) then  !   use default number concentration
131           
132                local_Np = p1 ! total number concentration / pa --- units kg^-1
133            else
134            print *, 'Error: Must specify Np or default value ', &
135                 '(p1=Dm [um] or p2=Np [1/kg]) for ', &
136                 'Modified Gamma distribution'
137                stop
138            endif
139        else
140            local_Np=Np;   
141    endif
142   
143    D0 = 1E6 * ( Q*1E-3*gamma(vu)/(apm*local_Np*gamma(vu+bpm)) )**(1/bpm)  ! units = microns
144
145    endif 
146     
147    Re = 0.5*D0*gamma(vu+3)/gamma(vu+2)
148   
149   
150! ---------------------------------------------------------!
151! // exponential                                           !
152! ---------------------------------------------------------!
153! :: N0 = intercept parameter (m^-4)
154! :: ld = slope parameter (um)
155
156  case(2)
157 
158    ! Np not specified (see if statement above)
159 
160    if((abs(p1+1) > 1E-8) ) then   ! N0 has been specified, determine ld
161   
162        N0 = p1
163    tmp1 = 1./(1.+bpm)
164    ld = ((apm*gamma(1.+bpm)*N0)/(rho_a*Q*1E-3))**tmp1
165    ld = ld/1E6                     ! set units to microns^-1
166       
167    elseif (abs(p2+1) > 1E-8) then  ! lambda=ld has been specified as default
168
169        ld = p2     ! should have units of microns^-1
170   
171    else
172   
173    print *, 'Error: Must specify Np or default value ', &
174         '(p1=No or p2=lambda) for Exponential distribution'
175        stop
176       
177    endif
178
179    Re = 1.5/ld ;
180 
181! ---------------------------------------------------------!
182! // power law                                             !
183! ---------------------------------------------------------!
184! :: ahp = Ar parameter (m^-4 mm^-bhp)
185! :: bhp = br parameter
186! :: dmin_mm = lower bound (mm)
187! :: dmax_mm = upper bound (mm)
188
189  case(3)
190
191    Re=0;  ! Not supporting LUT approach for power-law ...
192
193    if(Np>0) then
194   
195        print *, 'Variable Np not supported for ', &
196         'Power Law distribution'
197        stop
198    endif
199! ---------------------------------------------------------!
200! // monodisperse                                          !
201! ---------------------------------------------------------!
202! :: D0 = particle diameter (um) == Re
203
204  case(4)
205 
206        Re = p1
207   
208        if(Np>0) then
209        print *, 'Variable Np not supported for ', &
210         'Monodispersed distribution'
211        stop
212    endif
213   
214! ---------------------------------------------------------!
215! // lognormal                                             !
216! ---------------------------------------------------------!
217! :: N0 = total number concentration (m^-3)
218! :: np = fixed number concentration (kg^-1)
219! :: rg = mean radius (um)
220! :: log_sigma_g = ln(geometric standard deviation)
221
222  case(5)
223 
224        if( abs(local_p3+1) > 1E-8 ) then
225
226            !set natural log width
227            log_sigma_g = local_p3
228        else
229            print *, 'Error: Must specify a value for sigma_g ', &
230             'when using a Log-Normal distribution'
231            stop
232        endif
233     
234    ! get rg ...
235    if( Np.eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
236   
237            rg = p2
238             
239        else
240            if(Np>0) then
241           
242                local_Np=Np;
243               
244            elseif(abs(p2+1) < 1E-8) then
245           
246                local_Np=p1
247            else
248                print *, 'Error: Must specify Np or default value ', &
249                 '(p2=Rg or p1=Np) for Log-Normal distribution'
250            endif
251           
252            log_sigma_g = p3
253            tmp1 = (Q*1E-3)/(2.**bpm*apm*local_Np)
254            tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2.     
255            rg = ((tmp1/tmp2)**(1/bpm))*1E6
256        endif
257               
258        Re = rg*exp(+2.5*(log_sigma_g**2))
259       
260  end select
261 
262  end subroutine calc_Re
Note: See TracBrowser for help on using the repository browser.