Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

quickbeam_optics.F90 @ 5018

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

Integration of version 2 of the COSP simulator in LMDZ
This line, and those below, will be ignored--

M makegcm
M makelmdz
M makelmdz_fcm
M libf/phylmd/physiq_mod.F90
A libf/phylmd/cospv2
A libf/phylmd/cospv2/mo_rng.F90
A libf/phylmd/cospv2/quickbeam_optics.F90
A libf/phylmd/cospv2/cosp_cloudsat_interface.F90
A libf/phylmd/cospv2/cosp_config.F90
A libf/phylmd/cospv2/lidar_simulator.F90
A libf/phylmd/cospv2/prec_scops.F90
A libf/phylmd/cospv2/mrgrnk.F90
A libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
A libf/phylmd/cospv2/cosp_atlid_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
A libf/phylmd/cospv2/cosp_math_constants.F90
A libf/phylmd/cospv2/MISR_simulator.F90
A libf/phylmd/cospv2/modis_simulator.F90
A libf/phylmd/cospv2/math_lib.F90
A libf/phylmd/cospv2/cosp_grLidar532_interface.F90
A libf/phylmd/cospv2/cosp_errorHandling.F90
A libf/phylmd/cospv2/cosp_stats.F90
A libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
A libf/phylmd/cospv2/cosp_utils.F90
A libf/phylmd/cospv2/cosp_optics.F90
A libf/phylmd/cospv2/icarus.F90
A libf/phylmd/cospv2/scops.F90
A libf/phylmd/cospv2/optics_lib.F90
A libf/phylmd/cospv2/cosp_kinds.F90
A libf/phylmd/cospv2/cosp_calipso_interface.F90
A libf/phylmd/cospv2/quickbeam.F90
A libf/phylmd/cospv2/parasol.F90
A libf/phylmd/cospv2/cosp_phys_constants.F90
A libf/phylmd/cospv2/cosp.F90
A libf/phylmd/cospv2/array_lib.F90
A libf/phylmd/cospv2/cosp_isccp_interface.F90
A libf/phylmd/cospv2/cosp_parasol_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_interface.F90
A libf/phylmd/cospv2/cosp_misr_interface.F90
A libf/phylmd/cospv2/cosp_modis_interface.F90

File size: 64.3 KB

Rev	Line
[3491]	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	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	33	module 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
	63	contains
	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/pisd%apm(tp)/rhoice) * one_third ) * ( (Di1E-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)(Di1E-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(10log10(ze) - 10log10(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_wpNprho_a1E-6))2 ! units of Nt = Nprhoa = #/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 * ( Q1E-3gamma(vu)/(apmlocal_Npgamma(vu+bpm)) )**(1/bpm) ! units = microns
	458	endif
	459	Re = 0.5_wpD0gamma(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 = ((apmgamma(1.+bpm)N0)/(rho_aQ1E-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 = (Q1E-3)/(2._wpbpmapm*local_Np)
	539	tmp2 = exp(0.5_wpbpmbpm(log_sigma_g))exp(0.5_wpbpmbpm*(log_sigma_g))
	540	rg = ((tmp1/tmp2)*(1._wp/bpm))1E6
	541	endif
	542	Re = rgexp(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_wpNprho_a1E-6))2._wp ! units of Nt = Nprhoa = #/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._wpRegamma(vu+2)/gamma(vu+3)
	684	fc = (((D1E-6)(vu-1)exp(-1*D/D0)) / &
	685	(apm((D01E-6)*(vu+bpm))gamma(vu+bpm))) * 1E-12
	686	N = fcrho_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 = (((D1E-6)(vu-1)exp(-1*D/D0)) / &
	691	(apm((D01E-6)*(vu+bpm))gamma(vu+bpm))) * 1E-12
	692	N = fcrho_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 = (Q1E-3)*(1./bpm)
	696	fc = (D1E-6 / (gamma(vu)/(apmlocal_npgamma(vu+bpm)))(1._wp/bpm))*vu
	697	N = ((rho_alocal_npfc(D1E-6)*(-1._wp))/(gamma(vu)tmp1*vu) &
	698	exp(-1._wpfc(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 = (ld1E6)(1.+bpm)/(apmgamma(1+bpm))exp(-1._wp(ld1E6)(D1E-6))1E-12
	714	N = fcrho_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 = ((apmgamma(1.+bpm)N0)*tmp1)(D*1E-6)
	720	N = (N0exp(-1._wpfc(1._wp/(rho_aQ1E-3))tmp1)) 1E-12
	721	elseif (abs(p2+1) > 1E-8) then
	722	! Use default value for lambda
	723	ld = p2
	724	fc = (ld1E6)(1._wp+bpm)/(apmgamma(1+bpm))exp(-1._wp(ld1E6)(D1E-6))1E-12
	725	N = fcrho_a(Q*1E-3)
	726	else
	727	! ld "parameterized" from temperature (carry over from original Quickbeam).
	728	ld = 1220._wp10._wp(-0.0245_wptc)*1E-6
	729	N0 = ((ld1E6)(1._wp+bpm)Q1E-3rho_a)/(apm*gamma(1+bpm))
	730	N = (N0exp(-ldD)) * 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 = 6rho_a(bhp+4)
	785	! tmp2 = pirc(dmax_mm*(bhp+4))(1-(dmin_mm/dmax_mm)**(bhp+4))
	786	! ahp = (Q1E-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(Q1E-3)*(bhp+bpm+1)
	791	tmp2 = apm(dmax_mmbhpdmax(bpm+1)-dmin_mmbhpdmin*(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(D01E-6)*(bpm-3)
	813	fc(1) = (6._wp/(piD0D0D0rho_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 = (bpmlog_sigma_g)(bpm*log_sigma_g)
	828	if(Re.le.0) then
	829	rg = p2
	830	else
	831	!rg = Reexp(-2.5(log_sigma_g*log_sigma_g))
	832	rg =Reexp(-2.5_wp(log_sigma_g**2))
	833
	834	endif
	835
	836	fc = 0.5_wp((1._wp/((2._wprg1E-6)(bpm)apm(2._wppi)*(0.5_wp) &
	837	log_sigma_gD0.5_wp1E-6))exp(-0.5_wp((log(0.5_wpD/rg)/log_sigma_g)*2._wp+tmp2)))1E-12
	838	N = fcrho_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(Q1E-3))/(2._wp*bpmapm*N0)
	851	tmp2 = exp(0.5_wpbpmbpm(log_sigma_g))exp(0.5_wpbpmbpm*(log_sigma_g))
	852	rg = ((tmp1/tmp2)*(1/bpm))1E6
	853
	854	N = 0.5_wp(N0 / ((2._wppi)*(0.5_wp)log_sigma_gD0.5_wp1E-6) &
	855	exp((-0.5_wp(log(0.5_wpD/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+(m0m0)+2f((m0m0)-1))/(2+(m0m0)+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.25sum[qbscapiD^2N(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 = qbscaN0D0*D0
	1001	call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),eta_sum)
	1002	endif
	1003
	1004	eta_mie = eta_sum0.25_wppi
	1005	const = ((wlwlwlwl)/(pipipipipi))(1._wp/k2)
	1006
	1007	z0_eff = const*eta_mie
	1008
	1009	! kr = 0.25crsum[qextpiD^2N(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 = qextN0D0*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_sum0.25_wppi(1000._wpcr)
	1025
	1026	! z_ray = sum[D^6N(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 = N0D0D0D0D0D0D0
	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_eff1E18 ! 10.alog10(z0_eff*1E18)
	1039	z_ray = z0_ray1E18 ! 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 = (RHththththth)/(41.45_wp10*(9.834_wpth-10)) ! kPa
	1172	!e = (RHththththth)/(41.45_wp10*(9.834_wpth-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 = pthth*th
	1180	eth35 = eth*(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)(pth**(0.8_wp-a4(i))+aux1)
	1190	gm2 = gm*gm
	1191	delt = a5(i)pth**a6(i)
	1192	x = aux3*aux3+gm2
	1193	y = aux4*aux4+gm2
	1194	fpp_o2 = (((1._wp/x)+(1._wp/y))(gmaux2) - (deltaux2)((aux3/(x))-(aux4/(x))))
	1195	s_o2 = a1(i)pth3exp(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_wpe)th*(0.8_wp)
	1202	a0 = 3.07E-4_wp
	1203	ap = 1.4_wp(1-1.2_wpf*(1.5_wp)1E-5)*1E-10
	1204	term2 = (2a0(gm0(1+(f/gm0)(f/gm0))(1+(f/60._wp)2))(-1) + appth(2.5_wp))fpth*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)(pth**(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))(gmaux2) ! - (deltaux2)((aux3/(x))-(aux4/(x)))
	1218	s_h2o = b1(i)eth35exp(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 = (bfp+beeththth)feth**(2.5_wp)
	1227
	1228	! Summation and result
	1229	npp = term1 + term2 + term3 + term4
	1230	gases = 0.182_wpfnpp
	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
	1380	end module mod_quickbeam_optics

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format