Changeset 1929 for palm/trunk/SOURCE/lpm_advec.f90
- Timestamp:
- Jun 9, 2016 4:25:25 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_advec.f90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r1889 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Put stochastic equation in an extra subroutine. 22 ! Set flag for stochastic equation to communicate whether a particle is near 23 ! topography. This case, memory and drift term are disabled in the Weil equation. 24 ! 25 ! Enable vertical logarithmic interpolation also above topography. This case, 26 ! set a lower limit for the friction velocity, as it can become very small 27 ! in narrow street canyons, leading to too large particle speeds. 22 28 ! 23 29 ! Former revisions: … … 110 116 ONLY: block_offset, c_0, dt_min_part, grid_particles, & 111 117 iran_part, log_z_z0, number_of_particles, number_of_sublayers, & 112 particles, particle_groups, offset_ocean_nzt, sgs_wfu_part, & 113 sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles, & 114 vertical_particle_advection, z0_av_global 118 particles, particle_groups, offset_ocean_nzt, sgs_wf_part, & 119 use_sgs_for_particles, vertical_particle_advection, z0_av_global 115 120 116 121 USE statistics, & … … 119 124 IMPLICIT NONE 120 125 121 INTEGER(iwp) :: agp !< 122 INTEGER(iwp) :: gp_outside_of_building(1:8) !< 123 INTEGER(iwp) :: i !< 124 INTEGER(iwp) :: ip !< 125 INTEGER(iwp) :: j !< 126 INTEGER(iwp) :: jp !< 127 INTEGER(iwp) :: k !< 128 INTEGER(iwp) :: kp !< 129 INTEGER(iwp) :: kw !< 130 INTEGER(iwp) :: n !< 131 INTEGER(iwp) :: nb !< 132 INTEGER(iwp) :: num_gp !< 133 134 INTEGER(iwp), DIMENSION(0:7) :: start_index !< 135 INTEGER(iwp), DIMENSION(0:7) :: end_index !< 136 137 REAL(wp) :: aa !< 138 REAL(wp) :: bb !< 139 REAL(wp) :: cc !< 140 REAL(wp) :: d_sum !< 141 REAL(wp) :: d_z_p_z0 !< 142 REAL(wp) :: dd !< 143 REAL(wp) :: de_dx_int_l !< 144 REAL(wp) :: de_dx_int_u !< 145 REAL(wp) :: de_dy_int_l !< 146 REAL(wp) :: de_dy_int_u !< 147 REAL(wp) :: de_dt !< 148 REAL(wp) :: de_dt_min !< 149 REAL(wp) :: de_dz_int_l !< 150 REAL(wp) :: de_dz_int_u !< 126 INTEGER(iwp) :: agp !< loop variable 127 INTEGER(iwp) :: gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography 128 INTEGER(iwp) :: i !< index variable along x 129 INTEGER(iwp) :: ip !< index variable along x 130 INTEGER(iwp) :: ilog !< index variable along x 131 INTEGER(iwp) :: j !< index variable along y 132 INTEGER(iwp) :: jp !< index variable along y 133 INTEGER(iwp) :: jlog !< index variable along y 134 INTEGER(iwp) :: k !< index variable along z 135 INTEGER(iwp) :: kp !< index variable along z 136 INTEGER(iwp) :: kw !< index variable along z 137 INTEGER(iwp) :: n !< loop variable over all particles in a grid box 138 INTEGER(iwp) :: nb !< block number particles are sorted in 139 INTEGER(iwp) :: num_gp !< number of adjacent grid points inside topography 140 141 INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block 142 INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block 143 144 REAL(wp) :: aa !< dummy argument for horizontal particle interpolation 145 REAL(wp) :: bb !< dummy argument for horizontal particle interpolation 146 REAL(wp) :: cc !< dummy argument for horizontal particle interpolation 147 REAL(wp) :: d_sum !< dummy argument for horizontal particle interpolation in case of topography 148 REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation 149 REAL(wp) :: dd !< dummy argument for horizontal particle interpolation 150 REAL(wp) :: de_dx_int_l !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level 151 REAL(wp) :: de_dx_int_u !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level 152 REAL(wp) :: de_dy_int_l !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level 153 REAL(wp) :: de_dy_int_u !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level 154 REAL(wp) :: de_dt !< temporal derivative of TKE experienced by the particle 155 REAL(wp) :: de_dt_min !< lower level for temporal TKE derivative 156 REAL(wp) :: de_dz_int_l !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level 157 REAL(wp) :: de_dz_int_u !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level 151 158 REAL(wp) :: diameter !< diamter of droplet 152 REAL(wp) :: diss_int_l !< 153 REAL(wp) :: diss_int_u !< 154 REAL(wp) :: dt_gap !< 155 REAL(wp) :: dt_particle_m !< 156 REAL(wp) :: e_int_l !< 157 REAL(wp) :: e_int_u !< 158 REAL(wp) :: e_mean_int !< 159 REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level 160 REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level 161 REAL(wp) :: dt_gap !< remaining time until particle time integration reaches LES time 162 REAL(wp) :: dt_particle_m !< previous particle time step 163 REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level 164 REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level 165 REAL(wp) :: e_mean_int !< horizontal mean TKE at particle height 159 166 REAL(wp) :: exp_arg !< 160 167 REAL(wp) :: exp_term !< 161 REAL(wp) :: gg !< 162 REAL(wp) :: height_p !< 168 REAL(wp) :: gg !< dummy argument for horizontal particle interpolation 169 REAL(wp) :: height_p !< dummy argument for logarithmic interpolation 163 170 REAL(wp) :: lagr_timescale !< Lagrangian timescale 164 REAL(wp) :: location(1:30,1:3) !< 171 REAL(wp) :: location(1:30,1:3) !< wall locations 172 REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation 165 173 REAL(wp) :: random_gauss !< 166 174 REAL(wp) :: RL !< Lagrangian autocorrelation coefficient … … 169 177 REAL(wp) :: rg3 !< Gaussian distributed random number 170 178 REAL(wp) :: sigma !< velocity standard deviation 171 REAL(wp) :: u_int_l !< 172 REAL(wp) :: u_int_u !< 173 REAL(wp) :: us_int !< 174 REAL(wp) :: v_int_l !< 175 REAL(wp) :: v_int_u !< 179 REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level 180 REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level 181 REAL(wp) :: us_int !< friction velocity at particle grid box 182 REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level 183 REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level 176 184 REAL(wp) :: vv_int !< 177 REAL(wp) :: w_int_l !< 178 REAL(wp) :: w_int_u !< 185 REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level 186 REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level 179 187 REAL(wp) :: w_s !< terminal velocity of droplets 180 REAL(wp) :: x !< 181 REAL(wp) :: y !< 182 REAL(wp) :: z_p !< 188 REAL(wp) :: x !< dummy argument for horizontal particle interpolation 189 REAL(wp) :: y !< dummy argument for horizontal particle interpolation 190 REAL(wp) :: z_p !< surface layer height (0.5 dz) 183 191 184 192 REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity … … 189 197 REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter 190 198 191 REAL(wp), DIMENSION(1:30) :: d_gp_pl !< 192 REAL(wp), DIMENSION(1:30) :: de_dxi !< 193 REAL(wp), DIMENSION(1:30) :: de_dyi !< 194 REAL(wp), DIMENSION(1:30) :: de_dzi !< 195 REAL(wp), DIMENSION(1:30) :: dissi !< 196 REAL(wp), DIMENSION(1:30) :: ei !< 197 199 REAL(wp), DIMENSION(1:30) :: d_gp_pl !< dummy argument for particle interpolation scheme in case of topography 200 REAL(wp), DIMENSION(1:30) :: de_dxi !< horizontal TKE gradient along x at adjacent wall 201 REAL(wp), DIMENSION(1:30) :: de_dyi !< horizontal TKE gradient along y at adjacent wall 202 REAL(wp), DIMENSION(1:30) :: de_dzi !< horizontal TKE gradient along z at adjacent wall 203 REAL(wp), DIMENSION(1:30) :: dissi !< dissipation at adjacent wall 204 REAL(wp), DIMENSION(1:30) :: ei !< TKE at adjacent wall 205 206 REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not 198 207 REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !< 199 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< 200 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< 201 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< 202 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< 203 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< 204 REAL(wp), DIMENSION(number_of_particles) :: e_int !< 205 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< 206 REAL(wp), DIMENSION(number_of_particles) :: log_z_z0_int !< 207 REAL(wp), DIMENSION(number_of_particles) :: u_int !< 208 REAL(wp), DIMENSION(number_of_particles) :: v_int !< 209 REAL(wp), DIMENSION(number_of_particles) :: w_int !< 210 REAL(wp), DIMENSION(number_of_particles) :: xv !< 211 REAL(wp), DIMENSION(number_of_particles) :: yv !< 212 REAL(wp), DIMENSION(number_of_particles) :: zv !< 213 214 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< 208 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position 209 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position 210 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position 211 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position 212 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step 213 REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position 214 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed 215 REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed 216 REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed 217 REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed 218 REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position 219 REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position 220 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 221 222 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers 215 223 216 224 CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' ) … … 236 244 j = jp + block_offset(nb)%j_off 237 245 k = kp + block_offset(nb)%k_off 246 238 247 239 248 ! … … 249 258 !-- First, check if particle is located below first vertical grid level 250 259 !-- (Prandtl-layer height) 251 IF ( constant_flux_layer .AND. particles(n)%z < z_p ) THEN 260 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 261 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 262 263 IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 252 264 ! 253 265 !-- Resolved-scale horizontal particle velocity is zero below z0. 254 IF ( particles(n)%z< z0_av_global ) THEN266 IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN 255 267 u_int(n) = 0.0_wp 256 268 ELSE 257 269 ! 258 270 !-- Determine the sublayer. Further used as index. 259 height_p = ( particles(n)%z- z0_av_global ) &271 height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) & 260 272 * REAL( number_of_sublayers, KIND=wp ) & 261 273 * d_z_p_z0 … … 263 275 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 264 276 !-- interpolate linearly between precalculated logarithm. 265 log_z_z0_int (n) = log_z_z0(INT(height_p))&277 log_z_z0_int = log_z_z0(INT(height_p)) & 266 278 + ( height_p - INT(height_p) ) & 267 279 * ( log_z_z0(INT(height_p)+1) & 268 280 - log_z_z0(INT(height_p)) & 269 281 ) 282 ! 283 !-- Limit friction velocity. In narrow canyons or holes the 284 !-- friction velocity can become very small, resulting in a too 285 !-- large particle speed. 286 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ), & 287 0.01_wp ) 270 288 ! 271 289 !-- Neutral solution is applied for all situations, e.g. also for … … 275 293 !-- as sensitivity studies revealed no significant effect of 276 294 !-- using the neutral solution also for un/stable situations. 277 !-- Calculated left and bottom index on u grid. 278 us_int = 0.5_wp * ( us(j,i) + us(j,i-1) ) 279 280 u_int(n) = -usws(j,i) / ( us_int * kappa + 1E-10_wp ) & 281 * log_z_z0_int(n) - u_gtrans 282 295 u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) & 296 * log_z_z0_int - u_gtrans 297 283 298 ENDIF 284 299 ! … … 308 323 ( u_int_u - u_int_l ) 309 324 ENDIF 325 310 326 ENDIF 311 327 … … 319 335 DO n = start_index(nb), end_index(nb) 320 336 321 IF ( constant_flux_layer .AND. particles(n)%z < z_p ) THEN 322 323 IF ( particles(n)%z < z0_av_global ) THEN 337 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 338 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 339 IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 340 341 IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN 324 342 ! 325 343 !-- Resolved-scale horizontal particle velocity is zero below z0. 326 344 v_int(n) = 0.0_wp 327 345 ELSE 346 ! 347 !-- Determine the sublayer. Further used as index. Please note, 348 !-- logarithmus can not be reused from above, as in in case of 349 !-- topography particle on u-grid can be above surface-layer height, 350 !-- whereas it can be below on v-grid. 351 height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) & 352 * REAL( number_of_sublayers, KIND=wp ) & 353 * d_z_p_z0 354 ! 355 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 356 !-- interpolate linearly between precalculated logarithm. 357 log_z_z0_int = log_z_z0(INT(height_p)) & 358 + ( height_p - INT(height_p) ) & 359 * ( log_z_z0(INT(height_p)+1) & 360 - log_z_z0(INT(height_p)) & 361 ) 362 ! 363 !-- Limit friction velocity. In narrow canyons or holes the 364 !-- friction velocity can become very small, resulting in a too 365 !-- large particle speed. 366 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), & 367 0.01_wp ) 328 368 ! 329 369 !-- Neutral solution is applied for all situations, e.g. also for … … 333 373 !-- as sensitivity studies revealed no significant effect of 334 374 !-- using the neutral solution also for un/stable situations. 335 !-- Calculated left and bottom index on v grid. 336 us_int = 0.5_wp * ( us(j,i) + us(j-1,i) ) 337 338 v_int(n) = -vsws(j,i) / ( us_int * kappa + 1E-10_wp ) & 339 * log_z_z0_int(n) - v_gtrans 340 ENDIF 375 v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) & 376 * log_z_z0_int - v_gtrans 377 378 ENDIF 379 341 380 ELSE 342 381 x = xv(n) - i * dx … … 361 400 ( v_int_u - v_int_l ) 362 401 ENDIF 402 363 403 ENDIF 364 404 … … 367 407 i = ip + block_offset(nb)%i_off 368 408 j = jp + block_offset(nb)%j_off 369 k = kp -1409 k = kp - 1 370 410 ! 371 411 !-- Same procedure for interpolation of the w velocity-component … … 535 575 ENDIF 536 576 577 ! 578 !-- Set flag for stochastic equation. 579 term_1_2(n) = 1.0_wp 580 537 581 ENDDO 538 582 ENDDO … … 541 585 542 586 DO n = 1, number_of_particles 543 544 587 i = particles(n)%x * ddx 545 588 j = particles(n)%y * ddy … … 575 618 de_dzi(num_gp) = de_dz(k,j,i) 576 619 ENDIF 577 578 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 579 THEN 620 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 580 621 num_gp = num_gp + 1 581 622 gp_outside_of_building(2) = 1 … … 603 644 ENDIF 604 645 605 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 606 THEN 646 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 607 647 num_gp = num_gp + 1 608 648 gp_outside_of_building(4) = 1 … … 617 657 ENDIF 618 658 619 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 620 THEN 659 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 621 660 num_gp = num_gp + 1 622 661 gp_outside_of_building(5) = 1 … … 631 670 ENDIF 632 671 633 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) & 634 THEN 672 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) THEN 635 673 num_gp = num_gp + 1 636 674 gp_outside_of_building(6) = 1 … … 645 683 ENDIF 646 684 647 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 648 THEN 685 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 649 686 num_gp = num_gp + 1 650 687 gp_outside_of_building(7) = 1 … … 659 696 ENDIF 660 697 661 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)& 662 THEN 698 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0) THEN 663 699 num_gp = num_gp + 1 664 700 gp_outside_of_building(8) = 1 … … 672 708 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 673 709 ENDIF 674 675 710 ! 676 711 !-- If all gridpoints are situated outside of a building, then the … … 685 720 dd = ( dx - x )**2 + ( dy - y )**2 686 721 gg = aa + bb + cc + dd 687 722 688 723 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 689 724 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 690 725 ) / ( 3.0_wp * gg ) 691 726 692 727 IF ( k == nzt ) THEN 693 728 e_int(n) = e_int_l … … 699 734 ) / ( 3.0_wp * gg ) 700 735 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 701 ( e_int_u - e_int_l )702 ENDIF 703 ! 736 ( e_int_u - e_int_l ) 737 ENDIF 738 ! 704 739 !-- Needed to avoid NaN particle velocities (this might not be 705 740 !-- required any more) … … 757 792 ( gg - dd ) * de_dz(k,j+1,i+1) & 758 793 ) / ( 3.0_wp * gg ) 759 794 760 795 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 761 796 de_dz_int(n) = de_dz_int_l … … 778 813 ( gg - dd ) * diss(k,j+1,i+1) & 779 814 ) / ( 3.0_wp * gg ) 780 815 781 816 IF ( k == nzt ) THEN 782 817 diss_int(n) = diss_int_l … … 790 825 ( diss_int_u - diss_int_l ) 791 826 ENDIF 792 827 ! 828 !-- Set flag for stochastic equation. 829 term_1_2(n) = 1.0_wp 830 793 831 ELSE 794 832 795 833 ! 796 834 !-- If wall between gridpoint 1 and gridpoint 5, then … … 810 848 811 849 IF ( gp_outside_of_building(5) == 1 .AND. & 812 gp_outside_of_building(1) == 0 ) THEN850 gp_outside_of_building(1) == 0 ) THEN 813 851 num_gp = num_gp + 1 814 852 location(num_gp,1) = i * dx + 0.5_wp * dx … … 893 931 de_dxi(num_gp) = de_dx(k,j,i) 894 932 de_dyi(num_gp) = 0.0_wp 895 de_dzi(num_gp) = de_dz(k,j,i) 933 de_dzi(num_gp) = de_dz(k,j,i) 896 934 ENDIF 897 935 … … 1076 1114 ENDIF 1077 1115 1078 ! 1116 ! 1079 1117 !-- If wall between gridpoint 6 and gridpoint 8, then 1080 1118 !-- Neumann boundary condition has to be applied … … 1092 1130 de_dzi(num_gp) = 0.0_wp 1093 1131 ENDIF 1094 1132 1095 1133 ! 1096 1134 !-- Carry out the interpolation 1097 1135 IF ( num_gp == 1 ) THEN 1098 ! 1136 ! 1099 1137 !-- If only one of the gridpoints is situated outside of the 1100 1138 !-- building, it follows that the values at the particle 1101 1139 !-- location are the same as the gridpoint values 1102 e_int(n) = ei(num_gp) 1103 diss_int(n) = dissi(num_gp) 1140 e_int(n) = ei(num_gp) 1141 diss_int(n) = dissi(num_gp) 1104 1142 de_dx_int(n) = de_dxi(num_gp) 1105 1143 de_dy_int(n) = de_dyi(num_gp) 1106 1144 de_dz_int(n) = de_dzi(num_gp) 1145 ! 1146 !-- Set flag for stochastic equation which disables calculation 1147 !-- of drift and memory term near topography. 1148 term_1_2(n) = 0.0_wp 1107 1149 ELSE IF ( num_gp > 1 ) THEN 1108 1150 1109 1151 d_sum = 0.0_wp 1110 ! 1152 ! 1111 1153 !-- Evaluation of the distances between the gridpoints 1112 1154 !-- contributing to the interpolated values, and the particle … … 1118 1160 d_sum = d_sum + d_gp_pl(agp) 1119 1161 ENDDO 1120 1162 1121 1163 ! 1122 1164 !-- Finally the interpolation can be carried out 1123 1165 e_int(n) = 0.0_wp 1124 1166 diss_int(n) = 0.0_wp 1125 de_dx_int(n) = 0.0_wp 1126 de_dy_int(n) = 0.0_wp 1127 de_dz_int(n) = 0.0_wp 1167 de_dx_int(n) = 0.0_wp 1168 de_dy_int(n) = 0.0_wp 1169 de_dz_int(n) = 0.0_wp 1128 1170 DO agp = 1, num_gp 1129 1171 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & … … 1138 1180 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1139 1181 ENDDO 1140 1141 ENDIF 1142 1182 1183 ENDIF 1184 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1185 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1186 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1187 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1188 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1189 ! 1190 !-- Set flag for stochastic equation which disables calculation 1191 !-- of drift and memory term near topography. 1192 term_1_2(n) = 0.0_wp 1143 1193 ENDIF 1144 1194 ENDDO … … 1209 1259 ! 1210 1260 !-- Calculate the Lagrangian timescale according to Weil et al. (2004). 1211 lagr_timescale = ( 4.0_wp * e_int(n) ) / &1212 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )1261 lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / & 1262 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp ) 1213 1263 1214 1264 ! … … 1233 1283 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 1234 1284 !-- from becoming unrealistically large. 1235 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf u_part * e_int(n)) * &1285 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1236 1286 ( rg(n,1) - 1.0_wp ) 1237 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf v_part * e_int(n)) * &1287 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1238 1288 ( rg(n,2) - 1.0_wp ) 1239 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf w_part * e_int(n)) * &1289 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1240 1290 ( rg(n,3) - 1.0_wp ) 1241 1291 … … 1267 1317 ENDIF 1268 1318 1269 particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 * & 1270 diss_int(n) * particles(n)%rvar1 * dt_particle(n) / & 1271 ( 4.0_wp * sgs_wfu_part * e_int(n) ) + & 1272 ( 2.0_wp * sgs_wfu_part * de_dt * & 1273 particles(n)%rvar1 / & 1274 ( 2.0_wp * sgs_wfu_part * e_int(n) ) + & 1275 de_dx_int(n) & 1276 ) * dt_particle(n) / 2.0_wp + & 1277 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1278 ( rg(n,1) - 1.0_wp ) * & 1279 SQRT( dt_particle(n) ) 1280 1281 particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 * & 1282 diss_int(n) * particles(n)%rvar2 * dt_particle(n) / & 1283 ( 4.0_wp * sgs_wfv_part * e_int(n) ) + & 1284 ( 2.0_wp * sgs_wfv_part * de_dt * & 1285 particles(n)%rvar2 / & 1286 ( 2.0_wp * sgs_wfv_part * e_int(n) ) + & 1287 de_dy_int(n) & 1288 ) * dt_particle(n) / 2.0_wp + & 1289 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1290 ( rg(n,2) - 1.0_wp ) * & 1291 SQRT( dt_particle(n) ) 1292 1293 particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 * & 1294 diss_int(n) * particles(n)%rvar3 * dt_particle(n) / & 1295 ( 4.0_wp * sgs_wfw_part * e_int(n) ) + & 1296 ( 2.0_wp * sgs_wfw_part * de_dt * & 1297 particles(n)%rvar3 / & 1298 ( 2.0_wp * sgs_wfw_part * e_int(n) ) + & 1299 de_dz_int(n) & 1300 ) * dt_particle(n) / 2.0_wp + & 1301 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1302 ( rg(n,3) - 1.0_wp ) * & 1303 SQRT( dt_particle(n) ) 1319 CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n), & 1320 de_dx_int(n), de_dt, diss_int(n), & 1321 dt_particle(n), rg(n,1), term_1_2(n) ) 1322 1323 CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n), & 1324 de_dy_int(n), de_dt, diss_int(n), & 1325 dt_particle(n), rg(n,2), term_1_2(n) ) 1326 1327 CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n), & 1328 de_dz_int(n), de_dt, diss_int(n), & 1329 dt_particle(n), rg(n,3), term_1_2(n) ) 1304 1330 1305 1331 ENDIF 1332 1306 1333 u_int(n) = u_int(n) + particles(n)%rvar1 1307 1334 v_int(n) = v_int(n) + particles(n)%rvar2 1308 1335 w_int(n) = w_int(n) + particles(n)%rvar3 1309 1310 1336 ! 1311 1337 !-- Store the SGS TKE of the current timelevel which is needed for … … 1507 1533 CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' ) 1508 1534 1535 1509 1536 END SUBROUTINE lpm_advec 1537 1538 ! Description: 1539 ! ------------ 1540 !> Calculation of subgrid-scale particle speed using the stochastic model 1541 !> of Weil et al. (2004, JAS, 61, 2877-2887). 1542 !------------------------------------------------------------------------------! 1543 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n, & 1544 dt_n, rg_n, fac ) 1545 1546 USE kinds 1547 1548 USE particle_attributes, & 1549 ONLY: c_0, sgs_wf_part 1550 1551 IMPLICIT NONE 1552 1553 REAL(wp) :: a1 !< dummy argument 1554 REAL(wp) :: dedt_n !< time derivative of TKE at particle position 1555 REAL(wp) :: dedxi_n !< horizontal derivative of TKE at particle position 1556 REAL(wp) :: diss_n !< dissipation at particle position 1557 REAL(wp) :: dt_n !< particle timestep 1558 REAL(wp) :: e_n !< TKE at particle position 1559 REAL(wp) :: fac !< flag to identify adjacent topography 1560 REAL(wp) :: fs_n !< weighting factor to prevent that subgrid-scale particle speed becomes too large 1561 REAL(wp) :: sgs_w !< constant (1/3) 1562 REAL(wp) :: rg_n !< random number 1563 REAL(wp) :: term1 !< memory term 1564 REAL(wp) :: term2 !< drift correction term 1565 REAL(wp) :: term3 !< random term 1566 REAL(wp) :: v_sgs !< subgrid-scale velocity component 1567 1568 ! 1569 !-- Please note, terms 1 and 2 (drift and memory term, respectively) are 1570 !-- multiplied by a flag to switch of both terms near topography. 1571 !-- This is necessary, as both terms may cause a subgrid-scale velocity build up 1572 !-- if particles are trapped in regions with very small TKE, e.g. in narrow street 1573 !-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are 1574 !-- disabled if one of the adjacent grid points belongs to topography. 1575 !-- Moreover, in this case, the previous subgrid-scale component is also set 1576 !-- to zero. 1577 1578 a1 = fs_n * c_0 * diss_n 1579 ! 1580 !-- Memory term 1581 term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp ) & 1582 * fac 1583 ! 1584 !-- Drift correction term 1585 term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n & 1586 * fac 1587 ! 1588 !-- Random term 1589 term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n ) 1590 ! 1591 !-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous 1592 !-- subgrid-scale velocity component is set to zero, in order to prevent a 1593 !-- velocity build-up. 1594 1595 !-- This case, set also previous subgrid-scale component to zero. 1596 v_sgs = v_sgs * fac + term1 + term2 + term3 1597 1598 END SUBROUTINE weil_stochastic_eq
Note: See TracChangeset
for help on using the changeset viewer.
