Changeset 2232 for palm/trunk/SOURCE/lpm_advec.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_advec.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r2101 r2232 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 109 109 110 110 USE arrays_3d, & 111 ONLY: de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw111 ONLY: de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw 112 112 113 113 USE cpulog … … 123 123 124 124 USE indices, & 125 ONLY: nzb, nzb_ s_inner, nzt125 ONLY: nzb, nzb_max, nzt, wall_flags_0 126 126 127 127 USE kinds … … 136 136 ONLY: hom 137 137 138 USE surface_mod, & 139 ONLY: surf_def_h, surf_lsm_h, surf_usm_h 140 138 141 IMPLICIT NONE 139 142 … … 147 150 INTEGER(iwp) :: jlog !< index variable along y 148 151 INTEGER(iwp) :: k !< index variable along z 152 INTEGER(iwp) :: k_wall !< vertical index of topography top 149 153 INTEGER(iwp) :: kp !< index variable along z 150 154 INTEGER(iwp) :: kw !< index variable along z … … 152 156 INTEGER(iwp) :: nb !< block number particles are sorted in 153 157 INTEGER(iwp) :: num_gp !< number of adjacent grid points inside topography 158 INTEGER(iwp) :: surf_start !< Index on surface data-type for current grid box 154 159 155 160 INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block … … 194 199 REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level 195 200 REAL(wp) :: us_int !< friction velocity at particle grid box 201 REAL(wp) :: usws_int !< surface momentum flux (u component) at particle grid box 196 202 REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level 197 203 REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level 204 REAL(wp) :: vsws_int !< surface momentum flux (u component) at particle grid box 198 205 REAL(wp) :: vv_int !< 199 206 REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level … … 258 265 j = jp + block_offset(nb)%j_off 259 266 k = kp + block_offset(nb)%k_off 260 261 262 267 ! 263 268 !-- Interpolate u velocity-component … … 271 276 !-- Monin-Obukhov relations (if branch). 272 277 !-- First, check if particle is located below first vertical grid level 273 !-- (Prandtl-layer height)278 !-- above topography (Prandtl-layer height) 274 279 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 275 280 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 276 277 IF ( constant_flux_layer .AND. & 278 zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 281 ! 282 !-- Determine vertical index of topography top 283 k_wall = MAXLOC( & 284 MERGE( 1, 0, & 285 BTEST( wall_flags_0(nzb:nzb_max,jlog,ilog), 12 ) & 286 ), DIM = 1 & 287 ) - 1 288 289 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 279 290 ! 280 291 !-- Resolved-scale horizontal particle velocity is zero below z0. 281 IF ( zv(n) - zw( nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN292 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 282 293 u_int(n) = 0.0_wp 283 294 ELSE 284 295 ! 285 296 !-- Determine the sublayer. Further used as index. 286 height_p = ( zv(n) - zw( nzb_s_inner(jlog,ilog)) - z0_av_global ) &297 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 287 298 * REAL( number_of_sublayers, KIND=wp ) & 288 299 * d_z_p_z0 … … 296 307 ) 297 308 ! 298 !-- Limit friction velocity. In narrow canyons or holes the 299 !-- friction velocity can become very small, resulting in a too 300 !-- large particle speed. 301 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ), & 302 0.01_wp ) 309 !-- Get friction velocity and momentum flux from new surface data 310 !-- types. 311 IF ( surf_def_h(0)%start_index(jlog,ilog) <= & 312 surf_def_h(0)%end_index(jlog,ilog) ) THEN 313 surf_start = surf_def_h(0)%start_index(jlog,ilog) 314 !-- Limit friction velocity. In narrow canyons or holes the 315 !-- friction velocity can become very small, resulting in a too 316 !-- large particle speed. 317 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 318 usws_int = surf_def_h(0)%usws(surf_start) 319 ELSEIF ( surf_lsm_h%start_index(jlog,ilog) <= & 320 surf_lsm_h%end_index(jlog,ilog) ) THEN 321 surf_start = surf_lsm_h%start_index(jlog,ilog) 322 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 323 usws_int = surf_lsm_h%usws(surf_start) 324 ELSEIF ( surf_usm_h%start_index(jlog,ilog) <= & 325 surf_usm_h%end_index(jlog,ilog) ) THEN 326 surf_start = surf_usm_h%start_index(jlog,ilog) 327 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 328 usws_int = surf_usm_h%usws(surf_start) 329 ENDIF 330 303 331 ! 304 332 !-- Neutral solution is applied for all situations, e.g. also for … … 308 336 !-- as sensitivity studies revealed no significant effect of 309 337 !-- using the neutral solution also for un/stable situations. 310 u_int(n) = -usws (jlog,ilog) / ( us_int * kappa + 1E-10_wp )&338 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 311 339 * log_z_z0_int - u_gtrans 312 340 … … 352 380 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 353 381 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 354 IF ( constant_flux_layer .AND. & 355 zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 356 357 IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN 382 ! 383 !-- Determine vertical index of topography top 384 k_wall = MAXLOC( & 385 MERGE( 1, 0, & 386 BTEST( wall_flags_0(nzb:nzb_max,jlog,ilog), 12 ) & 387 ), DIM = 1 & 388 ) - 1 389 390 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 391 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 358 392 ! 359 393 !-- Resolved-scale horizontal particle velocity is zero below z0. … … 365 399 !-- topography particle on u-grid can be above surface-layer height, 366 400 !-- whereas it can be below on v-grid. 367 height_p = ( zv(n) - zw( nzb_s_inner(jlog,ilog)) - z0_av_global ) &401 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 368 402 * REAL( number_of_sublayers, KIND=wp ) & 369 403 * d_z_p_z0 … … 377 411 ) 378 412 ! 379 !-- Limit friction velocity. In narrow canyons or holes the 380 !-- friction velocity can become very small, resulting in a too 381 !-- large particle speed. 382 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), & 383 0.01_wp ) 413 !-- Get friction velocity and momentum flux from new surface data 414 !-- types. 415 IF ( surf_def_h(0)%start_index(jlog,ilog) <= & 416 surf_def_h(0)%end_index(jlog,ilog) ) THEN 417 surf_start = surf_def_h(0)%start_index(jlog,ilog) 418 !-- Limit friction velocity. In narrow canyons or holes the 419 !-- friction velocity can become very small, resulting in a too 420 !-- large particle speed. 421 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 422 vsws_int = surf_def_h(0)%usws(surf_start) 423 ELSEIF ( surf_lsm_h%start_index(jlog,ilog) <= & 424 surf_lsm_h%end_index(jlog,ilog) ) THEN 425 surf_start = surf_lsm_h%start_index(jlog,ilog) 426 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 427 vsws_int = surf_lsm_h%usws(surf_start) 428 ELSEIF ( surf_usm_h%start_index(jlog,ilog) <= & 429 surf_usm_h%end_index(jlog,ilog) ) THEN 430 surf_start = surf_usm_h%start_index(jlog,ilog) 431 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 432 vsws_int = surf_usm_h%usws(surf_start) 433 ENDIF 384 434 ! 385 435 !-- Neutral solution is applied for all situations, e.g. also for … … 389 439 !-- as sensitivity studies revealed no significant effect of 390 440 !-- using the neutral solution also for un/stable situations. 391 v_int(n) = -vsws (jlog,ilog) / ( us_int * kappa + 1E-10_wp )&441 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 392 442 * log_z_z0_int - v_gtrans 393 443 … … 622 672 num_gp = 0 623 673 624 IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 674 ! 675 !-- Determine vertical index of topography top at (j,i) 676 k_wall = MAXLOC( & 677 MERGE( 1, 0, & 678 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) & 679 ), DIM = 1 & 680 ) - 1 681 ! 682 !-- To do: Reconsider order of computations in order to avoid 683 !-- unnecessary index calculations. 684 IF ( k > k_wall .OR. k_wall == 0 ) THEN 625 685 num_gp = num_gp + 1 626 686 gp_outside_of_building(1) = 1 … … 634 694 de_dzi(num_gp) = de_dz(k,j,i) 635 695 ENDIF 636 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 696 697 ! 698 !-- Determine vertical index of topography top at (j+1,i) 699 k_wall = MAXLOC( & 700 MERGE( 1, 0, & 701 BTEST( wall_flags_0(nzb:nzb_max,j+1,i), 12 ) & 702 ), DIM = 1 & 703 ) - 1 704 IF ( k > k_wall .OR. k_wall == 0 ) THEN 637 705 num_gp = num_gp + 1 638 706 gp_outside_of_building(2) = 1 … … 647 715 ENDIF 648 716 649 IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 717 ! 718 !-- Determine vertical index of topography top at (j,i) 719 k_wall = MAXLOC( & 720 MERGE( 1, 0, & 721 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) & 722 ), DIM = 1 & 723 ) - 1 724 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 650 725 num_gp = num_gp + 1 651 726 gp_outside_of_building(3) = 1 … … 660 735 ENDIF 661 736 662 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 737 ! 738 !-- Determine vertical index of topography top at (j+1,i) 739 k_wall = MAXLOC( & 740 MERGE( 1, 0, & 741 BTEST( wall_flags_0(nzb:nzb_max,j+1,i), 12 ) & 742 ), DIM = 1 & 743 ) - 1 744 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 663 745 num_gp = num_gp + 1 664 746 gp_outside_of_building(4) = 1 … … 673 755 ENDIF 674 756 675 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 757 ! 758 !-- Determine vertical index of topography top at (j,i+1) 759 k_wall = MAXLOC( & 760 MERGE( 1, 0, & 761 BTEST( wall_flags_0(nzb:nzb_max,j,i+1), 12 ) & 762 ), DIM = 1 & 763 ) - 1 764 IF ( k > k_wall .OR. k_wall == 0 ) THEN 676 765 num_gp = num_gp + 1 677 766 gp_outside_of_building(5) = 1 … … 686 775 ENDIF 687 776 688 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) THEN 777 ! 778 !-- Determine vertical index of topography top at (j+1,i+1) 779 k_wall = MAXLOC( & 780 MERGE( 1, 0, & 781 BTEST( wall_flags_0(nzb:nzb_max,j+1,i+1), 12 )& 782 ), DIM = 1 & 783 ) - 1 784 IF ( k > k_wall .OR. k_wall == 0 ) THEN 689 785 num_gp = num_gp + 1 690 786 gp_outside_of_building(6) = 1 … … 699 795 ENDIF 700 796 701 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 797 ! 798 !-- Determine vertical index of topography top at (j,i+1) 799 k_wall = MAXLOC( & 800 MERGE( 1, 0, & 801 BTEST( wall_flags_0(nzb:nzb_max,j,i+1), 12 ) & 802 ), DIM = 1 & 803 ) - 1 804 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 702 805 num_gp = num_gp + 1 703 806 gp_outside_of_building(7) = 1 … … 712 815 ENDIF 713 816 714 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0) THEN 817 ! 818 !-- Determine vertical index of topography top at (j+1,i+1) 819 k_wall = MAXLOC( & 820 MERGE( 1, 0, & 821 BTEST( wall_flags_0(nzb:nzb_max,j+1,i+1), 12 )& 822 ), DIM = 1 & 823 ) - 1 824 IF ( k+1 > k_wall .OR. k_wall == 0) THEN 715 825 num_gp = num_gp + 1 716 826 gp_outside_of_building(8) = 1
Note: See TracChangeset
for help on using the changeset viewer.
