source: LMDZ6/trunk/libf/phylmd/cosp2/quickbeam_optics.F90 @ 3358

Last change on this file since 3358 was 3358, checked in by idelkadi, 6 years ago

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