My Project
 All Classes Files Functions Variables Macros
radar_simulator.F90
Go to the documentation of this file.
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 = ((apm*gamma(1.+bpm)*hp%p1(tp))/(rho_a*hm_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)*(Di*1E-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/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * &
306  ( (di*1e-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