source: LMDZ6/trunk/libf/phylmd/cospv2/quickbeam_optics.f90 @ 5268

Last change on this file since 5268 was 5268, checked in by abarral, 8 days ago

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