source: LMDZ4/branches/LMDZ4V5.0-LF/libf/cosp/radar_simulator.F90 @ 4082

Last change on this file since 4082 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 = ((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 
Note: See TracBrowser for help on using the repository browser.