subroutine radar_simulator( & hp, & ! hydrometeor and radar simulator properties nprof,ngate, & ! size of input matricies hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix,rh_matrix, & ! inputs Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, & ! outputs g_to_vol_in,g_to_vol_out) ! optional in/out use m_mrgrnk use array_lib use math_lib use optics_lib use radar_simulator_types use scale_LUTs_io implicit none ! Purpose: ! Simulates a vertical profile of radar reflectivity ! Part of QuickBeam v1.04 by John Haynes & Roger Marchand ! ! Inputs: ! ! [hp] structure that defines hydrometeor types and other radar properties ! ! [nprof] number of hydrometeor profiles ! [ngate] number of vertical layers ! ! (The following 5 arrays must be in order from closest to the radar ! to farthest...) ! ! [hgt_matrix] height of hydrometeors (km) ! [p_matrix] pressure profile (hPa) ! [t_matrix] temperature profile (C) ! [rh_matrix] relative humidity profile (%) -- only needed if gaseous aborption calculated. ! ! [hm_matrix] table of hydrometeor mixing rations (g/kg) ! [re_matrix] OPTIONAL table of hydrometeor effective radii (microns)! ! ! Outputs: ! ! [Ze_non] radar reflectivity without attenuation (dBZ) ! [Ze_ray] Rayleigh reflectivity (dBZ) ! [h_atten_to_vol] attenuation by hydromets, radar to vol (dB) ! [g_atten_to_vol] gaseous atteunation, radar to vol (dB) ! [dBZe] effective radar reflectivity factor (dBZ) ! ! Optional: ! [g_to_vol_in] integrated atten due to gases, r>v (dB). ! If present then is used as gaseous absorption, independently of the ! value in use_gas_abs ! [g_to_vol_out] integrated atten due to gases, r>v (dB). ! If present then gaseous absorption for each profile is returned here. ! ! Created: ! 11/28/2005 John Haynes (haynes@atmos.colostate.edu) ! Modified: ! 09/2006 placed into subroutine form, scaling factors (Roger Marchand,JMH) ! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand) ! 01/2008 'Do while' to determine if hydrometeor(s) present in volume ! changed for vectorization purposes (A. Bodas-Salcedo) ! ! 06/2010 V3.0 ... Modified to load or save scale factors to disk as a LUT. ! ... All hydrometeor and radar simulator properties now included in hp structure ! ... hp structure should be initialized by call to radar_simulator_init prior to calling ! ... this subroutine. ! Also ... Support of Morrison 2-moment style microphyscis added ! ... Changes implement by Roj Marchand following work by Laura Fowler ! ! ! ----- INPUTS ----- type(class_param), intent(inout) :: hp integer, intent(in) :: nprof,ngate real*8, dimension(nprof,ngate), intent(in) :: & hgt_matrix, p_matrix,t_matrix,rh_matrix real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix ! ----- OUTPUTS ----- real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, & g_atten_to_vol,dBZe,h_atten_to_vol ! ----- OPTIONAL ----- real*8, optional, dimension(ngate,nprof) :: & g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input ! the same gaseous absorption in different calls. Optional to allow compatibility ! with original version. A. Bodas April 2008. ! real*8, dimension(nprof,ngate) :: kr_matrix ! ----- INTERNAL ----- integer :: & phase, & ! 0=liquid, 1=ice ns ! number of discrete drop sizes integer*4, dimension(ngate) :: & hydro ! 1=hydrometeor in vol, 0=none real*8 :: & rho_a, & ! air density (kg m^-3) gases ! function: 2-way gas atten (dB/km) real*8, dimension(:), allocatable :: & Di, Deq, & ! discrete drop sizes (um) Ni, Ntemp, & ! discrete concentrations (cm^-3 um^-1) rhoi ! discrete densities (kg m^-3) real*8, dimension(ngate) :: & z_vol, & ! effective reflectivity factor (mm^6/m^3) z_ray, & ! reflectivity factor, Rayleigh only (mm^6/m^3) kr_vol, & ! attenuation coefficient hydro (dB/km) g_vol, & ! attenuation coefficient gases (dB/km) a_to_vol, & ! integrated atten due to hydometeors, r>v (dB) g_to_vol ! integrated atten due to gases, r>v (dB) integer,parameter :: KR8 = selected_real_kind(15,300) real*8, parameter :: xx = -1.0_KR8 real*8, dimension(:), allocatable :: xxa real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2,apm,bpm integer*4 :: tp, i, j, k, pr, itt, iff real*8 step,base integer*4 iRe_type,n,max_bin logical :: g_to_vol_in_present, g_to_vol_out_present ! Logicals to avoid calling present within the loops g_to_vol_in_present = present(g_to_vol_in) g_to_vol_out_present = present(g_to_vol_out) ! ! load scaling matricies from disk -- but only the first time this subroutine is called ! if(hp%load_scale_LUTs) then call load_scale_LUTs(hp) hp%load_scale_LUTs=.false. hp%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run endif pi = acos(-1.0) ! // loop over each profile (nprof) do pr=1,nprof ! ----- calculations for each volume ----- z_vol(:) = 0 z_ray(:) = 0 kr_vol(:) = 0 hydro(:) = 0 ! // loop over eacho range gate (ngate) do k=1,ngate ! :: determine if hydrometeor(s) present in volume hydro(k) = 0 do j=1,hp%nhclass ! Do while changed for vectorization purposes (A. B-S) if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then hydro(k) = 1 exit endif enddo ! :: if there is hydrometeor in the volume if (hydro(k) == 1) then rho_a = (p_matrix(pr,k)*100.)/(287*(t_matrix(pr,k)+273.15)) ! :: loop over hydrometeor type do tp=1,hp%nhclass if (hm_matrix(tp,pr,k) <= 1E-12) cycle phase = hp%phase(tp) if(phase==0) then itt = infind(hp%mt_ttl,t_matrix(pr,k)) else itt = infind(hp%mt_tti,t_matrix(pr,k)) endif ! calculate Re if we have an exponential distribution with fixed No ... precipitation type particle if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8) then apm=hp%apm(tp) bpm=hp%bpm(tp) if ((hp%rho(tp) > 0) .and. (apm < 0)) then apm = (pi/6)*hp%rho(tp) bpm = 3. endif tmp1 = 1./(1.+bpm) ld = ((apm*gamma(1.+bpm)*hp%p1(tp))/(rho_a*hm_matrix(tp,pr,k)*1E-3))**tmp1 Re = 1.5E6/ld re_matrix(tp,pr,k) = Re; endif if(re_matrix(tp,pr,k).eq.0) then iRe_type=1 Re=0 else ! determine index in scale LUT iRe_type=1 Re=re_matrix(tp,pr,k) ! distance between Re points (defined by "base" and "step") in the LUT ! depends on the value of RE n=floor(Re/Re_BIN_LENGTH) if(n==0) then if(Re<25) then step=0.5 base=0 else step=1 base=25 endif else if(n>Re_MAX_BIN) then n=Re_MAX_BIN endif step=hp%step_list(n) base=hp%base_list(n) endif iRe_type=floor(Re/step) if(iRe_type.lt.1) then iRe_type=1 endif Re=step*(iRe_type+0.5) ! set value of Re to closest value in allowed in LUT. iRe_type=iRe_type+base-floor(n*Re_BIN_LENGTH/step) ! make sure iRe_type is within bounds if(iRe_type.ge.nRe_types) then ! print *, tp, re_matrix(tp,pr,k), Re, iRe_type ! no scaling allowed Re=re_matrix(tp,pr,k) iRe_type=nRe_types hp%Z_scale_flag(tp,itt,iRe_type)=.false. hp%N_scale_flag(tp,iRe_type)=.false. endif endif ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them ! if not we will calculate Ze, Zr, and Kr from the distribution parameters if( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then ! :: create a distribution of hydrometeors within volume select case(hp%dtype(tp)) case(4) ns = 1 allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) Di = hp%p1(tp) Ni = 0. case default ns = nd ! constant defined in radar_simulator_types.f90 allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) Di = hp%D Ni = 0. end select ! :: create a DSD (using scaling factor if applicable) ! hp%N_scale_flag(tp,iRe_type)=.false. ! turn off N scaling call dsd(hm_matrix(tp,pr,k),Re,Di,Ni,ns,hp%dtype(tp),rho_a, & t_matrix(pr,k),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), & hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),hp%fc(tp,1:ns,iRe_type), & hp%N_scale_flag(tp,iRe_type)) ! :: calculate particle density ! if ((hp%rho_eff(tp,1,iRe_type) < 0) .and. (phase == 1)) then if (phase == 1) then if (hp%rho(tp) < 0) then ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3) !MG Mie approach ! as the particle size gets small it is possible that the mass to size relationship of ! (given by power law in hclass.data) can produce impossible results ! where the mass is larger than a solid sphere of ice. ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere. ! do i=1,ns ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then ! hp%rho_eff(tp,i,iRe_type) = 917 !endif !enddo ! alternative is to use equivalent volume spheres. hp%rho_eff(tp,1:ns,iRe_type) = 917 ! solid ice == equivalent volume approach Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * & ( (Di*1E-6) ** (hp%bpm(tp)/3.0) ) * 1E6 ! Di now really Deq in microns. Di=Deq; else ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp) !MG Mie approach ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3). hp%rho_eff(tp,1:ns,iRe_type) = 917 Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0)) Di=Deq; endif endif rhoi = hp%rho_eff(tp,1:ns,iRe_type) ! :: calculate effective reflectivity factor of volume xxa = -9.9 call zeff(hp%freq,Di,Ni,ns,hp%k2,t_matrix(pr,k),phase,hp%do_ray, & ze,zr,kr,xxa,xxa,rhoi) deallocate(Di,Ni,rhoi,xxa,Deq) else ! can use z scaling if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then ze = hp%Ze_scaled(tp,itt,iRe_type) zr = hp%Zr_scaled(tp,itt,iRe_type) kr = hp%kr_scaled(tp,itt,iRe_type) else scale_factor=rho_a*hm_matrix(tp,pr,k) zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor endif endif ! end z_scaling ! kr=0 kr_vol(k) = kr_vol(k) + kr z_vol(k) = z_vol(k) + ze z_ray(k) = z_ray(k) + zr ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can if( .not. hp%Z_scale_flag(tp,itt,iRe_type)) then if( ( (hp%dtype(tp)==1 .or. hp%dtype(tp)==5 .or. hp%dtype(tp)==2) .and. abs(hp%p1(tp)+1) < 1E-8 ) .or. & ( hp%dtype(tp)==3 .or. hp%dtype(tp)==4 ) & ) then scale_factor=rho_a*hm_matrix(tp,pr,k) hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor hp%Z_scale_flag(tp,itt,iRe_type)=.True. hp%Z_scale_added_flag(tp,itt,iRe_type)=.True. elseif( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then hp%Ze_scaled(tp,itt,iRe_type) = ze hp%Zr_scaled(tp,itt,iRe_type) = zr hp%kr_scaled(tp,itt,iRe_type) = kr hp%Z_scale_flag(tp,itt,iRe_type)=.True. hp%Z_scale_added_flag(tp,itt,iRe_type)=.True. endif endif enddo ! end loop of tp (hydrometeor type) else ! :: volume is hydrometeor-free kr_vol(k) = 0 z_vol(k) = -999 z_ray(k) = -999 endif ! :: attenuation due to hydrometeors between radar and volume a_to_vol(k) = 2*path_integral(kr_vol,hgt_matrix(pr,:),1,k-1) ! :: attenuation due to gaseous absorption between radar and volume if (g_to_vol_in_present) then g_to_vol(k) = g_to_vol_in(k,pr) else if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .and. (pr == 1)) ) then g_vol(k) = gases(p_matrix(pr,k),t_matrix(pr,k)+273.15, & rh_matrix(pr,k),hp%freq) g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1) elseif (hp%use_gas_abs == 0) then g_to_vol(k) = 0 endif endif ! :: store results in matrix for return to calling program h_atten_to_vol(pr,k)=a_to_vol(k) g_atten_to_vol(pr,k)=g_to_vol(k) if ((hp%do_ray == 1) .and. (z_ray(k) > 0)) then Ze_ray(pr,k) = 10*log10(z_ray(k)) else Ze_ray(pr,k) = -999 endif if (z_vol(k) > 0) then dBZe(pr,k) = 10*log10(z_vol(k))-a_to_vol(k)-g_to_vol(k) Ze_non(pr,k) = 10*log10(z_vol(k)) else dBZe(pr,k) = -999 Ze_non(pr,k) = -999 endif enddo ! end loop of k (range gate) ! Output array with gaseous absorption if (g_to_vol_out_present) g_to_vol_out(:,pr) = g_to_vol enddo ! end loop over pr (profile) ! save any updates made if(hp%update_scale_LUTs) then call save_scale_LUTs(hp) endif end subroutine radar_simulator