source: LMDZ5/trunk/libf/phylmd/cosp/radar_simulator.F90 @ 3365

Last change on this file since 3365 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 17.9 KB
Line 
1  subroutine radar_simulator( &
2    hp, &
3    nprof,ngate, &
4    undef, &
5    hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
6    p_matrix,t_matrix,rh_matrix, &
7    Ze_non,Ze_ray,a_to_vol,g_to_vol,dBZe, &
8    g_to_vol_in,g_to_vol_out)
9 
10  use m_mrgrnk
11  use array_lib
12  use math_lib
13  use optics_lib
14  use radar_simulator_types
15  use scale_LUTs_io
16  implicit none
17 
18! Purpose:
19!
20!   Simulates a vertical profile of radar reflectivity
21!   Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand.
22!   but has been substantially modified since that time by
23!   Laura Fowler and Roger Marchand (see modifications below).
24!
25! Inputs:
26!
27!   [hp]              structure that defines hydrometeor types and other radar properties
28!
29!   [nprof]           number of hydrometeor profiles
30!   [ngate]           number of vertical layers
31!
32!   [undef]           missing data value
33!   (The following 5 arrays must be in order from closest to the radar
34!    to farthest...)
35!
36!   [hgt_matrix]      height of hydrometeors (km)
37!   [p_matrix]        pressure profile (hPa)
38!   [t_matrix]        temperature profile (K)
39!   [rh_matrix]       relative humidity profile (%) -- only needed if gaseous aborption calculated.
40!
41!   [hm_matrix]       table of hydrometeor mixing rations (g/kg)
42!   [re_matrix]       table of hydrometeor effective radii.  0 ==> use defaults. (units=microns)   
43!   [Np_matrix]       table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
44!
45! Outputs:
46!
47!   [Ze_non]          radar reflectivity without attenuation (dBZ)
48!   [Ze_ray]          Rayleigh reflectivity (dBZ)
49!   [h_atten_to_vol]  attenuation by hydromets, radar to vol (dB)
50!   [g_atten_to_vol]  gaseous atteunation, radar to vol (dB)
51!   [dBZe]            effective radar reflectivity factor (dBZ)
52!
53! Optional:
54!   [g_to_vol_in]     integrated atten due to gases, r>v (dB).
55!                     If present then is used as gaseous absorption, independently of the
56!                     value in use_gas_abs
57!   [g_to_vol_out]    integrated atten due to gases, r>v (dB).
58!                     If present then gaseous absorption for each profile is returned here.
59!
60! Created:
61!   11/28/2005  John Haynes (haynes@atmos.colostate.edu)
62!
63! Modified:
64!   09/2006  placed into subroutine form (Roger Marchand,JMH)
65!   08/2007  added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
66!   01/2008  'Do while' to determine if hydrometeor(s) present in volume
67!             changed for vectorization purposes (A. Bodas-Salcedo)
68!
69!   07/2010  V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT)
70!  ... All hydrometeor and radar simulator properties now included in hp structure
71!  ... hp structure should be initialized by call to radar_simulator_init prior
72!  ... to calling this subroutine. 
73!     Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added
74!  ... Changes implement by Roj Marchand following work by Laura Fowler
75!
76!   10/2011  Modified ngate loop to go in either direction depending on flag
77!     hp%radar_at_layer_one.  This affects the direction in which attenuation is summed.
78!
79!     Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple
80!     summation. (Roger Marchand)
81!
82!
83! ----- INPUTS ----- 
84 
85  logical, parameter  ::  DO_LUT_TEST = .false.
86  logical, parameter  ::  DO_NP_TEST = .false.
87
88  type(class_param), intent(inout) :: hp
89
90  integer, intent(in) ::  nprof,ngate
91
92  real undef
93  real*8, dimension(nprof,ngate), intent(in) :: &
94    hgt_matrix, p_matrix,t_matrix,rh_matrix
95
96  real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
97  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
98  real*8, dimension(hp%nhclass,nprof,ngate), intent(in)    :: Np_matrix
99
100! ----- OUTPUTS -----
101  real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
102       g_to_vol,dBZe,a_to_vol
103
104! ----- OPTIONAL -----
105  real*8, optional, dimension(nprof,ngate) :: &
106  g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
107                           ! the same gaseous absorption in different calls. Optional to allow compatibility
108                           ! with original version. A. Bodas April 2008.
109!  real*8, dimension(nprof,ngate) :: kr_matrix
110
111! ----- INTERNAL -----
112
113  real, parameter :: one_third = 1.0/3.0
114  real*8 :: t_kelvin
115  integer :: &
116  phase, & ! 0=liquid, 1=ice
117  ns       ! number of discrete drop sizes
118
119  logical :: hydro      ! true=hydrometeor in vol, false=none
120  real*8 :: &
121  rho_a, &   ! air density (kg m^-3)
122  gases      ! function: 2-way gas atten (dB/km)
123
124  real*8, dimension(:), allocatable :: &
125  Di, Deq, &   ! discrete drop sizes (um)
126  Ni, &        ! discrete concentrations (cm^-3 um^-1)
127  rhoi         ! discrete densities (kg m^-3)
128
129  real*8, dimension(nprof, ngate) :: &
130  z_vol, &      ! effective reflectivity factor (mm^6/m^3)
131  z_ray, &                      ! reflectivity factor, Rayleigh only (mm^6/m^3)
132  kr_vol, &     ! attenuation coefficient hydro (dB/km)
133  g_vol         ! attenuation coefficient gases (dB/km)
134
135
136  integer,parameter :: KR8 = selected_real_kind(15,300)
137  real*8, parameter :: xx = -1.0_KR8
138  real*8,  dimension(:), allocatable :: xxa
139  real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
140  real*8 :: half_a_atten_current,half_a_atten_above
141  real*8 :: half_g_atten_current,half_g_atten_above
142  integer*4 :: tp, i, j, k, pr, itt, iff
143
144  real*8    step,base, Np
145  integer*4 iRe_type,n,max_bin
146
147  integer   start_gate,end_gate,d_gate
148
149  logical :: g_to_vol_in_present, g_to_vol_out_present
150
151  ! Logicals to avoid calling present within the loops
152  g_to_vol_in_present  = present(g_to_vol_in)
153  g_to_vol_out_present = present(g_to_vol_out)
154
155  !
156  ! load scaling matricies from disk -- but only the first time this subroutine is called
157  !
158  if(hp%load_scale_LUTs) then
159    call load_scale_LUTs(hp)
160    hp%load_scale_LUTs=.false.
161    hp%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
162  endif
163
164  pi = acos(-1.0)
165
166!   ----- Initialisation -----
167  g_to_vol = 0.0
168  a_to_vol = 0.0
169  z_vol    = 0.0
170  z_ray    = 0.0
171  kr_vol   = 0.0
172
173!   // loop over each range gate (ngate) ... starting with layer closest to the radar !
174  if(hp%radar_at_layer_one) then
175    start_gate=1
176    end_gate=ngate
177    d_gate=1
178  else
179    start_gate=ngate
180    end_gate=1
181    d_gate=-1
182  endif
183  do k=start_gate,end_gate,d_gate
184  ! // loop over each profile (nprof)
185    do pr=1,nprof
186      t_kelvin = t_matrix(pr,k)
187!     :: determine if hydrometeor(s) present in volume
188      hydro = .false.
189      do j=1,hp%nhclass
190        if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
191          hydro = .true.
192          exit
193        endif
194      enddo
195
196!     :: if there is hydrometeor in the volume
197      if (hydro) then
198
199        rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin))
200!       :: loop over hydrometeor type
201        do tp=1,hp%nhclass
202          if (hm_matrix(tp,pr,k) <= 1E-12) cycle
203          phase = hp%phase(tp)
204          if (phase==0) then
205            itt = infind(hp%mt_ttl,t_kelvin)
206          else
207            itt = infind(hp%mt_tti,t_kelvin)
208          endif
209          if (re_matrix(tp,pr,k).eq.0) then
210            call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, &
211              hp%dtype(tp),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
212              hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),Re)
213            re_matrix(tp,pr,k)=Re
214          else
215            if (Np_matrix(tp,pr,k)>0) then
216               print *, 'Warning: Re and Np set for the same ', &
217                        'volume & hydrometeor type.  Np is being ignored.'
218            endif
219            Re = re_matrix(tp,pr,k)
220          endif
221
222          iRe_type=1
223          if(Re.gt.0) then
224            ! determine index in to scale LUT
225            !
226            ! distance between Re points (defined by "base" and "step") for
227            ! each interval of size Re_BIN_LENGTH
228            ! Integer asignment, avoids calling floor intrinsic
229            n=Re/Re_BIN_LENGTH
230            if (n>=Re_MAX_BIN) n=Re_MAX_BIN-1
231            step=hp%step_list(n+1)
232            base=hp%base_list(n+1)
233            iRe_type=Re/step
234            if (iRe_type.lt.1) iRe_type=1
235
236            Re=step*(iRe_type+0.5)      ! set value of Re to closest value allowed in LUT.
237            iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)
238
239            ! make sure iRe_type is within bounds
240            if (iRe_type.ge.nRe_types) then
241!               write(*,*) 'Warning: size of Re exceed value permitted ', &
242!                    'in Look-Up Table (LUT).  Will calculate. '
243               ! no scaling allowed
244               iRe_type=nRe_types
245               hp%Z_scale_flag(tp,itt,iRe_type)=.false.
246            else
247               ! set value in re_matrix to closest values in LUT
248              if (.not. DO_LUT_TEST) re_matrix(tp,pr,k)=Re
249            endif
250          endif
251          ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
252          ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
253          if( (.not. hp%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST)  then
254            ! :: create a distribution of hydrometeors within volume
255            select case(hp%dtype(tp))
256              case(4)
257                ns = 1
258                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
259                Di = hp%p1(tp)
260                Ni = 0.
261              case default
262                ns = nd   ! constant defined in radar_simulator_types.f90
263                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
264                Di = hp%D
265                Ni = 0.
266            end select
267            call dsd(hm_matrix(tp,pr,k),re_matrix(tp,pr,k),Np_matrix(tp,pr,k), &
268                     Di,Ni,ns,hp%dtype(tp),rho_a,t_kelvin, &
269                     hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
270                     hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp))
271
272            ! calculate particle density
273            if (phase == 1) then
274              if (hp%rho(tp) < 0) then
275                ! Use equivalent volume spheres.
276                hp%rho_eff(tp,1:ns,iRe_type) = 917                              ! solid ice == equivalent volume approach
277                Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6
278                ! alternative is to comment out above two lines and use the following block
279                ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density
280                !
281                ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3)   !MG Mie approach
282
283                ! as the particle size gets small it is possible that the mass to size relationship of
284                ! (given by power law in hclass.data) can produce impossible results
285                ! where the mass is larger than a solid sphere of ice. 
286                ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
287                ! do i=1,ns
288                ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
289                ! hp%rho_eff(tp,i,iRe_type) = 917
290                ! endif
291                ! enddo
292              else
293                ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
294                hp%rho_eff(tp,1:ns,iRe_type) = 917
295                Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0))
296                ! alternative ... coment out above two lines and use the following for MG-Mie
297                ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp)   !MG Mie approach
298              endif
299            else
300              ! I assume here that water phase droplets are spheres.
301              ! hp%rho should be ~ 1000  or hp%apm=524 .and. hp%bpm=3
302              Deq = Di
303            endif
304
305            ! calculate effective reflectivity factor of volume
306            xxa = -9.9
307            rhoi = hp%rho_eff(tp,1:ns,iRe_type)
308            call zeff(hp%freq,Deq,Ni,ns,hp%k2,t_kelvin,phase,hp%do_ray, &
309                      ze,zr,kr,xxa,xxa,rhoi)
310
311            ! test code ... compare Np value input to routine with sum of DSD
312            ! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation
313            ! not just the DSD representation given by Ni
314            if(Np_matrix(tp,pr,k)>0 .and. DO_NP_TEST ) then
315              Np = path_integral(Ni,Di,1,ns-1)/rho_a*1E6
316              ! Note: Representation is not great or small Re < 2
317              if( (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)>0.1 ) then
318                write(*,*) 'Error: Np input does not match sum(N)'
319                write(*,*) tp,pr,k,Re,Ni(1),Ni(ns),10*log10(ze)
320                write(*,*) Np_matrix(tp,pr,k),Np,(Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)
321                write(*,*)
322              endif
323            endif
324
325            deallocate(Di,Ni,rhoi,xxa,Deq)
326
327            ! LUT test code
328            ! This segment of code compares full calculation to scaling result
329            if ( hp%Z_scale_flag(tp,itt,iRe_type) .and. DO_LUT_TEST )  then
330              scale_factor=rho_a*hm_matrix(tp,pr,k)
331              ! if more than 2 dBZe difference print error message/parameters.
332              if ( abs(10*log10(ze) - 10*log10(hp%Ze_scaled(tp,itt,iRe_type) * &
333                   scale_factor)) > 2 ) then
334                write(*,*) 'Roj Error: ',tp,itt,iRe_type,hp%Z_scale_flag(tp,itt,iRe_type),n,step,base
335                write(*,*) 10*log10(ze),10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor)
336                write(*,*) hp%Ze_scaled(tp,itt,iRe_type),scale_factor
337                write(*,*) re_matrix(tp,pr,k),Re
338                write(*,*)
339              endif
340            endif
341
342          else ! can use z scaling
343            scale_factor=rho_a*hm_matrix(tp,pr,k)
344            zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor
345            ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
346            kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor
347          endif  ! end z_scaling
348
349          kr_vol(pr,k) = kr_vol(pr,k) + kr
350          z_vol(pr,k)  = z_vol(pr,k)  + ze
351          z_ray(pr,k)  = z_ray(pr,k)  + zr
352
353          ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
354          if ( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then
355            if (iRe_type>1) then
356              scale_factor=rho_a*hm_matrix(tp,pr,k)
357              hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
358              hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
359              hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
360              hp%Z_scale_flag(tp,itt,iRe_type) = .true.
361              hp%Z_scale_added_flag(tp,itt,iRe_type)=.true.
362            endif
363          endif
364
365        enddo   ! end loop of tp (hydrometeor type)
366
367      else
368!     :: volume is hydrometeor-free     
369        kr_vol(pr,k) = 0
370        z_vol(pr,k)  = undef
371        z_ray(pr,k)  = undef
372      endif
373
374      !     :: attenuation due to hydrometeors between radar and volume
375      !
376      ! NOTE old scheme integrates attenuation only for the layers ABOVE
377      ! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
378      ! which may be a problem.   ROJ
379      ! in the new scheme I assign half the attenuation to the current layer
380      if(d_gate==1) then
381        ! dheight calcuations assumes hgt_matrix points are the cell mid-points.
382        if (k>2) then
383          ! add to previous value to half of above layer + half of current layer
384          a_to_vol(pr,k)=  a_to_vol(pr,k-1) + &
385             (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
386        else
387          a_to_vol(pr,k)=  kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
388        endif
389      else   ! d_gate==-1
390        if(k<ngate) then
391          ! add to previous value half of above layer + half of current layer
392          a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
393              (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
394        else
395          a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
396        endif
397      endif
398
399      !     :: attenuation due to gaseous absorption between radar and volume
400      if (g_to_vol_in_present) then
401        g_to_vol(pr,k) = g_to_vol_in(pr,k)
402      else
403        if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .and. (pr == 1)) ) then
404          g_vol(pr,k) = gases(p_matrix(pr,k),t_kelvin,rh_matrix(pr,k),hp%freq)
405          if (d_gate==1) then
406            if (k>1) then
407              ! add to previous value to half of above layer + half of current layer
408              g_to_vol(pr,k) =  g_to_vol(pr,k-1) + &
409                  0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
410            else
411              g_to_vol(pr,k)=  0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
412            endif
413          else   ! d_gate==-1
414            if (k<ngate) then
415              ! add to previous value to half of above layer + half of current layer
416              g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
417                 0.5*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
418            else
419              g_to_vol(pr,k)= 0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
420            endif
421          endif
422        elseif(hp%use_gas_abs == 2) then
423          ! using value calculated for the first column
424          g_to_vol(pr,k) = g_to_vol(1,k)
425        elseif (hp%use_gas_abs == 0) then
426          g_to_vol(pr,k) = 0
427        endif
428      endif
429
430      ! Compute Rayleigh reflectivity, and full, attenuated reflectivity
431      if ((hp%do_ray == 1) .and. (z_ray(pr,k) > 0)) then
432        Ze_ray(pr,k) = 10*log10(z_ray(pr,k))
433      else
434        Ze_ray(pr,k) = undef
435      endif
436      if (z_vol(pr,k) > 0) then
437        Ze_non(pr,k) = 10*log10(z_vol(pr,k))
438        dBZe(pr,k) = Ze_non(pr,k)-a_to_vol(pr,k)-g_to_vol(pr,k)
439      else
440        dBZe(pr,k) = undef
441        Ze_non(pr,k) = undef
442      endif
443
444    enddo   ! end loop over pr (profile)
445
446  enddo ! end loop of k (range gate)
447
448  ! Output array with gaseous absorption
449  if (g_to_vol_out_present) g_to_vol_out = g_to_vol
450
451  ! save any updates made
452  if (hp%update_scale_LUTs) call save_scale_LUTs(hp)
453
454  end subroutine radar_simulator
Note: See TracBrowser for help on using the repository browser.