| 1 | subroutine radar_simulator( & |
|---|
| 2 | hp, & |
|---|
| 3 | nprof,ngate, & |
|---|
| 4 | undef, & |
|---|
| 5 | hgt_matrix,hm_matrix,re_matrix,Np_matrix, & |
|---|
| 6 | p_matrix,t_matrix,rh_matrix, & |
|---|
| 7 | Ze_non,Ze_ray,a_to_vol,g_to_vol,dBZe, & |
|---|
| 8 | g_to_vol_in,g_to_vol_out) |
|---|
| 9 | |
|---|
| 10 | use m_mrgrnk |
|---|
| 11 | use array_lib |
|---|
| 12 | use math_lib |
|---|
| 13 | use optics_lib |
|---|
| 14 | use radar_simulator_types |
|---|
| 15 | use scale_LUTs_io |
|---|
| 16 | implicit none |
|---|
| 17 | |
|---|
| 18 | ! Purpose: |
|---|
| 19 | |
|---|
| 20 | ! Simulates a vertical profile of radar reflectivity |
|---|
| 21 | ! Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand. |
|---|
| 22 | ! but has been substantially modified since that time by |
|---|
| 23 | ! Laura Fowler and Roger Marchand (see modifications below). |
|---|
| 24 | |
|---|
| 25 | ! Inputs: |
|---|
| 26 | |
|---|
| 27 | ! [hp] structure that defines hydrometeor types and other radar properties |
|---|
| 28 | |
|---|
| 29 | ! [nprof] number of hydrometeor profiles |
|---|
| 30 | ! [ngate] number of vertical layers |
|---|
| 31 | |
|---|
| 32 | ! [undef] missing data value |
|---|
| 33 | ! (The following 5 arrays must be in order from closest to the radar |
|---|
| 34 | ! to farthest...) |
|---|
| 35 | |
|---|
| 36 | ! [hgt_matrix] height of hydrometeors (km) |
|---|
| 37 | ! [p_matrix] pressure profile (hPa) |
|---|
| 38 | ! [t_matrix] temperature profile (K) |
|---|
| 39 | ! [rh_matrix] relative humidity profile (%) -- only needed if gaseous aborption calculated. |
|---|
| 40 | |
|---|
| 41 | ! [hm_matrix] table of hydrometeor mixing rations (g/kg) |
|---|
| 42 | ! [re_matrix] table of hydrometeor effective radii. 0 ==> use defaults. (units=microns) |
|---|
| 43 | ! [Np_matrix] table of hydrometeor number concentration. 0 ==> use defaults. (units = 1/kg) |
|---|
| 44 | |
|---|
| 45 | ! Outputs: |
|---|
| 46 | |
|---|
| 47 | ! [Ze_non] radar reflectivity without attenuation (dBZ) |
|---|
| 48 | ! [Ze_ray] Rayleigh reflectivity (dBZ) |
|---|
| 49 | ! [h_atten_to_vol] attenuation by hydromets, radar to vol (dB) |
|---|
| 50 | ! [g_atten_to_vol] gaseous atteunation, radar to vol (dB) |
|---|
| 51 | ! [dBZe] effective radar reflectivity factor (dBZ) |
|---|
| 52 | |
|---|
| 53 | ! Optional: |
|---|
| 54 | ! [g_to_vol_in] integrated atten due to gases, r>v (dB). |
|---|
| 55 | ! If present then is used as gaseous absorption, independently of the |
|---|
| 56 | ! value in use_gas_abs |
|---|
| 57 | ! [g_to_vol_out] integrated atten due to gases, r>v (dB). |
|---|
| 58 | ! If present then gaseous absorption for each profile is returned here. |
|---|
| 59 | |
|---|
| 60 | ! Created: |
|---|
| 61 | ! 11/28/2005 John Haynes (haynes@atmos.colostate.edu) |
|---|
| 62 | |
|---|
| 63 | ! Modified: |
|---|
| 64 | ! 09/2006 placed into subroutine form (Roger Marchand,JMH) |
|---|
| 65 | ! 08/2007 added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand) |
|---|
| 66 | ! 01/2008 'Do while' to determine if hydrometeor(s) present in volume |
|---|
| 67 | ! changed for vectorization purposes (A. Bodas-Salcedo) |
|---|
| 68 | |
|---|
| 69 | ! 07/2010 V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT) |
|---|
| 70 | ! ... All hydrometeor and radar simulator properties now included in hp structure |
|---|
| 71 | ! ... hp structure should be initialized by call to radar_simulator_init prior |
|---|
| 72 | ! ... to calling this subroutine. |
|---|
| 73 | ! Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added |
|---|
| 74 | ! ... Changes implement by Roj Marchand following work by Laura Fowler |
|---|
| 75 | |
|---|
| 76 | ! 10/2011 Modified ngate loop to go in either direction depending on flag |
|---|
| 77 | ! hp%radar_at_layer_one. This affects the direction in which attenuation is summed. |
|---|
| 78 | |
|---|
| 79 | ! Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple |
|---|
| 80 | ! summation. (Roger Marchand) |
|---|
| 81 | |
|---|
| 82 | |
|---|
| 83 | ! ----- INPUTS ----- |
|---|
| 84 | |
|---|
| 85 | logical, parameter :: DO_LUT_TEST = .false. |
|---|
| 86 | logical, parameter :: DO_NP_TEST = .false. |
|---|
| 87 | |
|---|
| 88 | type(class_param), intent(inout) :: hp |
|---|
| 89 | |
|---|
| 90 | integer, intent(in) :: nprof,ngate |
|---|
| 91 | |
|---|
| 92 | real undef |
|---|
| 93 | real*8, dimension(nprof,ngate), intent(in) :: & |
|---|
| 94 | hgt_matrix, p_matrix,t_matrix,rh_matrix |
|---|
| 95 | |
|---|
| 96 | real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix |
|---|
| 97 | real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix |
|---|
| 98 | real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: Np_matrix |
|---|
| 99 | |
|---|
| 100 | ! ----- OUTPUTS ----- |
|---|
| 101 | real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, & |
|---|
| 102 | g_to_vol,dBZe,a_to_vol |
|---|
| 103 | |
|---|
| 104 | ! ----- OPTIONAL ----- |
|---|
| 105 | real*8, optional, dimension(nprof,ngate) :: & |
|---|
| 106 | g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input |
|---|
| 107 | ! the same gaseous absorption in different calls. Optional to allow compatibility |
|---|
| 108 | ! with original version. A. Bodas April 2008. |
|---|
| 109 | ! real*8, dimension(nprof,ngate) :: kr_matrix |
|---|
| 110 | |
|---|
| 111 | ! ----- INTERNAL ----- |
|---|
| 112 | |
|---|
| 113 | real, parameter :: one_third = 1.0/3.0 |
|---|
| 114 | real*8 :: t_kelvin |
|---|
| 115 | integer :: & |
|---|
| 116 | phase, & ! 0=liquid, 1=ice |
|---|
| 117 | ns ! number of discrete drop sizes |
|---|
| 118 | |
|---|
| 119 | logical :: hydro ! true=hydrometeor in vol, false=none |
|---|
| 120 | real*8 :: & |
|---|
| 121 | rho_a, & ! air density (kg m^-3) |
|---|
| 122 | gases ! function: 2-way gas atten (dB/km) |
|---|
| 123 | |
|---|
| 124 | real*8, dimension(:), allocatable :: & |
|---|
| 125 | Di, Deq, & ! discrete drop sizes (um) |
|---|
| 126 | Ni, & ! discrete concentrations (cm^-3 um^-1) |
|---|
| 127 | rhoi ! discrete densities (kg m^-3) |
|---|
| 128 | |
|---|
| 129 | real*8, dimension(nprof, ngate) :: & |
|---|
| 130 | z_vol, & ! effective reflectivity factor (mm^6/m^3) |
|---|
| 131 | z_ray, & ! reflectivity factor, Rayleigh only (mm^6/m^3) |
|---|
| 132 | kr_vol, & ! attenuation coefficient hydro (dB/km) |
|---|
| 133 | g_vol ! attenuation coefficient gases (dB/km) |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | integer,parameter :: KR8 = selected_real_kind(15,300) |
|---|
| 137 | real*8, parameter :: xx = -1.0_KR8 |
|---|
| 138 | real*8, dimension(:), allocatable :: xxa |
|---|
| 139 | real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm |
|---|
| 140 | real*8 :: half_a_atten_current,half_a_atten_above |
|---|
| 141 | real*8 :: half_g_atten_current,half_g_atten_above |
|---|
| 142 | integer*4 :: tp, i, j, k, pr, itt, iff |
|---|
| 143 | |
|---|
| 144 | real*8 step,base, Np |
|---|
| 145 | integer*4 iRe_type,n,max_bin |
|---|
| 146 | |
|---|
| 147 | integer start_gate,end_gate,d_gate |
|---|
| 148 | |
|---|
| 149 | logical :: g_to_vol_in_present, g_to_vol_out_present |
|---|
| 150 | |
|---|
| 151 | ! Logicals to avoid calling present within the loops |
|---|
| 152 | g_to_vol_in_present = present(g_to_vol_in) |
|---|
| 153 | g_to_vol_out_present = present(g_to_vol_out) |
|---|
| 154 | |
|---|
| 155 | ! load scaling matricies from disk -- but only the first time this subroutine is called |
|---|
| 156 | |
|---|
| 157 | if(hp%load_scale_LUTs) then |
|---|
| 158 | call load_scale_LUTs(hp) |
|---|
| 159 | hp%load_scale_LUTs=.false. |
|---|
| 160 | hp%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run |
|---|
| 161 | endif |
|---|
| 162 | |
|---|
| 163 | pi = acos(-1.0) |
|---|
| 164 | |
|---|
| 165 | ! ----- Initialisation ----- |
|---|
| 166 | g_to_vol = 0.0 |
|---|
| 167 | a_to_vol = 0.0 |
|---|
| 168 | z_vol = 0.0 |
|---|
| 169 | z_ray = 0.0 |
|---|
| 170 | kr_vol = 0.0 |
|---|
| 171 | |
|---|
| 172 | ! // loop over each range gate (ngate) ... starting with layer closest to the radar ! |
|---|
| 173 | if(hp%radar_at_layer_one) then |
|---|
| 174 | start_gate=1 |
|---|
| 175 | end_gate=ngate |
|---|
| 176 | d_gate=1 |
|---|
| 177 | else |
|---|
| 178 | start_gate=ngate |
|---|
| 179 | end_gate=1 |
|---|
| 180 | d_gate=-1 |
|---|
| 181 | endif |
|---|
| 182 | DO k=start_gate,end_gate,d_gate |
|---|
| 183 | ! // loop over each profile (nprof) |
|---|
| 184 | DO pr=1,nprof |
|---|
| 185 | t_kelvin = t_matrix(pr,k) |
|---|
| 186 | ! :: determine if hydrometeor(s) present in volume |
|---|
| 187 | hydro = .false. |
|---|
| 188 | DO j=1,hp%nhclass |
|---|
| 189 | if ((hm_matrix(j,pr,k) > 1E-12) .AND. (hp%dtype(j) > 0)) then |
|---|
| 190 | hydro = .true. |
|---|
| 191 | exit |
|---|
| 192 | endif |
|---|
| 193 | enddo |
|---|
| 194 | |
|---|
| 195 | ! :: if there is hydrometeor in the volume |
|---|
| 196 | if (hydro) then |
|---|
| 197 | |
|---|
| 198 | rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin)) |
|---|
| 199 | ! :: loop over hydrometeor type |
|---|
| 200 | DO tp=1,hp%nhclass |
|---|
| 201 | if (hm_matrix(tp,pr,k) <= 1E-12) cycle |
|---|
| 202 | phase = hp%phase(tp) |
|---|
| 203 | if (phase==0) then |
|---|
| 204 | itt = infind(hp%mt_ttl,t_kelvin) |
|---|
| 205 | else |
|---|
| 206 | itt = infind(hp%mt_tti,t_kelvin) |
|---|
| 207 | endif |
|---|
| 208 | if (re_matrix(tp,pr,k).eq.0) then |
|---|
| 209 | call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, & |
|---|
| 210 | hp%dtype(tp),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), & |
|---|
| 211 | hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),Re) |
|---|
| 212 | re_matrix(tp,pr,k)=Re |
|---|
| 213 | else |
|---|
| 214 | if (Np_matrix(tp,pr,k)>0) then |
|---|
| 215 | PRINT *, 'Warning: Re and Np set for the same ', & |
|---|
| 216 | 'volume & hydrometeor type. Np is being ignored.' |
|---|
| 217 | endif |
|---|
| 218 | Re = re_matrix(tp,pr,k) |
|---|
| 219 | endif |
|---|
| 220 | |
|---|
| 221 | iRe_type=1 |
|---|
| 222 | if(Re.gt.0) then |
|---|
| 223 | ! determine index in to scale LUT |
|---|
| 224 | |
|---|
| 225 | ! distance between Re points (defined by "base" and "step") for |
|---|
| 226 | ! each interval of size Re_BIN_LENGTH |
|---|
| 227 | ! Integer asignment, avoids calling floor intrinsic |
|---|
| 228 | n=Re/Re_BIN_LENGTH |
|---|
| 229 | if (n>=Re_MAX_BIN) n=Re_MAX_BIN-1 |
|---|
| 230 | step=hp%step_list(n+1) |
|---|
| 231 | base=hp%base_list(n+1) |
|---|
| 232 | iRe_type=Re/step |
|---|
| 233 | if (iRe_type.lt.1) iRe_type=1 |
|---|
| 234 | |
|---|
| 235 | Re=step*(iRe_type+0.5) ! set value of Re to closest value allowed in LUT. |
|---|
| 236 | iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step) |
|---|
| 237 | |
|---|
| 238 | ! make sure iRe_type is within bounds |
|---|
| 239 | if (iRe_type.ge.nRe_types) then |
|---|
| 240 | ! write(*,*) 'Warning: size of Re exceed value permitted ', & |
|---|
| 241 | ! 'in Look-Up Table (LUT). Will calculate. ' |
|---|
| 242 | ! no scaling allowed |
|---|
| 243 | iRe_type=nRe_types |
|---|
| 244 | hp%Z_scale_flag(tp,itt,iRe_type)=.false. |
|---|
| 245 | else |
|---|
| 246 | ! set value in re_matrix to closest values in LUT |
|---|
| 247 | if (.not. DO_LUT_TEST) re_matrix(tp,pr,k)=Re |
|---|
| 248 | endif |
|---|
| 249 | endif |
|---|
| 250 | ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them |
|---|
| 251 | ! if not we will calculate Ze, Zr, and Kr from the distribution parameters |
|---|
| 252 | if( (.not. hp%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST) then |
|---|
| 253 | ! :: create a distribution of hydrometeors within volume |
|---|
| 254 | select case(hp%dtype(tp)) |
|---|
| 255 | case(4) |
|---|
| 256 | ns = 1 |
|---|
| 257 | allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) |
|---|
| 258 | Di = hp%p1(tp) |
|---|
| 259 | Ni = 0. |
|---|
| 260 | case default |
|---|
| 261 | ns = nd ! constant defined in radar_simulator_types.f90 |
|---|
| 262 | allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns)) |
|---|
| 263 | Di = hp%D |
|---|
| 264 | Ni = 0. |
|---|
| 265 | end select |
|---|
| 266 | call dsd(hm_matrix(tp,pr,k),re_matrix(tp,pr,k),Np_matrix(tp,pr,k), & |
|---|
| 267 | Di,Ni,ns,hp%dtype(tp),rho_a,t_kelvin, & |
|---|
| 268 | hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), & |
|---|
| 269 | hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp)) |
|---|
| 270 | |
|---|
| 271 | ! calculate particle density |
|---|
| 272 | if (phase == 1) then |
|---|
| 273 | if (hp%rho(tp) < 0) then |
|---|
| 274 | ! Use equivalent volume spheres. |
|---|
| 275 | hp%rho_eff(tp,1:ns,iRe_type) = 917 ! solid ice == equivalent volume approach |
|---|
| 276 | Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) ) * 1E6 |
|---|
| 277 | ! alternative is to comment out above two lines and use the following block |
|---|
| 278 | ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density |
|---|
| 279 | |
|---|
| 280 | ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3) !MG Mie approach |
|---|
| 281 | |
|---|
| 282 | ! as the particle size gets small it is possible that the mass to size relationship of |
|---|
| 283 | ! (given by power law in hclass.data) can produce impossible results |
|---|
| 284 | ! where the mass is larger than a solid sphere of ice. |
|---|
| 285 | ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere. |
|---|
| 286 | ! do i=1,ns |
|---|
| 287 | ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then |
|---|
| 288 | ! hp%rho_eff(tp,i,iRe_type) = 917 |
|---|
| 289 | ! endif |
|---|
| 290 | ! enddo |
|---|
| 291 | else |
|---|
| 292 | ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3). |
|---|
| 293 | hp%rho_eff(tp,1:ns,iRe_type) = 917 |
|---|
| 294 | Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0)) |
|---|
| 295 | ! alternative ... coment out above two lines and use the following for MG-Mie |
|---|
| 296 | ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp) !MG Mie approach |
|---|
| 297 | endif |
|---|
| 298 | else |
|---|
| 299 | ! I assume here that water phase droplets are spheres. |
|---|
| 300 | ! hp%rho should be ~ 1000 or hp%apm=524 .AND. hp%bpm=3 |
|---|
| 301 | Deq = Di |
|---|
| 302 | endif |
|---|
| 303 | |
|---|
| 304 | ! calculate effective reflectivity factor of volume |
|---|
| 305 | xxa = -9.9 |
|---|
| 306 | rhoi = hp%rho_eff(tp,1:ns,iRe_type) |
|---|
| 307 | call zeff(hp%freq,Deq,Ni,ns,hp%k2,t_kelvin,phase,hp%do_ray, & |
|---|
| 308 | ze,zr,kr,xxa,xxa,rhoi) |
|---|
| 309 | |
|---|
| 310 | ! test code ... compare Np value input to routine with sum of DSD |
|---|
| 311 | ! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation |
|---|
| 312 | ! not just the DSD representation given by Ni |
|---|
| 313 | if(Np_matrix(tp,pr,k)>0 .AND. DO_NP_TEST ) then |
|---|
| 314 | Np = path_integral(Ni,Di,1,ns-1)/rho_a*1E6 |
|---|
| 315 | ! Note: Representation is not great or small Re < 2 |
|---|
| 316 | if( (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)>0.1 ) then |
|---|
| 317 | write(*,*) 'Error: Np input does not match sum(N)' |
|---|
| 318 | write(*,*) tp,pr,k,Re,Ni(1),Ni(ns),10*log10(ze) |
|---|
| 319 | write(*,*) Np_matrix(tp,pr,k),Np,(Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k) |
|---|
| 320 | write(*,*) |
|---|
| 321 | endif |
|---|
| 322 | endif |
|---|
| 323 | |
|---|
| 324 | deallocate(Di,Ni,rhoi,xxa,Deq) |
|---|
| 325 | |
|---|
| 326 | ! LUT test code |
|---|
| 327 | ! This segment of code compares full calculation to scaling result |
|---|
| 328 | if ( hp%Z_scale_flag(tp,itt,iRe_type) .AND. DO_LUT_TEST ) then |
|---|
| 329 | scale_factor=rho_a*hm_matrix(tp,pr,k) |
|---|
| 330 | ! if more than 2 dBZe difference print error message/parameters. |
|---|
| 331 | if ( abs(10*log10(ze) - 10*log10(hp%Ze_scaled(tp,itt,iRe_type) * & |
|---|
| 332 | scale_factor)) > 2 ) then |
|---|
| 333 | write(*,*) 'Roj Error: ',tp,itt,iRe_type,hp%Z_scale_flag(tp,itt,iRe_type),n,step,base |
|---|
| 334 | write(*,*) 10*log10(ze),10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor) |
|---|
| 335 | write(*,*) hp%Ze_scaled(tp,itt,iRe_type),scale_factor |
|---|
| 336 | write(*,*) re_matrix(tp,pr,k),Re |
|---|
| 337 | write(*,*) |
|---|
| 338 | endif |
|---|
| 339 | endif |
|---|
| 340 | |
|---|
| 341 | else ! can use z scaling |
|---|
| 342 | scale_factor=rho_a*hm_matrix(tp,pr,k) |
|---|
| 343 | zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 344 | ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 345 | kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor |
|---|
| 346 | endif ! end z_scaling |
|---|
| 347 | |
|---|
| 348 | kr_vol(pr,k) = kr_vol(pr,k) + kr |
|---|
| 349 | z_vol(pr,k) = z_vol(pr,k) + ze |
|---|
| 350 | z_ray(pr,k) = z_ray(pr,k) + zr |
|---|
| 351 | |
|---|
| 352 | ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can |
|---|
| 353 | if ( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then |
|---|
| 354 | if (iRe_type>1) then |
|---|
| 355 | scale_factor=rho_a*hm_matrix(tp,pr,k) |
|---|
| 356 | hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor |
|---|
| 357 | hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor |
|---|
| 358 | hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor |
|---|
| 359 | hp%Z_scale_flag(tp,itt,iRe_type) = .true. |
|---|
| 360 | hp%Z_scale_added_flag(tp,itt,iRe_type)=.true. |
|---|
| 361 | endif |
|---|
| 362 | endif |
|---|
| 363 | |
|---|
| 364 | enddo ! end loop of tp (hydrometeor type) |
|---|
| 365 | |
|---|
| 366 | else |
|---|
| 367 | ! :: volume is hydrometeor-free |
|---|
| 368 | kr_vol(pr,k) = 0 |
|---|
| 369 | z_vol(pr,k) = undef |
|---|
| 370 | z_ray(pr,k) = undef |
|---|
| 371 | endif |
|---|
| 372 | |
|---|
| 373 | ! :: attenuation due to hydrometeors between radar and volume |
|---|
| 374 | |
|---|
| 375 | ! NOTE old scheme integrates attenuation only for the layers ABOVE |
|---|
| 376 | ! the current layer ... i.e. 1 to k-1 rather than 1 to k ... |
|---|
| 377 | ! which may be a problem. ROJ |
|---|
| 378 | ! in the new scheme I assign half the attenuation to the current layer |
|---|
| 379 | if(d_gate==1) then |
|---|
| 380 | ! dheight calcuations assumes hgt_matrix points are the cell mid-points. |
|---|
| 381 | if (k>2) then |
|---|
| 382 | ! add to previous value to half of above layer + half of current layer |
|---|
| 383 | a_to_vol(pr,k)= a_to_vol(pr,k-1) + & |
|---|
| 384 | (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) |
|---|
| 385 | else |
|---|
| 386 | a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1)) |
|---|
| 387 | endif |
|---|
| 388 | else ! d_gate==-1 |
|---|
| 389 | if(k<ngate) then |
|---|
| 390 | ! add to previous value half of above layer + half of current layer |
|---|
| 391 | a_to_vol(pr,k) = a_to_vol(pr,k+1) + & |
|---|
| 392 | (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) |
|---|
| 393 | else |
|---|
| 394 | a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)) |
|---|
| 395 | endif |
|---|
| 396 | endif |
|---|
| 397 | |
|---|
| 398 | ! :: attenuation due to gaseous absorption between radar and volume |
|---|
| 399 | if (g_to_vol_in_present) then |
|---|
| 400 | g_to_vol(pr,k) = g_to_vol_in(pr,k) |
|---|
| 401 | else |
|---|
| 402 | if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .AND. (pr == 1)) ) then |
|---|
| 403 | g_vol(pr,k) = gases(p_matrix(pr,k),t_kelvin,rh_matrix(pr,k),hp%freq) |
|---|
| 404 | if (d_gate==1) then |
|---|
| 405 | if (k>1) then |
|---|
| 406 | ! add to previous value to half of above layer + half of current layer |
|---|
| 407 | g_to_vol(pr,k) = g_to_vol(pr,k-1) + & |
|---|
| 408 | 0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) |
|---|
| 409 | else |
|---|
| 410 | g_to_vol(pr,k)= 0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1)) |
|---|
| 411 | endif |
|---|
| 412 | else ! d_gate==-1 |
|---|
| 413 | if (k<ngate) then |
|---|
| 414 | ! add to previous value to half of above layer + half of current layer |
|---|
| 415 | g_to_vol(pr,k) = g_to_vol(pr,k+1) + & |
|---|
| 416 | 0.5*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) |
|---|
| 417 | else |
|---|
| 418 | g_to_vol(pr,k)= 0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)) |
|---|
| 419 | endif |
|---|
| 420 | endif |
|---|
| 421 | elseif(hp%use_gas_abs == 2) then |
|---|
| 422 | ! using value calculated for the first column |
|---|
| 423 | g_to_vol(pr,k) = g_to_vol(1,k) |
|---|
| 424 | elseif (hp%use_gas_abs == 0) then |
|---|
| 425 | g_to_vol(pr,k) = 0 |
|---|
| 426 | endif |
|---|
| 427 | endif |
|---|
| 428 | |
|---|
| 429 | ! Compute Rayleigh reflectivity, and full, attenuated reflectivity |
|---|
| 430 | if ((hp%do_ray == 1) .AND. (z_ray(pr,k) > 0)) then |
|---|
| 431 | Ze_ray(pr,k) = 10*log10(z_ray(pr,k)) |
|---|
| 432 | else |
|---|
| 433 | Ze_ray(pr,k) = undef |
|---|
| 434 | endif |
|---|
| 435 | if (z_vol(pr,k) > 0) then |
|---|
| 436 | Ze_non(pr,k) = 10*log10(z_vol(pr,k)) |
|---|
| 437 | dBZe(pr,k) = Ze_non(pr,k)-a_to_vol(pr,k)-g_to_vol(pr,k) |
|---|
| 438 | else |
|---|
| 439 | dBZe(pr,k) = undef |
|---|
| 440 | Ze_non(pr,k) = undef |
|---|
| 441 | endif |
|---|
| 442 | |
|---|
| 443 | enddo ! end loop over pr (profile) |
|---|
| 444 | |
|---|
| 445 | enddo ! end loop of k (range gate) |
|---|
| 446 | |
|---|
| 447 | ! Output array with gaseous absorption |
|---|
| 448 | if (g_to_vol_out_present) g_to_vol_out = g_to_vol |
|---|
| 449 | |
|---|
| 450 | ! save any updates made |
|---|
| 451 | if (hp%update_scale_LUTs) call save_scale_LUTs(hp) |
|---|
| 452 | |
|---|
| 453 | end subroutine radar_simulator |
|---|
