Line data Source code
1 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 : ! Copyright (c) 2015, Regents of the University of Colorado
3 : ! All rights reserved.
4 : !
5 : ! Redistribution and use in source and binary forms, with or without modification, are
6 : ! permitted provided that the following conditions are met:
7 : !
8 : ! 1. Redistributions of source code must retain the above copyright notice, this list of
9 : ! conditions and the following disclaimer.
10 : !
11 : ! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 : ! of conditions and the following disclaimer in the documentation and/or other
13 : ! materials provided with the distribution.
14 : !
15 : ! 3. Neither the name of the copyright holder nor the names of its contributors may be
16 : ! used to endorse or promote products derived from this software without specific prior
17 : ! written permission.
18 : !
19 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 : ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 : ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24 : ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 : ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 : ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 : ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 : !
29 : ! History:
30 : ! May 2015: Dustin Swales - Initial version
31 : !
32 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 : module mod_quickbeam_optics
34 : USE COSP_KINDS, ONLY: wp,dp
35 : USE array_lib, ONLY: infind
36 : USE math_lib, ONLY: path_integral,avint,gamma
37 : USE optics_lib, ONLY: m_wat,m_ice,MieInt
38 : USE cosp_math_constants, ONLY: pi
39 : USE cosp_phys_constants, ONLY: rhoice
40 : use quickbeam, ONLY: radar_cfg,dmin,dmax,Re_BIN_LENGTH, &
41 : Re_MAX_BIN,nRe_types,nd,maxhclass,save_scale_LUTs
42 : use mod_cosp_config, ONLY: N_HYDRO
43 : use mod_cosp_error, ONLY: errorMessage
44 : implicit none
45 :
46 : ! Derived type for particle size distribution
47 : TYPE size_distribution
48 : real(wp),dimension(maxhclass) :: p1,p2,p3,dmin,dmax,apm,bpm,rho
49 : integer, dimension(maxhclass) :: dtype,phase
50 : END TYPE size_distribution
51 :
52 : ! Parameters
53 : integer,parameter :: & !
54 : cnt_liq = 19, & ! Liquid temperature count
55 : cnt_ice = 20 ! Lce temperature count
56 :
57 : ! Initialization variables
58 : real(wp),dimension(cnt_ice) :: mt_tti
59 : real(wp),dimension(cnt_liq) :: mt_ttl
60 : real(wp),dimension(nd) :: D
61 :
62 : contains
63 : ! ######################################################################################
64 : ! SUBROUTINE quickbeam_optics_init
65 : ! ######################################################################################
66 1536 : subroutine quickbeam_optics_init()
67 : integer :: j
68 :
69 1536 : mt_tti = (/ ((j-1)*5-90 + 273.15, j = 1, cnt_ice) /)
70 1536 : mt_ttl = (/ ((j-1)*5-60 + 273.15, j = 1, cnt_liq) /)
71 1536 : D(1) = dmin
72 130560 : do j=2,nd
73 130560 : D(j) = D(j-1)*exp((log(dmax)-log(dmin))/(nd-1))
74 : enddo
75 1536 : end subroutine quickbeam_optics_init
76 :
77 : ! ######################################################################################
78 : ! SUBROUTINE QUICKBEAM_OPTICS
79 : ! ######################################################################################
80 92880 : subroutine quickbeam_optics(sd, rcfg, nprof, ngate, undef, hm_matrix, re_matrix, &
81 92880 : Np_matrix, p_matrix, t_matrix, sh_matrix,z_vol,kr_vol)
82 :
83 : ! INPUTS
84 : type(size_distribution),intent(inout) :: &
85 : sd !
86 : type(radar_cfg),intent(inout) :: &
87 : rcfg !
88 : integer,intent(in) :: &
89 : nprof, & ! Number of hydrometeor profiles
90 : ngate ! Number of vertical layers
91 : real(wp),intent(in) :: &
92 : undef ! Missing data value
93 : real(wp),intent(in),dimension(nprof,ngate) :: &
94 : p_matrix, & ! Pressure profile (hPa)
95 : t_matrix, & ! Temperature profile (K)
96 : sh_matrix ! Specific humidity profile (%) -- only needed if gaseous aborption calculated.
97 : real(wp),intent(in),dimension(nprof,ngate,rcfg%nhclass) :: &
98 : re_matrix, & ! Table of hydrometeor effective radii. 0 ==> use defaults. (units=microns)
99 : hm_matrix ! Table of hydrometeor mixing ratios (g/kg)
100 : real(wp),intent(inout),dimension(nprof,ngate,rcfg%nhclass) :: &
101 : Np_matrix ! Table of hydrometeor number concentration. 0 ==> use defaults. (units = 1/kg)
102 :
103 : ! OUTPUTS
104 : real(wp),intent(out), dimension(nprof, ngate) :: &
105 : z_vol, & ! Effective reflectivity factor (mm^6/m^3)
106 : kr_vol ! Attenuation coefficient hydro (dB/km)
107 :
108 : ! INTERNAL VARIABLES
109 : integer :: &
110 : phase, ns,tp,j,k,pr,itt,iRe_type,n
111 : logical :: &
112 : hydro
113 : real(wp) :: &
114 : t_kelvin,Re_internal
115 : real(wp) :: &
116 : rho_a,kr,ze,zr,scale_factor,Re,Np,base,step
117 :
118 : real(wp),dimension(:),allocatable :: &
119 92880 : Deq, & ! Discrete drop sizes (um)
120 92880 : Ni, & ! Discrete concentrations (cm^-3 um^-1)
121 92880 : rhoi, & ! Discrete densities (kg m^-3)
122 92880 : xxa, & !
123 92880 : Di ! Discrete drop sizes (um)
124 :
125 : real(wp), dimension(nprof, ngate) :: &
126 92880 : z_ray ! Reflectivity factor, Rayleigh only (mm^6/m^3)
127 :
128 : ! PARAMETERS
129 : logical, parameter :: & !
130 : DO_LUT_TEST = .false., & !
131 : DO_NP_TEST = .false. !
132 : real(wp), parameter :: &
133 : one_third = 1._wp/3._wp !
134 :
135 : ! Initialization
136 130366800 : z_vol = 0._wp
137 130366800 : z_ray = 0._wp
138 130366800 : kr_vol = 0._wp
139 :
140 7894800 : do k=1,ngate ! Loop over each profile (nprof)
141 130366800 : do pr=1,nprof
142 :
143 : ! Determine if hydrometeor(s) present in volume
144 122472000 : hydro = .false.
145 951717079 : do j=1,rcfg%nhclass
146 951717079 : if ((hm_matrix(pr,k,j) > 1E-12) .and. (sd%dtype(j) > 0)) then
147 : hydro = .true.
148 : exit
149 : endif
150 : enddo
151 :
152 122472000 : t_kelvin = t_matrix(pr,k)
153 : ! If there is hydrometeor in the volume
154 130273920 : if (hydro) then
155 40279215 : rho_a = (p_matrix(pr,k))/(287._wp*(t_kelvin))
156 :
157 : ! Loop over hydrometeor type
158 402792150 : do tp=1,rcfg%nhclass
159 362512935 : Re_internal = re_matrix(pr,k,tp)
160 :
161 362512935 : if (hm_matrix(pr,k,tp) <= 1E-12) cycle
162 :
163 : ! Index into temperature dimension of scaling tables
164 : ! These tables have regular steps -- exploit this and abandon infind
165 57569436 : phase = sd%phase(tp)
166 57569436 : if (phase==0) then
167 26172468 : itt = infind(mt_ttl,t_kelvin)
168 : else
169 31396968 : itt = infind(mt_tti,t_kelvin)
170 : endif
171 :
172 : ! Compute effective radius from number concentration and distribution parameters
173 57569436 : if (Re_internal .eq. 0) then
174 : call calc_Re(hm_matrix(pr,k,tp),Np_matrix(pr,k,tp),rho_a, &
175 0 : sd%dtype(tp),sd%apm(tp),sd%bpm(tp),sd%rho(tp),sd%p1(tp),sd%p2(tp),sd%p3(tp),Re)
176 0 : Re_internal=Re
177 : !re_matrix(pr,k,tp)=Re
178 : else
179 57569436 : if (Np_matrix(pr,k,tp) > 0) then
180 : call errorMessage('WARNING(optics/quickbeam_optics.f90): '//&
181 0 : 'Re and Np set for the same volume & hydrometeor type. Np is being ignored.')
182 : endif
183 57569436 : Re = Re_internal
184 : !Re = re_matrix(pr,k,tp)
185 : endif
186 :
187 : ! Index into particle size dimension of scaling tables
188 57569436 : iRe_type=1
189 57569436 : if(Re.gt.0) then
190 : ! Determine index in to scale LUT
191 : ! Distance between Re points (defined by "base" and "step") for
192 : ! each interval of size Re_BIN_LENGTH
193 : ! Integer asignment, avoids calling floor intrinsic
194 57569436 : n=Re/Re_BIN_LENGTH
195 57569436 : if (n>=Re_MAX_BIN) n=Re_MAX_BIN-1
196 57569436 : step = rcfg%step_list(n+1)
197 57569436 : base = rcfg%base_list(n+1)
198 57569436 : iRe_type=Re/step
199 57569436 : if (iRe_type.lt.1) iRe_type=1
200 57569436 : Re=step*(iRe_type+0.5_wp) ! set value of Re to closest value allowed in LUT.
201 57569436 : iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)
202 :
203 : ! Make sure iRe_type is within bounds
204 57569436 : if (iRe_type.ge.nRe_types) then
205 : !write(*,*) 'Warning: size of Re exceed value permitted ', &
206 : ! 'in Look-Up Table (LUT). Will calculate. '
207 : ! No scaling allowed
208 25118624 : iRe_type=nRe_types
209 25118624 : rcfg%Z_scale_flag(tp,itt,iRe_type)=.false.
210 : else
211 : ! Set value in re_matrix to closest values in LUT
212 32450812 : if (.not. DO_LUT_TEST) re_internal=Re
213 : !if (.not. DO_LUT_TEST) re_matrix(pr,k,tp)=Re
214 : endif
215 : endif
216 :
217 : ! Use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
218 : ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
219 : ! if( rcfg%Z_scale_flag(tp,itt,iRe_type) .and. .not. DO_LUT_TEST) then
220 : ! ! can use z scaling
221 : ! scale_factor=rho_a*hm_matrix(pr,k,tp)
222 : ! zr = rcfg%Zr_scaled(tp,itt,iRe_type) * scale_factor
223 : ! ze = rcfg%Ze_scaled(tp,itt,iRe_type) * scale_factor
224 : ! kr = rcfg%kr_scaled(tp,itt,iRe_type) * scale_factor
225 : ! else
226 57569436 : if( (.not. rcfg%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST) then
227 : ! Create a discrete distribution of hydrometeors within volume
228 27950468 : select case(sd%dtype(tp))
229 : case(4)
230 0 : ns = 1
231 0 : allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
232 0 : Di = sd%p1(tp)
233 0 : Ni = 0._wp
234 : case default
235 27950468 : ns = nd ! constant defined in simulator/quickbeam.f90
236 27950468 : allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
237 2431690716 : Di = D
238 2431690716 : Ni = 0._wp
239 : end select
240 : call dsd(hm_matrix(pr,k,tp),re_internal,Np_matrix(pr,k,tp), &
241 : Di,Ni,ns,sd%dtype(tp),rho_a,t_kelvin, &
242 : sd%dmin(tp),sd%dmax(tp),sd%apm(tp),sd%bpm(tp), &
243 27950468 : sd%rho(tp),sd%p1(tp),sd%p2(tp),sd%p3(tp))
244 :
245 : ! Calculate particle density
246 27950468 : if (phase == 1) then
247 13633725 : if (sd%rho(tp) < 0) then
248 : ! Use equivalent volume spheres.
249 17877078 : rcfg%rho_eff(tp,1:ns,iRe_type) = rhoice ! solid ice == equivalent volume approach
250 18084951 : Deq = ( ( 6/pi*sd%apm(tp)/rhoice) ** one_third ) * ( (Di*1E-6) ** (sd%bpm(tp)/3._wp) ) * 1E6
251 : ! alternative is to comment out above two lines and use the following block
252 : ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density
253 : !
254 : ! rcfg%rho_eff(tp,1:ns,iRe_type) = (6/pi)*sd%apm(tp)*(Di*1E-6)**(sd%bpm(tp)-3) !MG Mie approach
255 :
256 : ! as the particle size gets small it is possible that the mass to size relationship of
257 : ! (given by power law in hclass.data) can produce impossible results
258 : ! where the mass is larger than a solid sphere of ice.
259 : ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
260 : ! do i=1,ns
261 : ! if(rcfg%rho_eff(tp,i,iRe_type) > 917 ) then
262 : ! rcfg%rho_eff(tp,i,iRe_type) = 917
263 : ! endif
264 : ! enddo
265 : else
266 : ! Equivalent volume sphere (solid ice rhoice=917 kg/m^3).
267 1154623272 : rcfg%rho_eff(tp,1:ns,iRe_type) = rhoice
268 1168049124 : Deq=Di * ((sd%rho(tp)/rhoice)**one_third)
269 : ! alternative ... coment out above two lines and use the following for MG-Mie
270 : ! rcfg%rho_eff(tp,1:ns,iRe_type) = sd%rho(tp) !MG Mie approach
271 : endif
272 : else
273 : ! I assume here that water phase droplets are spheres.
274 : ! sd%rho should be ~ 1000 or sd%apm=524 .and. sd%bpm=3
275 1245556641 : Deq = Di
276 : endif
277 :
278 : ! Calculate effective reflectivity factor of volume
279 : ! xxa are unused (Mie scattering and extinction efficiencies)
280 2403740248 : xxa(1:ns) = -9.9_wp
281 2431690716 : rhoi = rcfg%rho_eff(tp,1:ns,iRe_type)
282 : call zeff(rcfg%freq,Deq,Ni,ns,rcfg%k2,t_kelvin,phase,rcfg%do_ray, &
283 27950468 : ze,zr,kr,xxa,xxa,rhoi)
284 :
285 : ! Test compares total number concentration with sum of discrete samples
286 : ! The second test, below, compares ab initio and "scaled" computations
287 : ! of reflectivity
288 : ! These should get broken out as a unit test that gets called on
289 : ! data. That routine could write to std out.
290 :
291 : ! Test code ... compare Np value input to routine with sum of DSD
292 : ! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation
293 : ! not just the DSD representation given by Ni
294 : if(Np_matrix(pr,k,tp)>0 .and. DO_NP_TEST ) then
295 : Np = path_integral(Ni,Di,1,ns-1)/rho_a*1.E6_wp
296 : ! Note: Representation is not great or small Re < 2
297 : if( (Np_matrix(pr,k,tp)-Np)/Np_matrix(pr,k,tp)>0.1 ) then
298 : call errorMessage('ERROR(optics/quickbeam_optics.f90): Error: Np input does not match sum(N)')
299 : endif
300 : endif
301 :
302 : ! Clean up space
303 27950468 : deallocate(Di,Ni,rhoi,xxa,Deq)
304 :
305 : ! LUT test code
306 : ! This segment of code compares full calculation to scaling result
307 : if ( rcfg%Z_scale_flag(tp,itt,iRe_type) .and. DO_LUT_TEST ) then
308 : scale_factor=rho_a*hm_matrix(pr,k,tp)
309 : ! if more than 2 dBZe difference print error message/parameters.
310 : if ( abs(10*log10(ze) - 10*log10(rcfg%Ze_scaled(tp,itt,iRe_type) * &
311 : scale_factor)) > 2 ) then
312 : call errorMessage('ERROR(optics/quickbeam_optics.f90): ERROR: Roj Error?')
313 : endif
314 : endif
315 : else
316 : ! Use z scaling
317 29618968 : scale_factor=rho_a*hm_matrix(pr,k,tp)
318 29618968 : zr = rcfg%Zr_scaled(tp,itt,iRe_type) * scale_factor
319 29618968 : ze = rcfg%Ze_scaled(tp,itt,iRe_type) * scale_factor
320 29618968 : kr = rcfg%kr_scaled(tp,itt,iRe_type) * scale_factor
321 : endif ! end z_scaling
322 :
323 57569436 : kr_vol(pr,k) = kr_vol(pr,k) + kr
324 57569436 : z_vol(pr,k) = z_vol(pr,k) + ze
325 57569436 : z_ray(pr,k) = z_ray(pr,k) + zr
326 :
327 : ! Construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
328 97848651 : if ( .not. rcfg%Z_scale_flag(tp,itt,iRe_type) ) then
329 27950468 : if (iRe_type>1) then
330 27950468 : scale_factor=rho_a*hm_matrix(pr,k,tp)
331 27950468 : rcfg%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
332 27950468 : rcfg%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
333 27950468 : rcfg%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
334 27950468 : rcfg%Z_scale_flag(tp,itt,iRe_type) = .true.
335 27950468 : rcfg%Z_scale_added_flag(tp,itt,iRe_type)=.true.
336 : endif
337 : endif
338 : enddo ! end loop of tp (hydrometeor type)
339 : endif
340 : enddo
341 : enddo
342 :
343 130366800 : where(kr_vol(:,:) <= EPSILON(kr_vol))
344 : ! Volume is hydrometeor-free
345 : !z_vol(:,:) = undef
346 : z_ray(:,:) = undef
347 : end where
348 :
349 : ! Save any updates made
350 92880 : if (rcfg%update_scale_LUTs) call save_scale_LUTs(rcfg)
351 :
352 92880 : end subroutine quickbeam_optics
353 : ! ##############################################################################################
354 : ! ##############################################################################################
355 0 : subroutine calc_Re(Q,Np,rho_a,dtype,apm,bpm,rho_c,p1,p2,p3,Re)
356 : ! ##############################################################################################
357 : ! Purpose:
358 : ! Calculates Effective Radius (1/2 distribution 3rd moment / 2nd moment).
359 : !
360 : ! For some distribution types, the total number concentration (per kg), Np
361 : ! may be optionally specified. Should be set to zero, otherwise.
362 : !
363 : ! Roj Marchand July 2010
364 : !
365 : ! Inputs:
366 : !
367 : ! [Q] hydrometeor mixing ratio (g/kg) ! not needed for some distribution types
368 : ! [Np] Optional Total number concentration (per kg). 0 = use defaults (p1, p2, p3)
369 : ! [rho_a] ambient air density (kg m^-3)
370 : !
371 : ! Distribution parameters as per quickbeam documentation.
372 : ! [dtype] distribution type
373 : ! [apm] a parameter for mass (kg m^[-bpm])
374 : ! [bmp] b params for mass
375 : ! [p1],[p2],[p3] distribution parameters
376 : !
377 : ! Outputs:
378 : ! [Re] Effective radius, 1/2 the 3rd moment/2nd moment (um)
379 : !
380 : ! Created:
381 : ! July 2010 Roj Marchand
382 : ! Modified:
383 : ! 12/18/14 Dustin Swales: Define type REALs as double precision (dustin.swales@noaa.gov)
384 : !
385 : ! ##############################################################################################
386 : ! ##############################################################################################
387 :
388 : ! Inputs
389 : real(wp), intent(in) :: Q,Np,rho_a,rho_c,p1,p2,p3
390 : integer, intent(in) :: dtype
391 : real(wp), intent(inout) :: apm,bpm
392 :
393 : ! Outputs
394 : real(wp), intent(out) :: Re
395 :
396 : ! Internal
397 : integer :: local_dtype
398 : real(wp) :: local_p3,local_Np,tmp1,tmp2
399 : real(wp) :: N0,D0,vu,dm,ld,rg,log_sigma_g ! gamma, exponential variables
400 :
401 :
402 : ! If density is constant, set equivalent values for apm and bpm
403 0 : if ((rho_c > 0) .and. (apm < 0)) then
404 0 : apm = (pi/6)*rho_c
405 0 : bpm = 3._wp
406 : endif
407 :
408 : ! Exponential is same as modified gamma with vu =1
409 : ! if Np is specified then we will just treat as modified gamma
410 0 : if(dtype .eq. 2 .and. Np .gt. 0) then
411 : local_dtype = 1
412 : local_p3 = 1
413 : else
414 0 : local_dtype = dtype
415 0 : local_p3 = p3
416 : endif
417 0 : select case(local_dtype)
418 :
419 : ! ---------------------------------------------------------!
420 : ! Modified gamma !
421 : ! Np = total number concentration (1/kg) = Nt / rho_a !
422 : ! D0 = characteristic diameter (um) !
423 : ! dm = mean diameter (um) - first moment over zeroth moment!
424 : ! vu = distribution width parameter !
425 : ! ---------------------------------------------------------!
426 : case(1)
427 :
428 0 : if( abs(local_p3+2) < 1E-8) then
429 0 : if(Np>1E-30) then
430 : ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1)
431 : ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane
432 0 : vu = (1/(0.2714_wp + 0.00057145_wp*Np*rho_a*1E-6))**2 ! units of Nt = Np*rhoa = #/cm^3
433 : else
434 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
435 0 : 'Must specify a value for Np in each volume with Morrison/Martin Scheme.')
436 0 : return
437 : endif
438 0 : elseif (abs(local_p3+1) > 1E-8) then
439 : ! vu is fixed in hp structure
440 0 : vu = local_p3
441 : else
442 : ! vu isn't specified
443 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
444 0 : 'Must specify a value for vu for Modified Gamma distribution')
445 0 : return
446 : endif
447 :
448 0 : if( Np.eq.0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default
449 0 : dm = p2 ! by definition, should have units of microns
450 0 : D0 = gamma(vu)/gamma(vu+1)*dm
451 : else ! use value of Np
452 0 : if(Np.eq.0) then
453 0 : if( abs(p1+1) > 1E-8 ) then ! use default number concentration
454 : local_Np = p1 ! total number concentration / pa --- units kg^-1
455 : else
456 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
457 0 : 'Must specify Np or default value (p1=Dm [um] or p2=Np [1/kg]) for Modified Gamma distribution')
458 0 : return
459 : endif
460 : else
461 : local_Np=Np;
462 : endif
463 0 : D0 = 1E6 * ( Q*1E-3*gamma(vu)/(apm*local_Np*gamma(vu+bpm)) )**(1/bpm) ! units = microns
464 : endif
465 0 : Re = 0.5_wp*D0*gamma(vu+3)/gamma(vu+2)
466 :
467 : ! ---------------------------------------------------------!
468 : ! Exponential !
469 : ! N0 = intercept parameter (m^-4) !
470 : ! ld = slope parameter (um) !
471 : ! ---------------------------------------------------------!
472 : case(2)
473 :
474 : ! Np not specified (see if statement above)
475 0 : if((abs(p1+1) > 1E-8) ) then ! N0 has been specified, determine ld
476 0 : N0 = p1
477 0 : tmp1 = 1._wp/(1._wp+bpm)
478 0 : ld = ((apm*gamma(1.+bpm)*N0)/(rho_a*Q*1E-3))**tmp1
479 0 : ld = ld/1E6 ! set units to microns^-1
480 0 : elseif (abs(p2+1) > 1E-8) then ! lambda=ld has been specified as default
481 : ld = p2 ! should have units of microns^-1
482 : else
483 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
484 0 : 'Must specify Np or default value (p1=No or p2=lambda) for Exponential distribution')
485 0 : return
486 : endif
487 0 : Re = 1.5_wp/ld
488 :
489 : ! ---------------------------------------------------------!
490 : ! Power law !
491 : ! ahp = Ar parameter (m^-4 mm^-bhp) !
492 : ! bhp = br parameter !
493 : ! dmin_mm = lower bound (mm) !
494 : ! dmax_mm = upper bound (mm) !
495 : ! ---------------------------------------------------------!
496 : case(3)
497 :
498 0 : Re=0._wp ! Not supporting LUT approach for power-law ...
499 0 : if(Np>0) then
500 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
501 0 : 'Variable Np not supported for Power Law distribution')
502 0 : return
503 : endif
504 :
505 : ! ---------------------------------------------------------!
506 : ! Monodisperse !
507 : ! D0 = particle diameter (um) == Re !
508 : ! ---------------------------------------------------------!
509 : case(4)
510 :
511 0 : Re = p1
512 0 : if(Np>0) then
513 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
514 0 : 'Variable Np not supported for Monodispersed distribution')
515 0 : return
516 : endif
517 :
518 : ! ---------------------------------------------------------!
519 : ! Lognormal !
520 : ! N0 = total number concentration (m^-3) !
521 : ! np = fixed number concentration (kg^-1) !
522 : ! rg = mean radius (um) !
523 : ! log_sigma_g = ln(geometric standard deviation) !
524 : ! ---------------------------------------------------------!
525 : case(5)
526 :
527 0 : if( abs(local_p3+1) > 1E-8 ) then
528 : !set natural log width
529 0 : log_sigma_g = local_p3
530 : else
531 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
532 0 : 'Must specify a value for sigma_g when using a Log-Normal distribution')
533 0 : return
534 : endif
535 :
536 : ! get rg ...
537 0 : if( Np.eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
538 : rg = p2
539 : else
540 0 : if(Np>0) then
541 : local_Np=Np;
542 0 : elseif(abs(p2+1) < 1E-8) then
543 0 : local_Np=p1
544 : else
545 : call errorMessage('ERROR(optics/quickbeam_optics.f90:Calc_Re): '//&
546 0 : 'Must specify Np or default value (p2=Rg or p1=Np) for Log-Normal distribution')
547 : endif
548 0 : log_sigma_g = p3
549 0 : tmp1 = (Q*1E-3)/(2._wp**bpm*apm*local_Np)
550 0 : tmp2 = exp(0.5_wp*bpm*bpm*(log_sigma_g))*exp(0.5_wp*bpm*bpm*(log_sigma_g))
551 0 : rg = ((tmp1/tmp2)**(1._wp/bpm))*1E6
552 : endif
553 0 : Re = rg*exp(2.5_wp*(log_sigma_g*log_sigma_g))
554 : end select
555 : end subroutine calc_Re
556 : ! ##############################################################################################
557 : ! ##############################################################################################
558 27950468 : subroutine dsd(Q,Re,Np,D,N,nsizes,dtype,rho_a,tk,dmin,dmax,apm,bpm,rho_c,p1,p2,p3)
559 : ! ##############################################################################################
560 : ! Purpose:
561 : ! Create a discrete drop size distribution
562 : !
563 : ! Starting with Quickbeam V3, this routine now allows input of
564 : ! both effective radius (Re) and total number concentration (Nt)
565 : ! Roj Marchand July 2010
566 : !
567 : ! The version in Quickbeam v.104 was modified to allow Re but not Nt
568 : ! This is a significantly modified form for the version
569 : !
570 : ! Originally Part of QuickBeam v1.03 by John Haynes
571 : ! http://reef.atmos.colostate.edu/haynes/radarsim
572 : !
573 : ! Inputs:
574 : !
575 : ! [Q] hydrometeor mixing ratio (g/kg)
576 : ! [Re] Optional Effective Radius (microns). 0 = use defaults (p1, p2, p3)
577 : !
578 : ! [D] array of discrete drop sizes (um) where we desire to know the number concentraiton n(D).
579 : ! [nsizes] number of elements of [D]
580 : !
581 : ! [dtype] distribution type
582 : ! [rho_a] ambient air density (kg m^-3)
583 : ! [tk] temperature (K)
584 : ! [dmin] minimum size cutoff (um)
585 : ! [dmax] maximum size cutoff (um)
586 : ! [rho_c] alternate constant density (kg m^-3)
587 : ! [p1],[p2],[p3] distribution parameters
588 : !
589 : ! Input/Output:
590 : ! [apm] a parameter for mass (kg m^[-bpm])
591 : ! [bmp] b params for mass
592 : !
593 : ! Outputs:
594 : ! [N] discrete concentrations (cm^-3 um^-1)
595 : ! or, for monodisperse, a constant (1/cm^3)
596 : !
597 : ! Requires:
598 : ! function infind
599 : !
600 : ! Created:
601 : ! 11/28/05 John Haynes (haynes@atmos.colostate.edu)
602 : ! Modified:
603 : ! 01/31/06 Port from IDL to Fortran 90
604 : ! 07/07/06 Rewritten for variable DSD's
605 : ! 10/02/06 Rewritten using scaling factors (Roger Marchand and JMH), Re added V1.04
606 : ! July 2020 "N Scale factors" (variable fc) removed (Roj Marchand).
607 : ! 12/18/14 Define type REALs as double precision (dustin.swales@noaa.gov)
608 : ! ##############################################################################################
609 :
610 : ! Inputs
611 : integer, intent(in) :: &
612 : nsizes,& ! Number of elements of [D]
613 : dtype ! distribution type
614 : real(wp),intent(in),dimension(nsizes) :: &
615 : D ! Array of discrete drop sizes (um) where we desire to know the number concentraiton n(D).
616 : real(wp),intent(in) :: &
617 : Q, & ! Hydrometeor mixing ratio (g/kg)
618 : Np, & !
619 : rho_a, & ! Ambient air density (kg m^-3)
620 : tk, & ! Temperature (K)
621 : dmin, & ! Minimum size cutoff (um)
622 : dmax, & ! Maximum size cutoff (um)
623 : rho_c, & ! Alternate constant density (kg m^-3)
624 : p1, & ! Distribution parameter 1
625 : p2, & ! Distribution parameter 2
626 : p3 ! Distribution parameter 3
627 : real(wp),intent(inout) :: &
628 : apm, & ! a parameter for mass (kg m^[-bpm])
629 : bpm, & ! b params for mass
630 : Re ! Optional Effective Radius (microns)
631 :
632 : ! Outputs
633 : real(wp),intent(out),dimension(nsizes) :: &
634 : N ! Discrete concentrations (cm^-3 um^-1)
635 : ! or, for monodisperse, a constant (1/cm^3)
636 :
637 : ! Internal Variables
638 : real(wp),dimension(nsizes) :: &
639 27950468 : fc
640 : real(wp) :: &
641 : N0,D0,vu,local_np,dm,ld, & ! gamma, exponential variables
642 : dmin_mm,dmax_mm,ahp,bhp, & ! power law variables
643 : rg,log_sigma_g, & ! lognormal variables
644 : rho_e, & ! particle density (kg m^-3)
645 : tmp1,tmp2,tc
646 : integer :: &
647 : k,lidx,uidx
648 :
649 : ! Convert temperature from Kelvin to Celsius
650 27950468 : tc = tk - 273.15_wp
651 :
652 : ! If density is constant, store equivalent values for apm and bpm
653 27950468 : if ((rho_c > 0) .and. (apm < 0)) then
654 7292 : apm = (pi/6)*rho_c
655 7292 : bpm = 3._wp
656 : endif
657 :
658 : ! Will preferentially use Re input over Np.
659 : ! if only Np given then calculate Re
660 : ! if neigher than use other defaults (p1,p2,p3) following quickbeam documentation
661 27950468 : if(Re==0 .and. Np>0) then
662 0 : call calc_Re(Q,Np,rho_a,dtype,apm,bpm,rho_c,p1,p2,p3,Re)
663 : endif
664 28158341 : select case(dtype)
665 :
666 : ! ---------------------------------------------------------!
667 : ! Modified gamma !
668 : ! np = total number concentration !
669 : ! D0 = characteristic diameter (um) !
670 : ! dm = mean diameter (um) - first moment over zeroth moment!
671 : ! vu = distribution width parameter !
672 : ! ---------------------------------------------------------!
673 : case(1)
674 :
675 207873 : if( abs(p3+2) < 1E-8) then
676 0 : if( Np>1E-30) then
677 : ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1)
678 : ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane
679 0 : vu = (1/(0.2714_wp + 0.00057145_wp*Np*rho_a*1E-6))**2._wp ! units of Nt = Np*rhoa = #/cm^3
680 : else
681 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
682 0 : 'Must specify a value for Np in each volume with Morrison/Martin Scheme.')
683 0 : return
684 : endif
685 207873 : elseif (abs(p3+1) > 1E-8) then
686 : ! vu is fixed in hp structure
687 207873 : vu = p3
688 : else
689 : ! vu isn't specified
690 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
691 0 : 'Must specify a value for vu for Modified Gamma distribution')
692 0 : return
693 : endif
694 :
695 207873 : if(Re>0) then
696 207873 : D0 = 2._wp*Re*gamma(vu+2)/gamma(vu+3)
697 : fc = (((D*1E-6)**(vu-1)*exp(-1*D/D0)) / &
698 17877078 : (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm))) * 1E-12
699 17877078 : N = fc*rho_a*(Q*1E-3)
700 0 : elseif( p2+1 > 1E-8) then ! use default value for MEAN diameter
701 0 : dm = p2
702 0 : D0 = gamma(vu)/gamma(vu+1)*dm
703 : fc = (((D*1E-6)**(vu-1)*exp(-1*D/D0)) / &
704 0 : (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm))) * 1E-12
705 0 : N = fc*rho_a*(Q*1E-3)
706 0 : elseif(abs(p3+1) > 1E-8) then! use default number concentration
707 0 : local_np = p1 ! total number concentration / pa check
708 0 : tmp1 = (Q*1E-3)**(1./bpm)
709 0 : fc = (D*1E-6 / (gamma(vu)/(apm*local_np*gamma(vu+bpm)))**(1._wp/bpm))**vu
710 : N = ((rho_a*local_np*fc*(D*1E-6)**(-1._wp))/(gamma(vu)*tmp1**vu) * &
711 0 : exp(-1._wp*fc**(1._wp/vu)/tmp1)) * 1E-12
712 : else
713 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
714 0 : 'No default value for Dm or Np provided!')
715 0 : return
716 : endif
717 :
718 : ! ---------------------------------------------------------!
719 : ! Exponential !
720 : ! N0 = intercept parameter (m^-4) !
721 : ! ld = slope parameter (um) !
722 : ! ---------------------------------------------------------!
723 : case(2)
724 :
725 27517235 : if(Re>0) then
726 27517235 : ld = 1.5_wp/Re ! units 1/um
727 2366482210 : fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))*exp(-1._wp*(ld*1E6)*(D*1E-6))*1E-12
728 2366482210 : N = fc*rho_a*(Q*1E-3)
729 0 : elseif (abs(p1+1) > 1E-8) then
730 : ! Use N0 default value
731 0 : N0 = p1
732 0 : tmp1 = 1._wp/(1._wp+bpm)
733 0 : fc = ((apm*gamma(1.+bpm)*N0)**tmp1)*(D*1E-6)
734 0 : N = (N0*exp(-1._wp*fc*(1._wp/(rho_a*Q*1E-3))**tmp1)) * 1E-12
735 0 : elseif (abs(p2+1) > 1E-8) then
736 : ! Use default value for lambda
737 0 : ld = p2
738 0 : fc = (ld*1E6)**(1._wp+bpm)/(apm*gamma(1+bpm))*exp(-1._wp*(ld*1E6)*(D*1E-6))*1E-12
739 0 : N = fc*rho_a*(Q*1E-3)
740 : else
741 : ! ld "parameterized" from temperature (carry over from original Quickbeam).
742 0 : ld = 1220._wp*10._wp**(-0.0245_wp*tc)*1E-6
743 0 : N0 = ((ld*1E6)**(1._wp+bpm)*Q*1E-3*rho_a)/(apm*gamma(1+bpm))
744 0 : N = (N0*exp(-ld*D)) * 1E-12
745 : endif
746 :
747 : ! ---------------------------------------------------------!
748 : ! Power law !
749 : ! ahp = Ar parameter (m^-4 mm^-bhp) !
750 : ! bhp = br parameter !
751 : ! dmin_mm = lower bound (mm) !
752 : ! dmax_mm = upper bound (mm) !
753 : ! ---------------------------------------------------------!
754 : case(3)
755 :
756 0 : if(Re>0) then
757 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
758 0 : 'Variable Re not supported for Power-Law distribution')
759 0 : return
760 0 : elseif(Np>0) then
761 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
762 0 : 'Variable Np not supported for Power-Law distribution')
763 0 : return
764 : endif
765 :
766 : ! br parameter
767 0 : if (abs(p1+2) < 1E-8) then
768 : ! if p1=-2, bhp is parameterized according to Ryan (2000),
769 : ! applicatable to cirrus clouds
770 0 : if (tc < -30) then
771 0 : bhp = -1.75_wp+0.09_wp*((tc+273._wp)-243.16_wp)
772 0 : elseif ((tc >= -30) .and. (tc < -9)) then
773 0 : bhp = -3.25_wp-0.06_wp*((tc+273._wp)-265.66_wp)
774 : else
775 : bhp = -2.15_wp
776 : endif
777 0 : elseif (abs(p1+3) < 1E-8) then
778 : ! if p1=-3, bhp is parameterized according to Ryan (2000),
779 : ! applicable to frontal clouds
780 0 : if (tc < -35) then
781 0 : bhp = -1.75_wp+0.09_wp*((tc+273._wp)-243.16_wp)
782 0 : elseif ((tc >= -35) .and. (tc < -17.5)) then
783 0 : bhp = -2.65_wp+0.09_wp*((tc+273._wp)-255.66_wp)
784 0 : elseif ((tc >= -17.5) .and. (tc < -9)) then
785 0 : bhp = -3.25_wp-0.06_wp*((tc+273._wp)-265.66_wp)
786 : else
787 : bhp = -2.15_wp
788 : endif
789 : else
790 : ! Otherwise the specified value is used
791 : bhp = p1
792 : endif
793 :
794 : ! Ar parameter
795 0 : dmin_mm = dmin*1E-3
796 0 : dmax_mm = dmax*1E-3
797 :
798 : ! Commented lines are original method with constant density
799 : ! rc = 500. ! (kg/m^3)
800 : ! tmp1 = 6*rho_a*(bhp+4)
801 : ! tmp2 = pi*rc*(dmax_mm**(bhp+4))*(1-(dmin_mm/dmax_mm)**(bhp+4))
802 : ! ahp = (Q*1E-3)*1E12*tmp1/tmp2
803 :
804 : ! New method is more consistent with the rest of the distributions
805 : ! and allows density to vary with particle size
806 0 : tmp1 = rho_a*(Q*1E-3)*(bhp+bpm+1)
807 0 : tmp2 = apm*(dmax_mm**bhp*dmax**(bpm+1)-dmin_mm**bhp*dmin**(bpm+1))
808 0 : ahp = tmp1/tmp2 * 1E24
809 : ! ahp = tmp1/tmp2
810 0 : lidx = infind(D,dmin)
811 0 : uidx = infind(D,dmax)
812 0 : do k=lidx,uidx
813 0 : N(k) = (ahp*(D(k)*1E-3)**bhp) * 1E-12
814 : enddo
815 :
816 : ! ---------------------------------------------------------!
817 : ! Monodisperse !
818 : ! D0 = particle diameter (um) !
819 : ! ---------------------------------------------------------!
820 : case(4)
821 :
822 0 : if (Re>0) then
823 : D0 = Re
824 : else
825 0 : D0 = p1
826 : endif
827 :
828 0 : rho_e = (6._wp/pi)*apm*(D0*1E-6)**(bpm-3)
829 0 : fc(1) = (6._wp/(pi*D0*D0*D0*rho_e))*1E12
830 0 : N(1) = fc(1)*rho_a*(Q*1E-3)
831 :
832 : ! ---------------------------------------------------------!
833 : ! Lognormal !
834 : ! N0 = total number concentration (m^-3) !
835 : ! np = fixed number concentration (kg^-1) !
836 : ! rg = mean radius (um) !
837 : ! og_sigma_g = ln(geometric standard deviation) !
838 : ! ---------------------------------------------------------!
839 : case(5)
840 27950468 : if (abs(p1+1) < 1E-8 .or. Re>0 ) then
841 : ! rg, log_sigma_g are given
842 225360 : log_sigma_g = p3
843 225360 : tmp2 = (bpm*log_sigma_g)*(bpm*log_sigma_g)
844 225360 : if(Re.le.0) then
845 0 : rg = p2
846 : else
847 : !rg = Re*exp(-2.5*(log_sigma_g*log_sigma_g))
848 225360 : rg =Re*exp(-2.5_wp*(log_sigma_g**2))
849 :
850 : endif
851 :
852 : fc = 0.5_wp*((1._wp/((2._wp*rg*1E-6)**(bpm)*apm*(2._wp*pi)**(0.5_wp) * &
853 19380960 : log_sigma_g*D*0.5_wp*1E-6))*exp(-0.5_wp*((log(0.5_wp*D/rg)/log_sigma_g)**2._wp+tmp2)))*1E-12
854 19380960 : N = fc*rho_a*(Q*1E-3)
855 :
856 0 : elseif (abs(p2+1) < 1E-8 .or. Np>0) then
857 : ! Np, log_sigma_g are given
858 0 : if(Np>0) then
859 : local_Np = Np
860 : else
861 0 : local_Np = p1
862 : endif
863 :
864 0 : log_sigma_g = p3
865 0 : N0 = local_np*rho_a
866 0 : tmp1 = (rho_a*(Q*1E-3))/(2._wp**bpm*apm*N0)
867 0 : tmp2 = exp(0.5_wp*bpm*bpm*(log_sigma_g))*exp(0.5_wp*bpm*bpm*(log_sigma_g))
868 0 : rg = ((tmp1/tmp2)**(1/bpm))*1E6
869 :
870 : N = 0.5_wp*(N0 / ((2._wp*pi)**(0.5_wp)*log_sigma_g*D*0.5_wp*1E-6) * &
871 0 : exp((-0.5_wp*(log(0.5_wp*D/rg)/log_sigma_g)**2._wp)))*1E-12
872 : else
873 : call errorMessage('FATAL ERROR(optics/quickbeam_optics.f90:dsd): '//&
874 0 : 'Must specify a value for sigma_g')
875 0 : return
876 : endif
877 : end select
878 : end subroutine dsd
879 : ! ##############################################################################################
880 : ! ##############################################################################################
881 27950468 : subroutine zeff(freq,D,N,nsizes,k2,tt,ice,xr,z_eff,z_ray,kr,qe,qs,rho_e)
882 : ! ##############################################################################################
883 : ! Purpose:
884 : ! Simulates radar return of a volume given DSD of spheres
885 : ! Part of QuickBeam v1.03 by John Haynes
886 : ! http://reef.atmos.colostate.edu/haynes/radarsim
887 : !
888 : ! Inputs:
889 : ! [freq] radar frequency (GHz)
890 : ! [D] discrete drop sizes (um)
891 : ! [N] discrete concentrations (cm^-3 um^-1)
892 : ! [nsizes] number of discrete drop sizes
893 : ! [k2] |K|^2, -1=use frequency dependent default
894 : ! [tt] hydrometeor temperature (K)
895 : ! [ice] indicates volume consists of ice
896 : ! [xr] perform Rayleigh calculations?
897 : ! [qe] if using a mie table, these contain ext/sca ...
898 : ! [qs] ... efficiencies; otherwise set to -1
899 : ! [rho_e] medium effective density (kg m^-3) (-1 = pure)
900 : !
901 : ! Outputs:
902 : ! [z_eff] unattenuated effective reflectivity factor (mm^6/m^3)
903 : ! [z_ray] reflectivity factor, Rayleigh only (mm^6/m^3)
904 : ! [kr] attenuation coefficient (db km^-1)
905 : !
906 : ! Created:
907 : ! 11/28/05 John Haynes (haynes@atmos.colostate.edu)
908 : ! Modified:
909 : ! 12/18/14 Dustin Swales: Define type REALs as double precision (dustin.swales@noaa.gov)
910 : ! ##############################################################################################
911 : ! Inputs
912 : integer, intent(in) :: &
913 : ice, & ! Indicates volume consists of ice
914 : xr, & ! Perform Rayleigh calculations?
915 : nsizes ! Number of discrete drop sizes
916 : real(wp), intent(in),dimension(nsizes) :: &
917 : D, & ! Discrete drop sizes (um)
918 : N, & ! Discrete concentrations (cm^-3 um^-1)
919 : rho_e, & ! Medium effective density (kg m^-3) (-1 = pure)
920 : qe, & ! Extinction efficiency, when using Mie tables
921 : qs ! Scatering efficiency, when using Mie tables
922 : real(wp),intent(in) :: &
923 : freq, & ! Radar frequency (GHz)
924 : tt ! Hydrometeor temperature (K)
925 : real(wp), intent(inout) :: &
926 : k2 ! |K|^2, -1=use frequency dependent default
927 :
928 : ! Outputs
929 : real(wp), intent(out) :: &
930 : z_eff, & ! Unattenuated effective reflectivity factor (mm^6/m^3)
931 : z_ray, & ! Reflectivity factor, Rayleigh only (mm^6/m^3)
932 : kr ! Attenuation coefficient (db km^-1)
933 :
934 : ! Internal Variables
935 : integer :: correct_for_rho ! Correct for density flag
936 : real(wp), dimension(nsizes) :: &
937 55900936 : D0, & ! D in (m)
938 55900936 : N0, & ! N in m^-3 m^-1
939 55900936 : sizep, & ! Size parameter
940 55900936 : qext, & ! Extinction efficiency
941 55900936 : qbsca, & ! Backscatter efficiency
942 55900936 : f, & ! Ice fraction
943 55900936 : xtemp !
944 : real(wp) :: &
945 : wl, cr,eta_sum,eta_mie,const,z0_eff,z0_ray,k_sum,n_r,n_i,dqv(1),dqsc,dg,dph(1)
946 : complex(wp) :: &
947 : m, & ! Complex index of refraction of bulk form
948 : Xs1(1), Xs2(1) !
949 : integer :: &
950 : i, err !
951 : integer, parameter :: &
952 : one=1 !
953 : real(wp),parameter :: &
954 : conv_d = 1e-6, & ! Conversion factor for drop sizes (to m)
955 : conv_n = 1e12, & ! Conversion factor for drop concentrations (to m^-3)
956 : conv_f = 0.299792458 ! Conversion for radar frequency (to m)
957 : complex(wp),dimension(nsizes) ::&
958 55900936 : m0 ! Complex index of refraction
959 :
960 : ! Initialize
961 27950468 : z0_ray = 0._wp
962 :
963 : ! Conversions
964 2403740248 : D0 = d*conv_d
965 2403740248 : N0 = n*conv_n
966 27950468 : wl = conv_f/freq
967 :
968 : ! // dielectric constant |k^2| defaults
969 27950468 : if (k2 < 0) then
970 9143 : k2 = 0.933_wp
971 9143 : if (abs(94.-freq) < 3.) k2=0.75_wp
972 9143 : if (abs(35.-freq) < 3.) k2=0.88_wp
973 9143 : if (abs(13.8-freq) < 3.) k2=0.925_wp
974 : endif
975 :
976 27950468 : if (qe(1) < -9) then
977 :
978 : ! Get the refractive index of the bulk hydrometeors
979 27950468 : if (ice == 0) then
980 14316743 : call m_wat(freq,tt,n_r,n_i)
981 : else
982 13633725 : call m_ice(freq,tt,n_r,n_i)
983 : endif
984 27950468 : m = cmplx(n_r,-n_i)
985 2403740248 : m0(1:nsizes) = m
986 :
987 : correct_for_rho = 0
988 2431690716 : if ((ice == 1) .and. (minval(rho_e) >= 0)) correct_for_rho = 1
989 :
990 : ! Correct refractive index for ice density if needed
991 : if (correct_for_rho == 1) then
992 1172500350 : f = rho_e/rhoice
993 1172500350 : m0 = sqrt((2+(m0*m0)+2*f*((m0*m0)-1))/(2+(m0*m0)+f*(1-(m0*m0))))
994 : endif
995 :
996 : ! Mie calculations
997 2403740248 : sizep = (pi*D0)/wl
998 27950468 : dqv(1) = 0._wp
999 2403740248 : do i=1,nsizes
1000 2375789780 : call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), &
1001 4779530028 : dg, xs1, xs2, dph, err)
1002 : end do
1003 :
1004 : else
1005 : ! Mie table used
1006 0 : qext = qe
1007 0 : qbsca = qs
1008 : endif
1009 :
1010 : ! eta_mie = 0.25*sum[qbsca*pi*D^2*N(D)*deltaD]
1011 : ! <--------- eta_sum --------->
1012 : ! z0_eff = (wl^4/!pi^5)*(1./k2)*eta_mie
1013 27950468 : eta_sum = 0._wp
1014 27950468 : if (size(D0) == 1) then
1015 0 : eta_sum = qbsca(1)*(n(1)*1E6)*D0(1)*D0(1)
1016 : else
1017 2403740248 : xtemp = qbsca*N0*D0*D0
1018 27950468 : call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),eta_sum)
1019 : endif
1020 :
1021 27950468 : eta_mie = eta_sum*0.25_wp*pi
1022 27950468 : const = ((wl*wl*wl*wl)/(pi*pi*pi*pi*pi))*(1._wp/k2)
1023 :
1024 27950468 : z0_eff = const*eta_mie
1025 :
1026 : ! kr = 0.25*cr*sum[qext*pi*D^2*N(D)*deltaD]
1027 : ! <---------- k_sum --------->
1028 27950468 : k_sum = 0._wp
1029 27950468 : if (size(D0) == 1) then
1030 0 : k_sum = qext(1)*(n(1)*1E6)*D0(1)*D0(1)
1031 : else
1032 2403740248 : xtemp = qext*N0*D0*D0
1033 27950468 : call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),k_sum)
1034 : endif
1035 : ! DS2014 START: Making this calculation in double precision results in a small
1036 : ! amount of very small errors in the COSP output field,dBZE94,
1037 : ! so it will be left as is.
1038 : !cr = 10._wp/log(10._wp)
1039 27950468 : cr = 10./log(10.)
1040 : ! DS2014 STOP
1041 27950468 : kr = k_sum*0.25_wp*pi*(1000._wp*cr)
1042 :
1043 : ! z_ray = sum[D^6*N(D)*deltaD]
1044 27950468 : if (xr == 1) then
1045 0 : z0_ray = 0._wp
1046 0 : if (size(D0) == 1) then
1047 0 : z0_ray = (n(1)*1E6)*D0(1)*D0(1)*D0(1)*D0(1)*D0(1)*D0(1)
1048 : else
1049 0 : xtemp = N0*D0*D0*D0*D0*D0*D0
1050 0 : call avint(xtemp,D0,nsizes,D0(1),D0(size(D0)),z0_ray)
1051 : endif
1052 : endif
1053 :
1054 : ! Convert to mm^6/m^3
1055 27950468 : z_eff = z0_eff*1E18 ! 10.*alog10(z0_eff*1E18)
1056 27950468 : z_ray = z0_ray*1E18 ! 10.*alog10(z0_ray*1E18)
1057 :
1058 27950468 : end subroutine zeff
1059 : ! ##############################################################################################
1060 : ! ##############################################################################################
1061 12247200 : function gases(PRES_mb,T,SH,f)
1062 : ! ##############################################################################################
1063 : ! Purpose:
1064 : ! Compute 2-way gaseous attenuation through a volume in microwave
1065 : !
1066 : ! Inputs:
1067 : ! [PRES_mb] pressure (mb) (hPa)
1068 : ! [T] temperature (K)
1069 : ! [RH] relative humidity (%)
1070 : ! [f] frequency (GHz), < 300 GHz
1071 : !
1072 : ! Returns:
1073 : ! 2-way gaseous attenuation (dB/km)
1074 : !
1075 : ! Reference:
1076 : ! Uses method of Liebe (1985)
1077 : !
1078 : ! Created:
1079 : ! 12/09/05 John Haynes (haynes@atmos.colostate.edu)
1080 : ! Modified:
1081 : ! 01/31/06 Port from IDL to Fortran 90
1082 : ! 12/19/14 Dustin Swales: Define type REALs as double precision (dustin.swales@noaa.gov)
1083 : ! ##############################################################################################
1084 :
1085 : ! INPUTS
1086 : real(wp), intent(in) :: & !
1087 : PRES_mb, & ! Pressure (mb) (hPa)
1088 : T, & ! Temperature (K)
1089 : SH, & ! Specific humidity
1090 : f ! Frequency (GHz), < 300 GHz
1091 :
1092 : ! PARAMETERS
1093 : integer, parameter :: & !
1094 : nbands_o2 = 48, & ! Number of O2 bands
1095 : nbands_h2o = 30 ! Number of h2o bands
1096 : ! LOCAL VARIABLES
1097 : real(wp) :: &
1098 : gases, th, e, p, sumo, gm0, a0, ap, term1, &
1099 : term2, term3, bf, be, term4, npp,e_th,one_th, &
1100 : pth3,eth35,aux1,aux2,aux3, aux4,gm,delt,x,y, &
1101 : gm2,fpp_o2,fpp_h2o,s_o2,s_h2o
1102 : integer :: i
1103 :
1104 : ! Table1 parameters v0, a1, a2, a3, a4, a5, a6
1105 : real(wp),dimension(nbands_o2),parameter :: &
1106 : v0 = (/49.4523790,49.9622570,50.4742380,50.9877480,51.5033500, &
1107 : 52.0214090,52.5423930,53.0669060,53.5957480,54.1299999,54.6711570, &
1108 : 55.2213650,55.7838000,56.2647770,56.3378700,56.9681000,57.6124810, &
1109 : 58.3238740,58.4465890,59.1642040,59.5909820,60.3060570,60.4347750, &
1110 : 61.1505580,61.8001520,62.4112120,62.4862530,62.9979740,63.5685150, &
1111 : 64.1277640,64.6789000,65.2240670,65.7647690,66.3020880,66.8368270, &
1112 : 67.3695950,67.9008620,68.4310010,68.9603060,69.4890210,70.0173420, &
1113 : 118.7503410,368.4983500,424.7631200,487.2493700,715.3931500, &
1114 : 773.8387300, 834.1453300/), &
1115 : a1 = (/0.0000001,0.0000003,0.0000009,0.0000025,0.0000061,0.0000141, &
1116 : 0.0000310,0.0000641,0.0001247,0.0002280,0.0003918,0.0006316,0.0009535, &
1117 : 0.0005489,0.0013440,0.0017630,0.0000213,0.0000239,0.0000146,0.0000240, &
1118 : 0.0000211,0.0000212,0.0000246,0.0000250,0.0000230,0.0000193,0.0000152, &
1119 : 0.0000150,0.0000109,0.0007335,0.0004635,0.0002748,0.0001530,0.0000801, &
1120 : 0.0000395,0.0000183,0.0000080,0.0000033,0.0000013,0.0000005,0.0000002, &
1121 : 0.0000094,0.0000679,0.0006380,0.0002350,0.0000996,0.0006710,0.0001800/),&
1122 : a2 = (/11.8300000,10.7200000,9.6900000,8.8900000,7.7400000,6.8400000, &
1123 : 6.0000000,5.2200000,4.4800000,3.8100000,3.1900000,2.6200000,2.1150000, &
1124 : 0.0100000,1.6550000,1.2550000,0.9100000,0.6210000,0.0790000,0.3860000, &
1125 : 0.2070000,0.2070000,0.3860000,0.6210000,0.9100000,1.2550000,0.0780000, &
1126 : 1.6600000,2.1100000,2.6200000,3.1900000,3.8100000,4.4800000,5.2200000, &
1127 : 6.0000000,6.8400000,7.7400000,8.6900000,9.6900000,10.7200000,11.8300000,&
1128 : 0.0000000,0.0200000,0.0110000,0.0110000,0.0890000,0.0790000,0.0790000/),&
1129 : a3 = (/0.0083000,0.0085000,0.0086000,0.0087000,0.0089000,0.0092000, &
1130 : 0.0094000,0.0097000,0.0100000,0.0102000,0.0105000,0.0107900,0.0111000, &
1131 : 0.0164600,0.0114400,0.0118100,0.0122100,0.0126600,0.0144900,0.0131900, &
1132 : 0.0136000,0.0138200,0.0129700,0.0124800,0.0120700,0.0117100,0.0146800, &
1133 : 0.0113900,0.0110800,0.0107800,0.0105000,0.0102000,0.0100000,0.0097000, &
1134 : 0.0094000,0.0092000,0.0089000,0.0087000,0.0086000,0.0085000,0.0084000, &
1135 : 0.0159200,0.0192000,0.0191600,0.0192000,0.0181000,0.0181000,0.0181000/),&
1136 : a4 = (/0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1137 : 0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1138 : 0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1139 : 0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1140 : 0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1141 : 0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
1142 : 0.0000000,0.6000000,0.6000000,0.6000000,0.6000000,0.6000000,0.6000000/),&
1143 : a5 = (/0.0056000,0.0056000,0.0056000,0.0055000,0.0056000,0.0055000, &
1144 : 0.0057000,0.0053000,0.0054000,0.0048000,0.0048000,0.0041700,0.0037500, &
1145 : 0.0077400,0.0029700,0.0021200,0.0009400,-0.0005500,0.0059700,-0.0024400,&
1146 : 0.0034400,-0.0041300,0.0013200,-0.0003600,-0.0015900,-0.0026600, &
1147 : -0.0047700,-0.0033400,-0.0041700,-0.0044800,-0.0051000,-0.0051000, &
1148 : -0.0057000,-0.0055000,-0.0059000,-0.0056000,-0.0058000,-0.0057000, &
1149 : -0.0056000,-0.0056000,-0.0056000,-0.0004400,0.0000000,0.0000000, &
1150 : 0.0000000,0.0000000,0.0000000,0.0000000/), &
1151 : a6 = (/1.7000000,1.7000000,1.7000000,1.7000000,1.8000000,1.8000000, &
1152 : 1.8000000,1.9000000,1.8000000,2.0000000,1.9000000,2.1000000,2.1000000, &
1153 : 0.9000000,2.3000000,2.5000000,3.7000000,-3.1000000,0.8000000,0.1000000, &
1154 : 0.5000000,0.7000000,-1.0000000,5.8000000,2.9000000,2.3000000,0.9000000, &
1155 : 2.2000000,2.0000000,2.0000000,1.8000000,1.9000000,1.8000000,1.8000000, &
1156 : 1.7000000,1.8000000,1.7000000,1.7000000,1.7000000,1.7000000,1.7000000, &
1157 : 0.9000000,1.0000000,1.0000000,1.0000000,1.0000000,1.0000000,1.0000000/)
1158 :
1159 : ! Table2 parameters v1, b1, b2, b3
1160 : real(wp),dimension(nbands_h2o),parameter :: &
1161 : v1 = (/22.2350800,67.8139600,119.9959400,183.3101170,321.2256440, &
1162 : 325.1529190,336.1870000,380.1973720,390.1345080,437.3466670,439.1508120, &
1163 : 443.0182950,448.0010750,470.8889740,474.6891270,488.4911330,503.5685320, &
1164 : 504.4826920,556.9360020,620.7008070,658.0065000,752.0332270,841.0735950, &
1165 : 859.8650000,899.4070000,902.5550000,906.2055240,916.1715820,970.3150220, &
1166 : 987.9267640/), &
1167 : b1 = (/0.1090000,0.0011000,0.0007000,2.3000000,0.0464000,1.5400000, &
1168 : 0.0010000,11.9000000,0.0044000,0.0637000,0.9210000,0.1940000,10.6000000, &
1169 : 0.3300000,1.2800000,0.2530000,0.0374000,0.0125000,510.0000000,5.0900000, &
1170 : 0.2740000,250.0000000,0.0130000,0.1330000,0.0550000,0.0380000,0.1830000, &
1171 : 8.5600000,9.1600000,138.0000000/), &
1172 : b2 = (/2.1430000,8.7300000,8.3470000,0.6530000,6.1560000,1.5150000, &
1173 : 9.8020000,1.0180000,7.3180000,5.0150000,3.5610000,5.0150000,1.3700000, &
1174 : 3.5610000,2.3420000,2.8140000,6.6930000,6.6930000,0.1140000,2.1500000, &
1175 : 7.7670000,0.3360000,8.1130000,7.9890000,7.8450000,8.3600000,5.0390000, &
1176 : 1.3690000,1.8420000,0.1780000/), &
1177 : b3 = (/0.0278400,0.0276000,0.0270000,0.0283500,0.0214000,0.0270000, &
1178 : 0.0265000,0.0276000,0.0190000,0.0137000,0.0164000,0.0144000,0.0238000, &
1179 : 0.0182000,0.0198000,0.0249000,0.0115000,0.0119000,0.0300000,0.0223000, &
1180 : 0.0300000,0.0286000,0.0141000,0.0286000,0.0286000,0.0264000,0.0234000, &
1181 : 0.0253000,0.0240000,0.0286000/)
1182 :
1183 : ! Conversions
1184 12247200 : th = 300._wp/T ! unitless
1185 :
1186 : ! DS2014 START: Using _wp for the exponential in the denominator results in slight errors
1187 : ! for dBze94. 0.01 % of values differ, relative range: 1.03e-05 to 1.78e-04
1188 : !e = (RH*th*th*th*th*th)/(41.45_wp*10**(9.834_wp*th-10)) ! kPa
1189 : !e = (RH*th*th*th*th*th)/(41.45_wp*10**(9.834_wp*th-10)) ! kPa
1190 12247200 : e = SH*PRES_mb/(SH+0.622_wp)/1000._wp !kPa
1191 : ! DS2014 END
1192 :
1193 12247200 : p = PRES_mb/1000._wp-e ! kPa
1194 12247200 : e_th = e*th
1195 12247200 : one_th = 1 - th
1196 12247200 : pth3 = p*th*th*th
1197 12247200 : eth35 = e*th**(3.5)
1198 :
1199 : ! Term1
1200 12247200 : sumo = 0._wp
1201 12247200 : aux1 = 1.1_wp*e_th
1202 600112800 : do i=1,nbands_o2
1203 587865600 : aux2 = f/v0(i)
1204 587865600 : aux3 = v0(i)-f
1205 587865600 : aux4 = v0(i)+f
1206 587865600 : gm = a3(i)*(p*th**(0.8_wp-a4(i))+aux1)
1207 587865600 : gm2 = gm*gm
1208 587865600 : delt = a5(i)*p*th**a6(i)
1209 587865600 : x = aux3*aux3+gm2
1210 587865600 : y = aux4*aux4+gm2
1211 587865600 : fpp_o2 = (((1._wp/x)+(1._wp/y))*(gm*aux2) - (delt*aux2)*((aux3/(x))-(aux4/(x))))
1212 587865600 : s_o2 = a1(i)*pth3*exp(a2(i)*one_th)
1213 600112800 : sumo = sumo + fpp_o2 * s_o2
1214 : enddo
1215 12247200 : term1 = sumo
1216 :
1217 : ! Term2
1218 12247200 : gm0 = 5.6E-3_wp*(p+1.1_wp*e)*th**(0.8_wp)
1219 12247200 : a0 = 3.07E-4_wp
1220 12247200 : ap = 1.4_wp*(1-1.2_wp*f**(1.5_wp)*1E-5)*1E-10
1221 12247200 : term2 = (2*a0*(gm0*(1+(f/gm0)*(f/gm0))*(1+(f/60._wp)**2))**(-1) + ap*p*th**(2.5_wp))*f*p*th*th
1222 :
1223 : ! Term3
1224 12247200 : sumo = 0._wp
1225 12247200 : aux1 = 4.8_wp*e_th
1226 379663200 : do i=1,nbands_h2o
1227 367416000 : aux2 = f/v1(i)
1228 367416000 : aux3 = v1(i)-f
1229 367416000 : aux4 = v1(i)+f
1230 367416000 : gm = b3(i)*(p*th**(0.8)+aux1)
1231 367416000 : gm2 = gm*gm
1232 367416000 : x = aux3*aux3+gm2
1233 367416000 : y = aux4*aux4+gm2
1234 367416000 : fpp_h2o = ((1._wp/x)+(1._wp/y))*(gm*aux2) ! - (delt*aux2)*((aux3/(x))-(aux4/(x)))
1235 367416000 : s_h2o = b1(i)*eth35*exp(b2(i)*one_th)
1236 379663200 : sumo = sumo + fpp_h2o * s_h2o
1237 : enddo
1238 12247200 : term3 = sumo
1239 :
1240 : ! Term4
1241 12247200 : bf = 1.4E-6_wp
1242 12247200 : be = 5.41E-5_wp
1243 12247200 : term4 = (bf*p+be*e*th*th*th)*f*e*th**(2.5_wp)
1244 :
1245 : ! Summation and result
1246 12247200 : npp = term1 + term2 + term3 + term4
1247 12247200 : gases = 0.182_wp*f*npp
1248 :
1249 12247200 : end function gases
1250 1536 : subroutine hydro_class_init(lsingle,ldouble,sd)
1251 : ! ##############################################################################################
1252 : ! Purpose:
1253 : !
1254 : ! Initialize variables used by the radar simulator.
1255 : ! Part of QuickBeam v3.0 by John Haynes and Roj Marchand
1256 : !
1257 : ! Inputs:
1258 : ! NAME SIZE DESCRIPTION
1259 : ! [lsingle] (1) Logical flag to use single moment
1260 : ! [ldouble] (1) Logical flag to use two moment
1261 : ! Outputs:
1262 : ! [sd] Structure that define hydrometeor types
1263 : !
1264 : ! Local variables:
1265 : ! [n_hydro] (1) Number of hydrometeor types
1266 : ! [hclass_type] (nhclass) Type of distribution (see quickbeam documentation)
1267 : ! [hclass_phase] (nhclass) 1==ice, 0=liquid
1268 : ! [hclass_dmin] (nhclass) Minimum diameter allowed is drop size distribution N(D<Dmin)=0
1269 : ! [hclass_dmax] (nhclass) Maximum diameter allowed is drop size distribution N(D>Dmax)=0
1270 : ! [hclass_apm] (nhclass) Density of partical apm*D^bpm or constant = rho
1271 : ! [hclass_bpm] (nhclass) Density of partical apm*D^bpm or constant = rho
1272 : ! [hclass_rho] (nhclass) Density of partical apm*D^bpm or constant = rho
1273 : ! [hclass_p1] (nhclass) Default values of DSD parameters (see quickbeam documentation)
1274 : ! [hclass_p2] (nhclass) Default values of DSD parameters (see quickbeam documentation)
1275 : ! [hclass_p3] (nhclass) Default values of DSD parameters (see quickbeam documentation)
1276 : ! Modified:
1277 : ! 08/23/2006 placed into subroutine form (Roger Marchand)
1278 : ! June 2010 New interface to support "radar_simulator_params" structure
1279 : ! 12/22/2014 Moved radar simulator (CLOUDSAT) configuration initialization to cloudsat_init
1280 : ! ##############################################################################################
1281 :
1282 : ! ####################################################################################
1283 : ! NOTES on HCLASS variables
1284 : !
1285 : ! TYPE - Set to
1286 : ! 1 for modified gamma distribution,
1287 : ! 2 for exponential distribution,
1288 : ! 3 for power law distribution,
1289 : ! 4 for monodisperse distribution,
1290 : ! 5 for lognormal distribution.
1291 : !
1292 : ! PHASE - Set to 0 for liquid, 1 for ice.
1293 : ! DMIN - The minimum drop size for this class (micron), ignored for monodisperse.
1294 : ! DMAX - The maximum drop size for this class (micron), ignored for monodisperse.
1295 : ! Important note: The settings for DMIN and DMAX are
1296 : ! ignored in the current version for all distributions except for power
1297 : ! law. Except when the power law distribution is used, particle size
1298 : ! is fixed to vary from zero to infinity, a restriction that is expected
1299 : ! to be lifted in future versions. A placeholder must still be specified
1300 : ! for each.
1301 : ! Density of particles is given by apm*D^bpm or a fixed value rho. ONLY specify ONE of these two!!
1302 : ! APM - The alpha_m coefficient in equation (1) (kg m**-beta_m )
1303 : ! BPM - The beta_m coefficient in equation (1), see section 4.1.
1304 : ! RHO - Hydrometeor density (kg m-3 ).
1305 : !
1306 : ! P1, P2, P3 - are default distribution parameters that depend on the type
1307 : ! of distribution (see quickmbeam documentation for more information)
1308 : !
1309 : ! Modified Gamma (must set P3 and one of P1 or P2)
1310 : ! P1 - Set to the total particle number concentration Nt /rho_a (kg-1 ), where
1311 : ! rho_a is the density of air in the radar volume.
1312 : ! P2 - Set to the particle mean diameter D (micron).
1313 : ! P3 - Set to the distribution width nu.
1314 : !
1315 : ! Exponetial (set one of)
1316 : ! P1 - Set to a constant intercept parameter N0 (m-4).
1317 : ! P2 - Set to a constant lambda (micron-1).
1318 : !
1319 : ! Power Law
1320 : ! P1 - Set this to the value of a constant power law parameter br
1321 : !
1322 : ! Monodisperse
1323 : ! P1 - Set to a constant diameter D0 (micron) = Re.
1324 : !
1325 : ! Log-normal (must set P3 and one of P1 or P2)
1326 : ! P1 - Set to the total particle number concentration Nt /rho_a (kg-1 )
1327 : ! P2 - Set to the geometric mean particle radius rg (micron).
1328 : ! P3 - Set to the natural logarithm of the geometric standard deviation.
1329 : ! ####################################################################################
1330 : ! INPUTS
1331 : logical,intent(in) :: &
1332 : lsingle, & ! True -> use single moment
1333 : ldouble ! True -> use two moment
1334 :
1335 : ! OUTPUTS
1336 : type(size_distribution),intent(out) ::&
1337 : sd !
1338 :
1339 : ! SINGLE MOMENT PARAMETERS
1340 : integer,parameter,dimension(N_HYDRO) :: &
1341 : ! LSL LSI LSR LSS CVL CVI CVR CVS LSG
1342 : HCLASS1_TYPE = (/5, 1, 2, 2, 5, 1, 2, 2, 2/), & !
1343 : HCLASS1_PHASE = (/0, 1, 0, 1, 0, 1, 0, 1, 1/) !
1344 : real(wp),parameter,dimension(N_HYDRO) ::&
1345 : ! LSL LSI LSR LSS CVL CVI CVR CVS LSG
1346 : HCLASS1_DMIN = (/ -1., -1., -1., -1., -1., -1., -1., -1., -1. /), &
1347 : HCLASS1_DMAX = (/ -1., -1., -1., -1., -1., -1., -1., -1., -1. /), &
1348 : HCLASS1_APM = (/524., 110.8, 524., -1., 524., 110.8, 524., -1., -1. /), &
1349 : HCLASS1_BPM = (/ 3., 2.91, 3., -1., 3., 2.91, 3., -1., -1. /), &
1350 : HCLASS1_RHO = (/ -1., -1., -1., 100., -1., -1., -1., 100., 400. /), &
1351 : HCLASS1_P1 = (/ -1., -1., 8.e6, 3.e6, -1., -1., 8.e6, 3.e6, 4.e6/), &
1352 : HCLASS1_P2 = (/ 6., 40., -1., -1., 6., 40., -1., -1., -1. /), &
1353 : HCLASS1_P3 = (/ 0.3, 2., -1., -1., 0.3, 2., -1., -1., -1. /)
1354 :
1355 : ! TWO MOMENT PARAMETERS
1356 : integer,parameter,dimension(N_HYDRO) :: &
1357 : ! LSL LSI LSR LSS CVL CVI CVR CVS LSG
1358 : HCLASS2_TYPE = (/ 1, 1, 1, 1, 1, 1, 1, 1, 1/), &
1359 : HCLASS2_PHASE = (/ 0, 1, 0, 1, 0, 1, 0, 1, 1/)
1360 :
1361 : real(wp),parameter,dimension(N_HYDRO) :: &
1362 : ! LSL LSI LSR LSS CVL CVI CVR CVS LSG
1363 : HCLASS2_DMIN = (/ -1, -1, -1, -1, -1, -1, -1, -1, -1/), &
1364 : HCLASS2_DMAX = (/ -1, -1, -1, -1, -1, -1, -1, -1, -1/), &
1365 : HCLASS2_APM = (/524, -1, 524, -1, 524, -1, 524, -1, -1/), &
1366 : HCLASS2_BPM = (/ 3, -1, 3, -1, 3, -1, 3, -1, -1/), &
1367 : HCLASS2_RHO = (/ -1, 500, -1, 100, -1, 500, -1, 100, 900/), &
1368 : HCLASS2_P1 = (/ -1, -1, -1, -1, -1, -1, -1, -1, -1/), &
1369 : HCLASS2_P2 = (/ -1, -1, -1, -1, -1, -1, -1, -1, -1/), &
1370 : HCLASS2_P3 = (/ -2, 1, 1, 1, -2, 1, 1, 1, 1/)
1371 :
1372 1536 : if (lsingle) then
1373 15360 : sd%dtype(1:N_HYDRO) = HCLASS1_TYPE(1:N_HYDRO)
1374 15360 : sd%phase(1:N_HYDRO) = HCLASS1_PHASE(1:N_HYDRO)
1375 15360 : sd%dmin(1:N_HYDRO) = HCLASS1_DMIN(1:N_HYDRO)
1376 15360 : sd%dmax(1:N_HYDRO) = HCLASS1_DMAX(1:N_HYDRO)
1377 15360 : sd%apm(1:N_HYDRO) = HCLASS1_APM(1:N_HYDRO)
1378 15360 : sd%bpm(1:N_HYDRO) = HCLASS1_BPM(1:N_HYDRO)
1379 15360 : sd%rho(1:N_HYDRO) = HCLASS1_RHO(1:N_HYDRO)
1380 15360 : sd%p1(1:N_HYDRO) = HCLASS1_P1(1:N_HYDRO)
1381 15360 : sd%p2(1:N_HYDRO) = HCLASS1_P2(1:N_HYDRO)
1382 15360 : sd%p3(1:N_HYDRO) = HCLASS1_P3(1:N_HYDRO)
1383 : endif
1384 1536 : if (ldouble) then
1385 0 : sd%dtype(1:N_HYDRO) = HCLASS2_TYPE(1:N_HYDRO)
1386 0 : sd%phase(1:N_HYDRO) = HCLASS2_PHASE(1:N_HYDRO)
1387 0 : sd%dmin(1:N_HYDRO) = HCLASS2_DMIN(1:N_HYDRO)
1388 0 : sd%dmax(1:N_HYDRO) = HCLASS2_DMAX(1:N_HYDRO)
1389 0 : sd%apm(1:N_HYDRO) = HCLASS2_APM(1:N_HYDRO)
1390 0 : sd%bpm(1:N_HYDRO) = HCLASS2_BPM(1:N_HYDRO)
1391 0 : sd%rho(1:N_HYDRO) = HCLASS2_RHO(1:N_HYDRO)
1392 0 : sd%p1(1:N_HYDRO) = HCLASS2_P1(1:N_HYDRO)
1393 0 : sd%p2(1:N_HYDRO) = HCLASS2_P2(1:N_HYDRO)
1394 0 : sd%p3(1:N_HYDRO) = HCLASS2_P3(1:N_HYDRO)
1395 : endif
1396 1536 : end subroutine hydro_class_init
1397 0 : end module mod_quickbeam_optics
|
