Context Navigation

quickbeam_optics.F90 @ 4744

Last change on this file since 4744 was 3358, checked in by idelkadi, 6 years ago
Implementation de la nouvelle version COSPv2 dans LMDZ. Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 64.1 KB

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

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

Download in other formats:

Original Format