#include "cosp_defs.H" #ifdef COSP_GFDL !--------------------------------------------------------------------- !------------ FMS version number and tagname for this file ----------- ! $Id: radar_simulator.F90,v 20.0 2013/12/13 23:16:51 fms Exp $ ! $Name: tikal $ ! cosp_version = 1.3.2 #endif subroutine radar_simulator(freq,k2,do_ray,use_gas_abs,use_mie_table,mt, & nhclass,hp,nprof,ngate,nsizes,D,hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix, & rh_matrix,Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, & g_to_vol_in,g_to_vol_out) ! rh_matrix,Ze_non,Ze_ray,kr_matrix,g_atten_to_vol,dBZe) use m_mrgrnk use array_lib use math_lib use optics_lib use radar_simulator_types implicit none ! Purpose: ! Simulates a vertical profile of radar reflectivity ! Part of QuickBeam v1.04 by John Haynes & Roger Marchand ! ! Inputs: ! [freq] radar frequency (GHz), can be anything unless ! use_mie_table=1, in which case one of 94,35,13.8,9.6,3 ! [k2] |K|^2, the dielectric constant, set to -1 to use the ! frequency dependent default ! [do_ray] 1=do Rayleigh calcs, 0=not ! [use_gas_abs] 1=do gaseous abs calcs, 0=not, ! 2=use same as first profile (undocumented) ! [use_mie_table] 1=use Mie tables, 0=not ! [mt] Mie look up table ! [nhclass] number of hydrometeor types ! [hp] structure that defines hydrometeor types ! [nprof] number of hydrometeor profiles ! [ngate] number of vertical layers ! [nsizes] number of discrete particles in [D] ! [D] array of discrete particles (um) ! ! (The following 5 arrays must be in order from closest to the radar ! to farthest...) ! [hgt_matrix] height of hydrometeors (km) ! [hm_matrix] table of hydrometeor mixing rations (g/kg) ! [re_matrix] OPTIONAL table of hydrometeor effective radii (microns) ! [p_matrix] pressure profile (hPa) ! [t_matrix] temperature profile (C) ! [rh_matrix] relative humidity profile (%) ! ! 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) ! ----- INPUTS ----- type(mie), intent(in) :: mt type(class_param), intent(inout) :: hp #ifdef COSP_GFDL real*8, intent(in) :: freq real*8, intent(inout) :: k2 #else real*8, intent(in) :: freq,k2 #endif integer, intent(in) :: do_ray,use_gas_abs,use_mie_table, & nhclass,nprof,ngate,nsizes real*8, dimension(nsizes), intent(in) :: D real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, & t_matrix,rh_matrix real*8, dimension(nhclass,nprof,ngate), intent(in) :: hm_matrix real*8, dimension(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) xx_rhoi ! to expand parameter xx to array parameter 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, Re, ld, tmp1, apm,bpm integer*4 :: tp, i, j, k, pr, itt, iff real*8 bin_length,step,base,step_list(25),base_list(25) 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) ! set up Re bins for z_scalling bin_length=50; max_bin=25 step_list(1)=1 base_list(1)=75 do j=2,max_bin step_list(j)=3*(j-1); if(step_list(j)>bin_length) then step_list(j)=bin_length; endif base_list(j)=base_list(j-1)+floor(bin_length/step_list(j-1)); enddo pi = acos(-1.0) if (use_mie_table == 1) iff = infind(mt%freq,freq,sort=1) ! // 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,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 (hydro(k) == 1) then ! :: if there is hydrometeor in the volume rho_a = (p_matrix(pr,k)*100.)/(287*(t_matrix(pr,k)+273.15)) ! :: loop over hydrometeor type do tp=1,nhclass if (hm_matrix(tp,pr,k) <= 1E-12) cycle phase = hp%phase(tp) if(phase==0) then itt = infind(mt_ttl,t_matrix(pr,k)) else itt = infind(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 iRe_type=1 Re=re_matrix(tp,pr,k) n=floor(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>max_bin) then n=max_bin endif step=step_list(n) base=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) iRe_type=iRe_type+base-floor(n*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_flag(tp,itt,iRe_type)=.false. hp%scaled(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 #ifdef COSP_GFDL !RSH: ! In order to get reproducibility across physics windows, the "else" branch ! of the following if block must be removed. However doing so increases ! total model time by over a factor of 2. Therefore I have left the loop ! in place for now, pending a response to this issue from the COSP ! developers. Instead, to avoid generating incorrect diagnostics, the ! following diagnostics must be removed from the diag_table: ! "cloudsatcfad_x", where x=1,15 ! "cloudsatcfad_mdl_x", where x=1,15 ! "cloudsatcfad_sat_x", where x=1,15 ! "clcalipso2_sat" ! "clcalipso2_mdl" ! "clcalipso2" ! ! If these diagnostics are required, they may be correctly obtained ! (at great cost) by removing the else block (see !RSH section below) and ! always executing the first part of this if loop. #endif if( .not. hp%z_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),xx_rhoi(ns),xxa(ns),Deq(ns)) if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns)) Di = hp%p1(tp) Ni = 0. case default ns = nsizes allocate(Di(ns),Ni(ns),rhoi(ns),xx_rhoi(ns),xxa(ns),Deq(ns)) if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns)) Di = D Ni = 0. end select xx_rhoi = xx ! :: create a DSD (using scaling factor if applicable) ! hp%scaled(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%scaled(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. 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)) endif ! if using equivalent volume spheres if (use_mie_table == 1) then Ntemp=Ni ! Find N(Di) from N(Deq) which we know do i=1,ns j=infind(Deq,Di(i)) Ni(i)=Ntemp(j) enddo else ! just use Deq and D variable input to mie code Di=Deq; endif endif rhoi = hp%rho_eff(tp,1:ns,iRe_type) ! :: calculate effective reflectivity factor of volume if (use_mie_table == 1) then #ifndef COSP_GFDL if ((hp%dtype(tp) == 4) .and. (hp%idd(tp) < 0)) then hp%idd(tp) = infind(mt%D,Di(1)) endif if (phase == 0) then ! itt = infind(mt_ttl,t_matrix(pr,k)) select case(hp%dtype(tp)) case(4) mt_qext(1) = mt%qext(hp%idd(tp),itt,1,iff) mt_qbsca(1) = mt%qbsca(hp%idd(tp),itt,1,iff) case default mt_qext = mt%qext(:,itt,1,iff) mt_qbsca = mt%qbsca(:,itt,1,iff) end select call zeff(freq,Di,Ni,ns,k2,mt_ttl(itt),0,do_ray, & ze,zr,kr,mt_qext,mt_qbsca,xx_rhoi) else ! itt = infind(mt_tti,t_matrix(pr,k)) select case(hp%dtype(tp)) case(4) if (hp%ifc(tp,1,iRe_type) < 0) then hp%ifc(tp,1,iRe_type) = infind(mt%f,rhoi(1)/917.) endif mt_qext(1) = & mt%qext(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff) mt_qbsca(1) = & mt%qbsca(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff) case default do i=1,ns if (hp%ifc(tp,i,iRe_type) < 0) then hp%ifc(tp,i,iRe_type) = infind(mt%f,rhoi(i)/917.) endif mt_qext(i) = mt%qext(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff) mt_qbsca(i) = mt%qbsca(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff) enddo end select call zeff(freq,Di,Ni,ns,k2,mt_tti(itt),1,do_ray, & ze,zr,kr,mt_qext,mt_qbsca,xx_rhoi) endif #endif else xxa = -9.9 call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, & ze,zr,kr,xxa,xxa,rhoi) endif ! end of use mie table ! xxa = -9.9 !call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, & ! ze2,zr,kr2,xxa,xxa,rhoi) ! if(abs(ze2-ze)/ze2 > 0.1) then ! if(abs(kr2-kr)/kr2 > 0.1) then ! write(*,*) pr,k,tp,ze2,ze2-ze,abs(ze2-ze)/ze2,itt+cnt_liq,iff ! write(*,*) pr,k,tp,ze2,kr2,kr2-kr,abs(kr2-kr)/kr2 ! stop !endif deallocate(Di,Ni,rhoi, xx_rhoi,xxa,Deq) if (use_mie_table == 1) deallocate(mt_qext,mt_qbsca,Ntemp) #ifdef COSP_GFDL !RSH ! the following lines need be removed for reproducible execution. !with these lines retained the diagnostics listed above will not reproduce !across physics_windows. #endif 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 #ifdef COSP_GFDL !RSH end of lines to remove #endif ! 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_flag(tp,itt,iRe_type) .and. 1.eq.1 ) 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_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_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 ( (use_gas_abs == 1) .or. ((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),freq) g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1) elseif (use_gas_abs == 0) then g_to_vol(k) = 0 endif endif ! kr_matrix(pr,:)=kr_vol ! :: 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 ((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) end subroutine radar_simulator