| 1 | subroutine radar_simulator(freq,k2,do_ray,use_gas_abs,use_mie_table,mt, & |
|---|
| 2 | nhclass,hp,nprof,ngate,nsizes,D,hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix, & |
|---|
| 3 | rh_matrix,Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, & |
|---|
| 4 | g_to_vol_in,g_to_vol_out) |
|---|
| 5 | |
|---|
| 6 | ! rh_matrix,Ze_non,Ze_ray,kr_matrix,g_atten_to_vol,dBZe) |
|---|
| 7 | |
|---|
| 8 | use m_mrgrnk |
|---|
| 9 | use array_lib |
|---|
| 10 | use math_lib |
|---|
| 11 | use optics_lib |
|---|
| 12 | use radar_simulator_types |
|---|
| 13 | implicit none |
|---|
| 14 | |
|---|
| 15 | ! Purpose: |
|---|
| 16 | ! Simulates a vertical profile of radar reflectivity |
|---|
| 17 | ! Part of QuickBeam v1.04 by John Haynes & Roger Marchand |
|---|
| 18 | ! |
|---|
| 19 | ! Inputs: |
|---|
| 20 | ! [freq] radar frequency (GHz), can be anything unless |
|---|
| 21 | ! use_mie_table=1, in which case one of 94,35,13.8,9.6,3 |
|---|
| 22 | ! [k2] |K|^2, the dielectric constant, set to -1 to use the |
|---|
| 23 | ! frequency dependent default |
|---|
| 24 | ! [do_ray] 1=do Rayleigh calcs, 0=not |
|---|
| 25 | ! [use_gas_abs] 1=do gaseous abs calcs, 0=not, |
|---|
| 26 | ! 2=use same as first profile (undocumented) |
|---|
| 27 | ! [use_mie_table] 1=use Mie tables, 0=not |
|---|
| 28 | ! [mt] Mie look up table |
|---|
| 29 | ! [nhclass] number of hydrometeor types |
|---|
| 30 | ! [hp] structure that defines hydrometeor types |
|---|
| 31 | ! [nprof] number of hydrometeor profiles |
|---|
| 32 | ! [ngate] number of vertical layers |
|---|
| 33 | ! [nsizes] number of discrete particles in [D] |
|---|
| 34 | ! [D] array of discrete particles (um) |
|---|
| 35 | ! |
|---|
| 36 | ! (The following 5 arrays must be in order from closest to the radar |
|---|
| 37 | ! to farthest...) |
|---|
| 38 | ! [hgt_matrix] height of hydrometeors (km) |
|---|
| 39 | ! [hm_matrix] table of hydrometeor mixing rations (g/kg) |
|---|
| 40 | ! [re_matrix] OPTIONAL table of hydrometeor effective radii (microns) |
|---|
| 41 | ! [p_matrix] pressure profile (hPa) |
|---|
| 42 | ! [t_matrix] temperature profile (C) |
|---|
| 43 | ! [rh_matrix] relative humidity profile (%) |
|---|
| 44 | ! |
|---|
| 45 | ! Outputs: |
|---|
| 46 | ! [Ze_non] radar reflectivity without attenuation (dBZ) |
|---|
| 47 | ! [Ze_ray] Rayleigh reflectivity (dBZ) |
|---|
| 48 | ! [h_atten_to_vol] attenuation by hydromets, radar to vol (dB) |
|---|
| 49 | ! [g_atten_to_vol] gaseous atteunation, radar to vol (dB) |
|---|
| 50 | ! [dBZe] effective radar reflectivity factor (dBZ) |
|---|
| 51 | ! |
|---|
| 52 | ! Optional: |
|---|
| 53 | ! [g_to_vol_in] integrated atten due to gases, r>v (dB). |
|---|
| 54 | ! If present then is used as gaseous absorption, independently of the |
|---|
| 55 | ! value in use_gas_abs |
|---|
| 56 | ! [g_to_vol_out] integrated atten due to gases, r>v (dB). |
|---|
| 57 | ! If present then gaseous absorption for each profile is returned here. |
|---|
| 58 | ! |
|---|
| 59 | ! Created: |
|---|
| 60 | ! 11/28/2005 John Haynes (haynes@atmos.colostate.edu) |
|---|
| 61 | ! Modified: |
|---|
| 62 | ! 09/2006 placed into subroutine form, scaling factors (Roger Marchand,JMH) |
|---|
| 63 | ! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand) |
|---|
| 64 | ! 01/2008 'Do while' to determine if hydrometeor(s) present in volume |
|---|
| 65 | ! changed for vectorization purposes (A. Bodas-Salcedo) |
|---|
| 66 | |
|---|
| 67 | ! ----- INPUTS ----- |
|---|
| 68 | type(mie), intent(in) :: mt |
|---|
| 69 | type(class_param), intent(inout) :: hp |
|---|
| 70 | real*8, intent(in) :: freq,k2 |
|---|
| 71 | integer, intent(in) :: do_ray,use_gas_abs,use_mie_table, & |
|---|
| 72 | nhclass,nprof,ngate,nsizes |
|---|
| 73 | real*8, dimension(nsizes), intent(in) :: D |
|---|
| 74 | real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, & |
|---|
| 75 | t_matrix,rh_matrix |
|---|
| 76 | real*8, dimension(nhclass,nprof,ngate), intent(in) :: hm_matrix |
|---|
| 77 | real*8, dimension(nhclass,nprof,ngate), intent(inout) :: re_matrix |
|---|
| 78 | |
|---|
| 79 | ! ----- OUTPUTS ----- |
|---|
| 80 | real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, & |
|---|
| 81 | g_atten_to_vol,dBZe,h_atten_to_vol |
|---|
| 82 | |
|---|
| 83 | ! ----- OPTIONAL ----- |
|---|
| 84 | real*8, optional, dimension(ngate,nprof) :: & |
|---|
| 85 | g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input |
|---|
| 86 | ! the same gaseous absorption in different calls. Optional to allow compatibility |
|---|
| 87 | ! with original version. A. Bodas April 2008. |
|---|
| 88 | |
|---|
| 89 | ! real*8, dimension(nprof,ngate) :: kr_matrix |
|---|
| 90 | |
|---|
| 91 | ! ----- INTERNAL ----- |
|---|
| 92 | integer :: & |
|---|
| 93 | phase, & ! 0=liquid, 1=ice |
|---|
| 94 | ns ! number of discrete drop sizes |
|---|
| 95 | |
|---|
| 96 | integer*4, dimension(ngate) :: & |
|---|
| 97 | hydro ! 1=hydrometeor in vol, 0=none |
|---|
| 98 | real*8 :: & |
|---|
| 99 | rho_a, & ! air density (kg m^-3) |
|---|
| 100 | gases ! function: 2-way gas atten (dB/km) |
|---|
| 101 | |
|---|
| 102 | real*8, dimension(:), allocatable :: & |
|---|
| 103 | Di, Deq, & ! discrete drop sizes (um) |
|---|
| 104 | Ni, Ntemp, & ! discrete concentrations (cm^-3 um^-1) |
|---|
| 105 | rhoi ! discrete densities (kg m^-3) |
|---|
| 106 | |
|---|
| 107 | real*8, dimension(ngate) :: & |
|---|
| 108 | z_vol, & ! effective reflectivity factor (mm^6/m^3) |
|---|
| 109 | z_ray, & ! reflectivity factor, Rayleigh only (mm^6/m^3) |
|---|
| 110 | kr_vol, & ! attenuation coefficient hydro (dB/km) |
|---|
| 111 | g_vol, & ! attenuation coefficient gases (dB/km) |
|---|
| 112 | a_to_vol, & ! integrated atten due to hydometeors, r>v (dB) |
|---|
| 113 | g_to_vol ! integrated atten due to gases, r>v (dB) |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | integer,parameter :: KR8 = selected_real_kind(15,300) |
|---|
| 117 | real*8, parameter :: xx = -1.0_KR8 |
|---|
| 118 | real*8, dimension(:), allocatable :: xxa |
|---|
| 119 | real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2,apm,bpm |
|---|
| 120 | integer*4 :: tp, i, j, k, pr, itt, iff |
|---|
| 121 | |
|---|
| 122 | real*8 bin_length,step,base,step_list(25),base_list(25) |
|---|
| 123 | integer*4 iRe_type,n,max_bin |
|---|
| 124 | |
|---|
| 125 | logical :: g_to_vol_in_present, g_to_vol_out_present |
|---|
| 126 | |
|---|
| 127 | ! Logicals to avoid calling present within the loops |
|---|
| 128 | g_to_vol_in_present = present(g_to_vol_in) |
|---|
| 129 | g_to_vol_out_present = present(g_to_vol_out) |
|---|
| 130 | |
|---|
| 131 | ! set up Re bins for z_scalling |
|---|
| 132 | bin_length=50; |
|---|
| 133 | max_bin=25 |
|---|
| 134 | |
|---|
| 135 | step_list(1)=1 |
|---|
| 136 | base_list(1)=75 |
|---|
| 137 | do j=2,max_bin |
|---|
| 138 | step_list(j)=3*(j-1); |
|---|
| 139 | if(step_list(j)>bin_length) then |
|---|
| 140 | step_list(j)=bin_length; |
|---|
| 141 | endif |
|---|
| 142 | base_list(j)=base_list(j-1)+floor(bin_length/step_list(j-1)); |
|---|
| 143 | enddo |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | pi = acos(-1.0) |
|---|
| 147 | if (use_mie_table == 1) iff = infind(mt%freq,freq,sort=1) |
|---|
| 148 | |
|---|
| 149 | |
|---|
| 150 | ! // loop over each profile (nprof) |
|---|
| 151 | do pr=1,nprof |
|---|
| 152 | |
|---|
| 153 | ! ----- calculations for each volume ----- |
|---|
| 154 | z_vol(:) = 0 |
|---|
| 155 | z_ray(:) = 0 |
|---|
| 156 | kr_vol(:) = 0 |
|---|
| 157 | hydro(:) = 0 |
|---|
| 158 | |
|---|
| 159 | ! // loop over eacho range gate (ngate) |
|---|
| 160 | do k=1,ngate |
|---|
| 161 | |
|---|
| 162 | ! :: determine if hydrometeor(s) present in volume |
|---|
| 163 | hydro(k) = 0 |
|---|
| 164 | do j=1,nhclass ! Do while changed for vectorization purposes (A. B-S) |
|---|
| 165 | if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then |
|---|
| 166 | hydro(k) = 1 |
|---|
| 167 | exit |
|---|
| 168 | endif |
|---|
| 169 | enddo |
|---|
| 170 | |
|---|
| 171 | if (hydro(k) == 1) then |
|---|
| 172 | ! :: if there is hydrometeor in the volume |
|---|
| 173 | |
|---|
| 174 | rho_a = (p_matrix(pr,k)*100.)/(287*(t_matrix(pr,k)+273.15)) |
|---|
| 175 | |
|---|
| 176 | ! :: loop over hydrometeor type |
|---|
| 177 | do tp=1,nhclass |
|---|
| 178 | |
|---|
| 179 | if (hm_matrix(tp,pr,k) <= 1E-12) cycle |
|---|
| 180 | |
|---|
| 181 | phase = hp%phase(tp) |
|---|
| 182 | if(phase==0) then |
|---|
| 183 | itt = infind(mt_ttl,t_matrix(pr,k)) |
|---|
| 184 | else |
|---|
| 185 | itt = infind(mt_tti,t_matrix(pr,k)) |
|---|
| 186 | endif |
|---|
| 187 | |
|---|
| 188 | ! calculate Re if we have an exponential distribution with fixed No ... precipitation type particle |
|---|
| 189 | if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8) then |
|---|
| 190 | |
|---|
| 191 | apm=hp%apm(tp) |
|---|
| 192 | bpm=hp%bpm(tp) |
|---|
| 193 | |
|---|
| 194 | if ((hp%rho(tp) > 0) .and. (apm < 0)) then |
|---|
| 195 | apm = (pi/6)*hp%rho(tp) |
|---|
| 196 | bpm = 3. |
|---|
| 197 | endif |
|---|
| 198 | |
|---|
| 199 | tmp1 = 1./(1.+bpm) |
|---|
| 200 | ld = ((apm*gamma(1.+bpm)*hp%p1(tp))/(rho_a*hm_matrix(tp,pr,k)*1E-3))**tmp1 |
|---|
| 201 | |
|---|
| 202 | Re = 1.5E6/ld |
|---|
| 203 | |
|---|
| 204 | re_matrix(tp,pr,k) = Re; |
|---|
| 205 | |
|---|
| 206 | endif |
|---|
| 207 | |
|---|
| 208 | if(re_matrix(tp,pr,k).eq.0) then |
|---|
| 209 | |
|---|
| 210 | iRe_type=1 |
|---|
| 211 | Re=0 |
|---|
| 212 | else |
|---|
| 213 | iRe_type=1 |
|---|
| 214 | Re=re_matrix(tp,pr,k) |
|---|
| 215 | |
|---|
| 216 | n=floor(Re/bin_length) |
|---|
| 217 | if(n==0) then |
|---|
| 218 | if(Re<25) then |
|---|
| 219 | step=0.5 |
|---|
| 220 | base=0 |
|---|
| 221 | else |
|---|
| 222 | step=1 |
|---|
| 223 | base=25 |
|---|
| 224 | endif |
|---|
| 225 | else |
|---|
| 226 | if(n>max_bin) then |
|---|
| 227 | n=max_bin |
|---|
| 228 | endif |
|---|
| 229 | |
|---|
| 230 | step=step_list(n) |
|---|
| 231 | base=base_list(n) |
|---|
| 232 | endif |
|---|
| 233 | |
|---|
| 234 | iRe_type=floor(Re/step) |
|---|
| 235 | |
|---|
| 236 | if(iRe_type.lt.1) then |
|---|
| 237 | iRe_type=1 |
|---|
| 238 | endif |
|---|
| 239 | |
|---|
| 240 | Re=step*(iRe_type+0.5) |
|---|
| 241 | iRe_type=iRe_type+base-floor(n*bin_length/step) |
|---|
| 242 | |
|---|
| 243 | ! make sure iRe_type is within bounds |
|---|
| 244 | if(iRe_type.ge.nRe_types) then |
|---|
| 245 | |
|---|
| 246 | ! print *, tp, re_matrix(tp,pr,k), Re, iRe_type |
|---|
| 247 | |
|---|
| 248 | ! no scaling allowed |
|---|
| 249 | Re=re_matrix(tp,pr,k) |
|---|
| 250 | |
|---|
| 251 | iRe_type=nRe_types |
|---|
| 252 | hp%z_flag(tp,itt,iRe_type)=.false. |
|---|
| 253 | hp%scaled(tp,iRe_type)=.false. |
|---|
| 254 | endif |
|---|
| 255 | endif |
|---|
| 256 | |
|---|
| 257 | ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them |
|---|
| 258 | ! if not we will calculate Ze, Zr, and Kr from the distribution parameters |
|---|
| 259 | if( .not. hp%z_flag(tp,itt,iRe_type) ) then |
|---|
| 260 | |
|---|
| 261 | ! :: create a distribution of hydrometeors within volume |
|---|
| 262 | select case(hp%dtype(tp)) |
|---|
| 263 | case(4) |
|---|
| 264 | ns = 1 |
|---|
| 265 | allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) |
|---|
| 266 | if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns)) |
|---|
| 267 | Di = hp%p1(tp) |
|---|
| 268 | Ni = 0. |
|---|
| 269 | case default |
|---|
| 270 | ns = nsizes |
|---|
| 271 | allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) |
|---|
| 272 | if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns)) |
|---|
| 273 | Di = D |
|---|
| 274 | Ni = 0. |
|---|
| 275 | end select |
|---|
| 276 | |
|---|
| 277 | ! :: create a DSD (using scaling factor if applicable) |
|---|
| 278 | ! hp%scaled(tp,iRe_type)=.false. ! turn off N scaling |
|---|
| 279 | |
|---|
| 280 | call dsd(hm_matrix(tp,pr,k),Re,Di,Ni,ns,hp%dtype(tp),rho_a, & |
|---|
| 281 | t_matrix(pr,k),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), & |
|---|
| 282 | hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),hp%fc(tp,1:ns,iRe_type), & |
|---|
| 283 | hp%scaled(tp,iRe_type)) |
|---|
| 284 | |
|---|
| 285 | ! :: calculate particle density |
|---|
| 286 | ! if ((hp%rho_eff(tp,1,iRe_type) < 0) .and. (phase == 1)) then |
|---|
| 287 | if (phase == 1) then |
|---|
| 288 | if (hp%rho(tp) < 0) then |
|---|
| 289 | |
|---|
| 290 | ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density |
|---|
| 291 | ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3) !MG Mie approach |
|---|
| 292 | |
|---|
| 293 | ! as the particle size gets small it is possible that the mass to size relationship of |
|---|
| 294 | ! (given by power law in hclass.data) can produce impossible results |
|---|
| 295 | ! where the mass is larger than a solid sphere of ice. |
|---|
| 296 | ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere. |
|---|
| 297 | ! do i=1,ns |
|---|
| 298 | ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then |
|---|
| 299 | ! hp%rho_eff(tp,i,iRe_type) = 917 |
|---|
| 300 | !endif |
|---|
| 301 | !enddo |
|---|
| 302 | |
|---|
| 303 | ! alternative is to use equivalent volume spheres. |
|---|
| 304 | hp%rho_eff(tp,1:ns,iRe_type) = 917 ! solid ice == equivalent volume approach |
|---|
| 305 | Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * & |
|---|
| 306 | ( (Di*1E-6) ** (hp%bpm(tp)/3.0) ) * 1E6 ! Di now really Deq in microns. |
|---|
| 307 | |
|---|
| 308 | else |
|---|
| 309 | |
|---|
| 310 | ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp) !MG Mie approach |
|---|
| 311 | |
|---|
| 312 | ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3). |
|---|
| 313 | hp%rho_eff(tp,1:ns,iRe_type) = 917 |
|---|
| 314 | Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0)) |
|---|
| 315 | |
|---|
| 316 | endif |
|---|
| 317 | |
|---|
| 318 | ! if using equivalent volume spheres |
|---|
| 319 | if (use_mie_table == 1) then |
|---|
| 320 | |
|---|
| 321 | Ntemp=Ni |
|---|
| 322 | |
|---|
| 323 | ! Find N(Di) from N(Deq) which we know |
|---|
| 324 | do i=1,ns |
|---|
| 325 | j=infind(Deq,Di(i)) |
|---|
| 326 | Ni(i)=Ntemp(j) |
|---|
| 327 | enddo |
|---|
| 328 | else |
|---|
| 329 | ! just use Deq and D variable input to mie code |
|---|
| 330 | Di=Deq; |
|---|
| 331 | endif |
|---|
| 332 | |
|---|
| 333 | endif |
|---|
| 334 | rhoi = hp%rho_eff(tp,1:ns,iRe_type) |
|---|
| 335 | |
|---|
| 336 | ! :: calculate effective reflectivity factor of volume |
|---|
| 337 | if (use_mie_table == 1) then |
|---|
| 338 | |
|---|
| 339 | if ((hp%dtype(tp) == 4) .and. (hp%idd(tp) < 0)) then |
|---|
| 340 | hp%idd(tp) = infind(mt%D,Di(1)) |
|---|
| 341 | endif |
|---|
| 342 | |
|---|
| 343 | if (phase == 0) then |
|---|
| 344 | |
|---|
| 345 | ! itt = infind(mt_ttl,t_matrix(pr,k)) |
|---|
| 346 | select case(hp%dtype(tp)) |
|---|
| 347 | case(4) |
|---|
| 348 | mt_qext(1) = mt%qext(hp%idd(tp),itt,1,iff) |
|---|
| 349 | mt_qbsca(1) = mt%qbsca(hp%idd(tp),itt,1,iff) |
|---|
| 350 | case default |
|---|
| 351 | mt_qext = mt%qext(:,itt,1,iff) |
|---|
| 352 | mt_qbsca = mt%qbsca(:,itt,1,iff) |
|---|
| 353 | end select |
|---|
| 354 | |
|---|
| 355 | call zeff(freq,Di,Ni,ns,k2,mt_ttl(itt),0,do_ray, & |
|---|
| 356 | ze,zr,kr,mt_qext,mt_qbsca,xx) |
|---|
| 357 | |
|---|
| 358 | else |
|---|
| 359 | |
|---|
| 360 | ! itt = infind(mt_tti,t_matrix(pr,k)) |
|---|
| 361 | select case(hp%dtype(tp)) |
|---|
| 362 | case(4) |
|---|
| 363 | if (hp%ifc(tp,1,iRe_type) < 0) then |
|---|
| 364 | hp%ifc(tp,1,iRe_type) = infind(mt%f,rhoi(1)/917.) |
|---|
| 365 | endif |
|---|
| 366 | mt_qext(1) = & |
|---|
| 367 | mt%qext(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff) |
|---|
| 368 | mt_qbsca(1) = & |
|---|
| 369 | mt%qbsca(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff) |
|---|
| 370 | case default |
|---|
| 371 | do i=1,ns |
|---|
| 372 | if (hp%ifc(tp,i,iRe_type) < 0) then |
|---|
| 373 | hp%ifc(tp,i,iRe_type) = infind(mt%f,rhoi(i)/917.) |
|---|
| 374 | endif |
|---|
| 375 | mt_qext(i) = mt%qext(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff) |
|---|
| 376 | mt_qbsca(i) = mt%qbsca(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff) |
|---|
| 377 | enddo |
|---|
| 378 | end select |
|---|
| 379 | |
|---|
| 380 | call zeff(freq,Di,Ni,ns,k2,mt_tti(itt),1,do_ray, & |
|---|
| 381 | ze,zr,kr,mt_qext,mt_qbsca,xx) |
|---|
| 382 | |
|---|
| 383 | endif |
|---|
| 384 | |
|---|
| 385 | else |
|---|
| 386 | |
|---|
| 387 | xxa = -9.9 |
|---|
| 388 | call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, & |
|---|
| 389 | ze,zr,kr,xxa,xxa,rhoi) |
|---|
| 390 | |
|---|
| 391 | |
|---|
| 392 | endif ! end of use mie table |
|---|
| 393 | |
|---|
| 394 | ! xxa = -9.9 |
|---|
| 395 | !call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, & |
|---|
| 396 | ! ze2,zr,kr2,xxa,xxa,rhoi) |
|---|
| 397 | |
|---|
| 398 | ! if(abs(ze2-ze)/ze2 > 0.1) then |
|---|
| 399 | ! if(abs(kr2-kr)/kr2 > 0.1) then |
|---|
| 400 | |
|---|
| 401 | ! write(*,*) pr,k,tp,ze2,ze2-ze,abs(ze2-ze)/ze2,itt+cnt_liq,iff |
|---|
| 402 | ! write(*,*) pr,k,tp,ze2,kr2,kr2-kr,abs(kr2-kr)/kr2 |
|---|
| 403 | ! stop |
|---|
| 404 | |
|---|
| 405 | !endif |
|---|
| 406 | |
|---|
| 407 | deallocate(Di,Ni,rhoi,xxa,Deq) |
|---|
| 408 | if (use_mie_table == 1) deallocate(mt_qext,mt_qbsca,Ntemp) |
|---|
| 409 | |
|---|
| 410 | else ! can use z scaling |
|---|
| 411 | |
|---|
| 412 | if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then |
|---|
| 413 | |
|---|
| 414 | ze = hp%Ze_scaled(tp,itt,iRe_type) |
|---|
| 415 | zr = hp%Zr_scaled(tp,itt,iRe_type) |
|---|
| 416 | kr = hp%kr_scaled(tp,itt,iRe_type) |
|---|
| 417 | |
|---|
| 418 | else |
|---|
| 419 | scale_factor=rho_a*hm_matrix(tp,pr,k) |
|---|
| 420 | |
|---|
| 421 | zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 422 | ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 423 | kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 424 | endif |
|---|
| 425 | |
|---|
| 426 | endif ! end z_scaling |
|---|
| 427 | |
|---|
| 428 | ! kr=0 |
|---|
| 429 | |
|---|
| 430 | kr_vol(k) = kr_vol(k) + kr |
|---|
| 431 | z_vol(k) = z_vol(k) + ze |
|---|
| 432 | z_ray(k) = z_ray(k) + zr |
|---|
| 433 | |
|---|
| 434 | ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can |
|---|
| 435 | if( .not. hp%z_flag(tp,itt,iRe_type) .and. 1.eq.1 ) then |
|---|
| 436 | |
|---|
| 437 | if( ( (hp%dtype(tp)==1 .or. hp%dtype(tp)==5 .or. hp%dtype(tp)==2) .and. abs(hp%p1(tp)+1) < 1E-8 ) .or. & |
|---|
| 438 | ( hp%dtype(tp)==3 .or. hp%dtype(tp)==4 ) & |
|---|
| 439 | ) then |
|---|
| 440 | |
|---|
| 441 | scale_factor=rho_a*hm_matrix(tp,pr,k) |
|---|
| 442 | |
|---|
| 443 | hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor |
|---|
| 444 | hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor |
|---|
| 445 | hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor |
|---|
| 446 | |
|---|
| 447 | hp%z_flag(tp,itt,iRe_type)=.True. |
|---|
| 448 | |
|---|
| 449 | elseif( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then |
|---|
| 450 | |
|---|
| 451 | hp%Ze_scaled(tp,itt,iRe_type) = ze |
|---|
| 452 | hp%Zr_scaled(tp,itt,iRe_type) = zr |
|---|
| 453 | hp%kr_scaled(tp,itt,iRe_type) = kr |
|---|
| 454 | |
|---|
| 455 | hp%z_flag(tp,itt,iRe_type)=.True. |
|---|
| 456 | endif |
|---|
| 457 | |
|---|
| 458 | endif |
|---|
| 459 | |
|---|
| 460 | enddo ! end loop of tp (hydrometeor type) |
|---|
| 461 | |
|---|
| 462 | else |
|---|
| 463 | ! :: volume is hydrometeor-free |
|---|
| 464 | |
|---|
| 465 | kr_vol(k) = 0 |
|---|
| 466 | z_vol(k) = -999 |
|---|
| 467 | z_ray(k) = -999 |
|---|
| 468 | |
|---|
| 469 | endif |
|---|
| 470 | |
|---|
| 471 | ! :: attenuation due to hydrometeors between radar and volume |
|---|
| 472 | a_to_vol(k) = 2*path_integral(kr_vol,hgt_matrix(pr,:),1,k-1) |
|---|
| 473 | |
|---|
| 474 | ! :: attenuation due to gaseous absorption between radar and volume |
|---|
| 475 | if (g_to_vol_in_present) then |
|---|
| 476 | g_to_vol(k) = g_to_vol_in(k,pr) |
|---|
| 477 | else |
|---|
| 478 | if ( (use_gas_abs == 1) .or. ((use_gas_abs == 2) .and. (pr == 1)) ) then |
|---|
| 479 | g_vol(k) = gases(p_matrix(pr,k),t_matrix(pr,k)+273.15, & |
|---|
| 480 | rh_matrix(pr,k),freq) |
|---|
| 481 | g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1) |
|---|
| 482 | elseif (use_gas_abs == 0) then |
|---|
| 483 | g_to_vol(k) = 0 |
|---|
| 484 | endif |
|---|
| 485 | endif |
|---|
| 486 | |
|---|
| 487 | ! kr_matrix(pr,:)=kr_vol |
|---|
| 488 | |
|---|
| 489 | ! :: store results in matrix for return to calling program |
|---|
| 490 | h_atten_to_vol(pr,k)=a_to_vol(k) |
|---|
| 491 | g_atten_to_vol(pr,k)=g_to_vol(k) |
|---|
| 492 | if ((do_ray == 1) .and. (z_ray(k) > 0)) then |
|---|
| 493 | Ze_ray(pr,k) = 10*log10(z_ray(k)) |
|---|
| 494 | else |
|---|
| 495 | Ze_ray(pr,k) = -999 |
|---|
| 496 | endif |
|---|
| 497 | if (z_vol(k) > 0) then |
|---|
| 498 | dBZe(pr,k) = 10*log10(z_vol(k))-a_to_vol(k)-g_to_vol(k) |
|---|
| 499 | Ze_non(pr,k) = 10*log10(z_vol(k)) |
|---|
| 500 | else |
|---|
| 501 | dBZe(pr,k) = -999 |
|---|
| 502 | Ze_non(pr,k) = -999 |
|---|
| 503 | endif |
|---|
| 504 | |
|---|
| 505 | enddo ! end loop of k (range gate) |
|---|
| 506 | ! Output array with gaseous absorption |
|---|
| 507 | if (g_to_vol_out_present) g_to_vol_out(:,pr) = g_to_vol |
|---|
| 508 | enddo ! end loop over pr (profile) |
|---|
| 509 | |
|---|
| 510 | end subroutine radar_simulator |
|---|
| 511 | |
|---|
