Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

radar_simulator.F90 @ 5429

Last change on this file since 5429 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
File size: 15.5 KB

Line
1	subroutine radar_simulator(freq,k2,do_ray,use_gas_abs,use_mie_table,mt, &
2	nhclass,hp,nprof,ngate,nsizes,D,hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix, &
3	rh_matrix,Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, &
4	g_to_vol_in,g_to_vol_out)
5
6	! rh_matrix,Ze_non,Ze_ray,kr_matrix,g_atten_to_vol,dBZe)
7
8	use m_mrgrnk
9	use array_lib
10	use math_lib
11	use optics_lib
12	use radar_simulator_types
13	implicit none
14
15	! Purpose:
16	! Simulates a vertical profile of radar reflectivity
17	! Part of QuickBeam v1.04 by John Haynes & Roger Marchand
18	!
19	! Inputs:
20	! [freq] radar frequency (GHz), can be anything unless
21	! use_mie_table=1, in which case one of 94,35,13.8,9.6,3
22	! [k2] \|K\|^2, the dielectric constant, set to -1 to use the
23	! frequency dependent default
24	! [do_ray] 1=do Rayleigh calcs, 0=not
25	! [use_gas_abs] 1=do gaseous abs calcs, 0=not,
26	! 2=use same as first profile (undocumented)
27	! [use_mie_table] 1=use Mie tables, 0=not
28	! [mt] Mie look up table
29	! [nhclass] number of hydrometeor types
30	! [hp] structure that defines hydrometeor types
31	! [nprof] number of hydrometeor profiles
32	! [ngate] number of vertical layers
33	! [nsizes] number of discrete particles in [D]
34	! [D] array of discrete particles (um)
35	!
36	! (The following 5 arrays must be in order from closest to the radar
37	! to farthest...)
38	! [hgt_matrix] height of hydrometeors (km)
39	! [hm_matrix] table of hydrometeor mixing rations (g/kg)
40	! [re_matrix] OPTIONAL table of hydrometeor effective radii (microns)
41	! [p_matrix] pressure profile (hPa)
42	! [t_matrix] temperature profile (C)
43	! [rh_matrix] relative humidity profile (%)
44	!
45	! Outputs:
46	! [Ze_non] radar reflectivity without attenuation (dBZ)
47	! [Ze_ray] Rayleigh reflectivity (dBZ)
48	! [h_atten_to_vol] attenuation by hydromets, radar to vol (dB)
49	! [g_atten_to_vol] gaseous atteunation, radar to vol (dB)
50	! [dBZe] effective radar reflectivity factor (dBZ)
51	!
52	! Optional:
53	! [g_to_vol_in] integrated atten due to gases, r>v (dB).
54	! If present then is used as gaseous absorption, independently of the
55	! value in use_gas_abs
56	! [g_to_vol_out] integrated atten due to gases, r>v (dB).
57	! If present then gaseous absorption for each profile is returned here.
58	!
59	! Created:
60	! 11/28/2005 John Haynes (haynes@atmos.colostate.edu)
61	! Modified:
62	! 09/2006 placed into subroutine form, scaling factors (Roger Marchand,JMH)
63	! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
64	! 01/2008 'Do while' to determine if hydrometeor(s) present in volume
65	! changed for vectorization purposes (A. Bodas-Salcedo)
66
67	! ----- INPUTS -----
68	type(mie), intent(in) :: mt
69	type(class_param), intent(inout) :: hp
70	real*8, intent(in) :: freq,k2
71	integer, intent(in) :: do_ray,use_gas_abs,use_mie_table, &
72	nhclass,nprof,ngate,nsizes
73	real*8, dimension(nsizes), intent(in) :: D
74	real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
75	t_matrix,rh_matrix
76	real*8, dimension(nhclass,nprof,ngate), intent(in) :: hm_matrix
77	real*8, dimension(nhclass,nprof,ngate), intent(inout) :: re_matrix
78
79	! ----- OUTPUTS -----
80	real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
81	g_atten_to_vol,dBZe,h_atten_to_vol
82
83	! ----- OPTIONAL -----
84	real*8, optional, dimension(ngate,nprof) :: &
85	g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
86	! the same gaseous absorption in different calls. Optional to allow compatibility
87	! with original version. A. Bodas April 2008.
88
89	! real*8, dimension(nprof,ngate) :: kr_matrix
90
91	! ----- INTERNAL -----
92	integer :: &
93	phase, & ! 0=liquid, 1=ice
94	ns ! number of discrete drop sizes
95
96	integer*4, dimension(ngate) :: &
97	hydro ! 1=hydrometeor in vol, 0=none
98	real*8 :: &
99	rho_a, & ! air density (kg m^-3)
100	gases ! function: 2-way gas atten (dB/km)
101
102	real*8, dimension(:), allocatable :: &
103	Di, Deq, & ! discrete drop sizes (um)
104	Ni, Ntemp, & ! discrete concentrations (cm^-3 um^-1)
105	rhoi ! discrete densities (kg m^-3)
106
107	real*8, dimension(ngate) :: &
108	z_vol, & ! effective reflectivity factor (mm^6/m^3)
109	z_ray, & ! reflectivity factor, Rayleigh only (mm^6/m^3)
110	kr_vol, & ! attenuation coefficient hydro (dB/km)
111	g_vol, & ! attenuation coefficient gases (dB/km)
112	a_to_vol, & ! integrated atten due to hydometeors, r>v (dB)
113	g_to_vol ! integrated atten due to gases, r>v (dB)
114
115
116	integer,parameter :: KR8 = selected_real_kind(15,300)
117	real*8, parameter :: xx = -1.0_KR8
118	real*8, dimension(:), allocatable :: xxa
119	real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2,apm,bpm
120	integer*4 :: tp, i, j, k, pr, itt, iff
121
122	real*8 bin_length,step,base,step_list(25),base_list(25)
123	integer*4 iRe_type,n,max_bin
124
125	logical :: g_to_vol_in_present, g_to_vol_out_present
126
127	! Logicals to avoid calling present within the loops
128	g_to_vol_in_present = present(g_to_vol_in)
129	g_to_vol_out_present = present(g_to_vol_out)
130
131	! set up Re bins for z_scalling
132	bin_length=50;
133	max_bin=25
134
135	step_list(1)=1
136	base_list(1)=75
137	do j=2,max_bin
138	step_list(j)=3*(j-1);
139	if(step_list(j)>bin_length) then
140	step_list(j)=bin_length;
141	endif
142	base_list(j)=base_list(j-1)+floor(bin_length/step_list(j-1));
143	enddo
144
145
146	pi = acos(-1.0)
147	if (use_mie_table == 1) iff = infind(mt%freq,freq,sort=1)
148
149
150	! // loop over each profile (nprof)
151	do pr=1,nprof
152
153	! ----- calculations for each volume -----
154	z_vol(:) = 0
155	z_ray(:) = 0
156	kr_vol(:) = 0
157	hydro(:) = 0
158
159	! // loop over eacho range gate (ngate)
160	do k=1,ngate
161
162	! :: determine if hydrometeor(s) present in volume
163	hydro(k) = 0
164	do j=1,nhclass ! Do while changed for vectorization purposes (A. B-S)
165	if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
166	hydro(k) = 1
167	exit
168	endif
169	enddo
170
171	if (hydro(k) == 1) then
172	! :: if there is hydrometeor in the volume
173
174	rho_a = (p_matrix(pr,k)100.)/(287(t_matrix(pr,k)+273.15))
175
176	! :: loop over hydrometeor type
177	do tp=1,nhclass
178
179	if (hm_matrix(tp,pr,k) <= 1E-12) cycle
180
181	phase = hp%phase(tp)
182	if(phase==0) then
183	itt = infind(mt_ttl,t_matrix(pr,k))
184	else
185	itt = infind(mt_tti,t_matrix(pr,k))
186	endif
187
188	! calculate Re if we have an exponential distribution with fixed No ... precipitation type particle
189	if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8) then
190
191	apm=hp%apm(tp)
192	bpm=hp%bpm(tp)
193
194	if ((hp%rho(tp) > 0) .and. (apm < 0)) then
195	apm = (pi/6)*hp%rho(tp)
196	bpm = 3.
197	endif
198
199	tmp1 = 1./(1.+bpm)
200	ld = ((apmgamma(1.+bpm)hp%p1(tp))/(rho_ahm_matrix(tp,pr,k)1E-3))**tmp1
201
202	Re = 1.5E6/ld
203
204	re_matrix(tp,pr,k) = Re;
205
206	endif
207
208	if(re_matrix(tp,pr,k).eq.0) then
209
210	iRe_type=1
211	Re=0
212	else
213	iRe_type=1
214	Re=re_matrix(tp,pr,k)
215
216	n=floor(Re/bin_length)
217	if(n==0) then
218	if(Re<25) then
219	step=0.5
220	base=0
221	else
222	step=1
223	base=25
224	endif
225	else
226	if(n>max_bin) then
227	n=max_bin
228	endif
229
230	step=step_list(n)
231	base=base_list(n)
232	endif
233
234	iRe_type=floor(Re/step)
235
236	if(iRe_type.lt.1) then
237	iRe_type=1
238	endif
239
240	Re=step*(iRe_type+0.5)
241	iRe_type=iRe_type+base-floor(n*bin_length/step)
242
243	! make sure iRe_type is within bounds
244	if(iRe_type.ge.nRe_types) then
245
246	! print *, tp, re_matrix(tp,pr,k), Re, iRe_type
247
248	! no scaling allowed
249	Re=re_matrix(tp,pr,k)
250
251	iRe_type=nRe_types
252	hp%z_flag(tp,itt,iRe_type)=.false.
253	hp%scaled(tp,iRe_type)=.false.
254	endif
255	endif
256
257	! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
258	! if not we will calculate Ze, Zr, and Kr from the distribution parameters
259	if( .not. hp%z_flag(tp,itt,iRe_type) ) then
260
261	! :: create a distribution of hydrometeors within volume
262	select case(hp%dtype(tp))
263	case(4)
264	ns = 1
265	allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
266	if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns))
267	Di = hp%p1(tp)
268	Ni = 0.
269	case default
270	ns = nsizes
271	allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
272	if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns))
273	Di = D
274	Ni = 0.
275	end select
276
277	! :: create a DSD (using scaling factor if applicable)
278	! hp%scaled(tp,iRe_type)=.false. ! turn off N scaling
279
280	call dsd(hm_matrix(tp,pr,k),Re,Di,Ni,ns,hp%dtype(tp),rho_a, &
281	t_matrix(pr,k),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
282	hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),hp%fc(tp,1:ns,iRe_type), &
283	hp%scaled(tp,iRe_type))
284
285	! :: calculate particle density
286	! if ((hp%rho_eff(tp,1,iRe_type) < 0) .and. (phase == 1)) then
287	if (phase == 1) then
288	if (hp%rho(tp) < 0) then
289
290	! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density
291	! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)hp%apm(tp)(Di1E-6)*(hp%bpm(tp)-3) !MG Mie approach
292
293	! as the particle size gets small it is possible that the mass to size relationship of
294	! (given by power law in hclass.data) can produce impossible results
295	! where the mass is larger than a solid sphere of ice.
296	! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
297	! do i=1,ns
298	! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
299	! hp%rho_eff(tp,i,iRe_type) = 917
300	!endif
301	!enddo
302
303	! alternative is to use equivalent volume spheres.
304	hp%rho_eff(tp,1:ns,iRe_type) = 917 ! solid ice == equivalent volume approach
305	Deq = ( ( 6/pihp%apm(tp)/917 ) * (1.0/3.0) ) * &
306	( (Di1E-6) * (hp%bpm(tp)/3.0) ) * 1E6 ! Di now really Deq in microns.
307
308	else
309
310	! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp) !MG Mie approach
311
312	! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
313	hp%rho_eff(tp,1:ns,iRe_type) = 917
314	Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0))
315
316	endif
317
318	! if using equivalent volume spheres
319	if (use_mie_table == 1) then
320
321	Ntemp=Ni
322
323	! Find N(Di) from N(Deq) which we know
324	do i=1,ns
325	j=infind(Deq,Di(i))
326	Ni(i)=Ntemp(j)
327	enddo
328	else
329	! just use Deq and D variable input to mie code
330	Di=Deq;
331	endif
332
333	endif
334	rhoi = hp%rho_eff(tp,1:ns,iRe_type)
335
336	! :: calculate effective reflectivity factor of volume
337	if (use_mie_table == 1) then
338
339	if ((hp%dtype(tp) == 4) .and. (hp%idd(tp) < 0)) then
340	hp%idd(tp) = infind(mt%D,Di(1))
341	endif
342
343	if (phase == 0) then
344
345	! itt = infind(mt_ttl,t_matrix(pr,k))
346	select case(hp%dtype(tp))
347	case(4)
348	mt_qext(1) = mt%qext(hp%idd(tp),itt,1,iff)
349	mt_qbsca(1) = mt%qbsca(hp%idd(tp),itt,1,iff)
350	case default
351	mt_qext = mt%qext(:,itt,1,iff)
352	mt_qbsca = mt%qbsca(:,itt,1,iff)
353	end select
354
355	call zeff(freq,Di,Ni,ns,k2,mt_ttl(itt),0,do_ray, &
356	ze,zr,kr,mt_qext,mt_qbsca,xx)
357
358	else
359
360	! itt = infind(mt_tti,t_matrix(pr,k))
361	select case(hp%dtype(tp))
362	case(4)
363	if (hp%ifc(tp,1,iRe_type) < 0) then
364	hp%ifc(tp,1,iRe_type) = infind(mt%f,rhoi(1)/917.)
365	endif
366	mt_qext(1) = &
367	mt%qext(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff)
368	mt_qbsca(1) = &
369	mt%qbsca(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff)
370	case default
371	do i=1,ns
372	if (hp%ifc(tp,i,iRe_type) < 0) then
373	hp%ifc(tp,i,iRe_type) = infind(mt%f,rhoi(i)/917.)
374	endif
375	mt_qext(i) = mt%qext(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff)
376	mt_qbsca(i) = mt%qbsca(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff)
377	enddo
378	end select
379
380	call zeff(freq,Di,Ni,ns,k2,mt_tti(itt),1,do_ray, &
381	ze,zr,kr,mt_qext,mt_qbsca,xx)
382
383	endif
384
385	else
386
387	xxa = -9.9
388	call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, &
389	ze,zr,kr,xxa,xxa,rhoi)
390
391
392	endif ! end of use mie table
393
394	! xxa = -9.9
395	!call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, &
396	! ze2,zr,kr2,xxa,xxa,rhoi)
397
398	! if(abs(ze2-ze)/ze2 > 0.1) then
399	! if(abs(kr2-kr)/kr2 > 0.1) then
400
401	! write(,) pr,k,tp,ze2,ze2-ze,abs(ze2-ze)/ze2,itt+cnt_liq,iff
402	! write(,) pr,k,tp,ze2,kr2,kr2-kr,abs(kr2-kr)/kr2
403	! stop
404
405	!endif
406
407	deallocate(Di,Ni,rhoi,xxa,Deq)
408	if (use_mie_table == 1) deallocate(mt_qext,mt_qbsca,Ntemp)
409
410	else ! can use z scaling
411
412	if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then
413
414	ze = hp%Ze_scaled(tp,itt,iRe_type)
415	zr = hp%Zr_scaled(tp,itt,iRe_type)
416	kr = hp%kr_scaled(tp,itt,iRe_type)
417
418	else
419	scale_factor=rho_a*hm_matrix(tp,pr,k)
420
421	zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor
422	ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
423	kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor
424	endif
425
426	endif ! end z_scaling
427
428	! kr=0
429
430	kr_vol(k) = kr_vol(k) + kr
431	z_vol(k) = z_vol(k) + ze
432	z_ray(k) = z_ray(k) + zr
433
434	! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
435	if( .not. hp%z_flag(tp,itt,iRe_type) .and. 1.eq.1 ) then
436
437	if( ( (hp%dtype(tp)==1 .or. hp%dtype(tp)==5 .or. hp%dtype(tp)==2) .and. abs(hp%p1(tp)+1) < 1E-8 ) .or. &
438	( hp%dtype(tp)==3 .or. hp%dtype(tp)==4 ) &
439	) then
440
441	scale_factor=rho_a*hm_matrix(tp,pr,k)
442
443	hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
444	hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
445	hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
446
447	hp%z_flag(tp,itt,iRe_type)=.True.
448
449	elseif( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then
450
451	hp%Ze_scaled(tp,itt,iRe_type) = ze
452	hp%Zr_scaled(tp,itt,iRe_type) = zr
453	hp%kr_scaled(tp,itt,iRe_type) = kr
454
455	hp%z_flag(tp,itt,iRe_type)=.True.
456	endif
457
458	endif
459
460	enddo ! end loop of tp (hydrometeor type)
461
462	else
463	! :: volume is hydrometeor-free
464
465	kr_vol(k) = 0
466	z_vol(k) = -999
467	z_ray(k) = -999
468
469	endif
470
471	! :: attenuation due to hydrometeors between radar and volume
472	a_to_vol(k) = 2*path_integral(kr_vol,hgt_matrix(pr,:),1,k-1)
473
474	! :: attenuation due to gaseous absorption between radar and volume
475	if (g_to_vol_in_present) then
476	g_to_vol(k) = g_to_vol_in(k,pr)
477	else
478	if ( (use_gas_abs == 1) .or. ((use_gas_abs == 2) .and. (pr == 1)) ) then
479	g_vol(k) = gases(p_matrix(pr,k),t_matrix(pr,k)+273.15, &
480	rh_matrix(pr,k),freq)
481	g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1)
482	elseif (use_gas_abs == 0) then
483	g_to_vol(k) = 0
484	endif
485	endif
486
487	! kr_matrix(pr,:)=kr_vol
488
489	! :: store results in matrix for return to calling program
490	h_atten_to_vol(pr,k)=a_to_vol(k)
491	g_atten_to_vol(pr,k)=g_to_vol(k)
492	if ((do_ray == 1) .and. (z_ray(k) > 0)) then
493	Ze_ray(pr,k) = 10*log10(z_ray(k))
494	else
495	Ze_ray(pr,k) = -999
496	endif
497	if (z_vol(k) > 0) then
498	dBZe(pr,k) = 10*log10(z_vol(k))-a_to_vol(k)-g_to_vol(k)
499	Ze_non(pr,k) = 10*log10(z_vol(k))
500	else
501	dBZe(pr,k) = -999
502	Ze_non(pr,k) = -999
503	endif
504
505	enddo ! end loop of k (range gate)
506	! Output array with gaseous absorption
507	if (g_to_vol_out_present) g_to_vol_out(:,pr) = g_to_vol
508	enddo ! end loop over pr (profile)
509
510	end subroutine radar_simulator
511

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format