Changeset 2417 for palm/trunk/SOURCE
- Timestamp:
- Sep 6, 2017 3:22:27 PM (7 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
lpm.f90 (modified) (7 diffs)
-
lpm_advec.f90 (modified) (11 diffs)
-
lpm_pack_arrays.f90 (modified) (5 diffs)
-
user_lpm_advec.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm.f90
r2263 r2417 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Major bugfixes in modeling SGS particle speeds (since revision 1359). 28 ! Particle sorting added to distinguish between already completed and 29 ! non-completed particles. 30 ! 31 ! 2263 2017-06-08 14:59:01Z schwenkel 27 32 ! Implemented splitting and merging algorithm 28 33 ! … … 140 145 141 146 USE lpm_pack_arrays_mod, & 142 ONLY: lpm_pack_all_arrays 147 ONLY: lpm_pack_all_arrays, lpm_sort 143 148 144 149 USE particle_attributes, & … … 259 264 CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' ) 260 265 261 grid_particles(:,:,:)%time_loop_done = .TRUE.262 266 ! 263 267 !-- If particle advection includes SGS velocity components, calculate the … … 265 269 !-- horizontally averaged profiles of the SGS TKE and the resolved-scale 266 270 !-- velocity variances) 267 268 271 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 269 272 CALL lpm_init_sgs_tke 270 273 ENDIF 274 275 ! 276 !-- In case SGS-particle speed is considered, particles may carry out 277 !-- several particle timesteps. In order to prevent unnecessary 278 !-- treatment of particles that already reached the final time level, 279 !-- particles are sorted into contiguous blocks of finished and 280 !-- not-finished particles, in addition to their already sorting 281 !-- according to their sub-boxes. 282 IF ( .NOT. first_loop_stride .AND. use_sgs_for_particles ) & 283 CALL lpm_sort 271 284 272 285 DO i = nxl, nxr … … 323 336 !-- optimization is still possible.) 324 337 IF ( topography /= 'flat' .AND. k < nzb_max + 2 ) THEN 325 CALL lpm_boundary_conds( 'walls' ) 338 CALL lpm_boundary_conds( 'walls' )!, i, j, k ) 326 339 ENDIF 327 340 ! … … 333 346 !-- the top or bottom boundary and delete those particles, which are 334 347 !-- older than allowed 335 CALL lpm_boundary_conds( 'bottom/top' ) 348 CALL lpm_boundary_conds( 'bottom/top' ) !, i, j, k ) 336 349 ! 337 350 !--- If not all particles of the actual grid cell have reached the 338 !-- LES timestep, this cell has to to another loop iteration. Due to 339 !-- the fact that particles can move into neighboring grid cell, 340 !-- these neighbor cells also have to perform another loop iteration 351 !-- LES timestep, this cell has to do another loop iteration. Due to 352 !-- the fact that particles can move into neighboring grid cells, 353 !-- these neighbor cells also have to perform another loop iteration. 354 !-- Please note, this realization does not work properly if 355 !-- particles move into another subdomain. 341 356 IF ( .NOT. dt_3d_reached_l ) THEN 342 ks = MAX(nzb+1,k )343 ke = MIN(nzt,k )344 js = MAX(nys,j )345 je = MIN(nyn,j )346 is = MAX(nxl,i )347 ie = MIN(nxr,i )357 ks = MAX(nzb+1,k-1) 358 ke = MIN(nzt,k+1) 359 js = MAX(nys,j-1) 360 je = MIN(nyn,j+1) 361 is = MAX(nxl,i-1) 362 ie = MIN(nxr,i+1) 348 363 grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE. 364 ELSE 365 grid_particles(k,j,i)%time_loop_done = .TRUE. 349 366 ENDIF 350 367 … … 392 409 !-- determine new number of particles 393 410 CALL lpm_pack_all_arrays 394 395 411 ! 396 412 !-- Initialize variables for the next (sub-) timestep, i.e., for marking -
palm/trunk/SOURCE/lpm_advec.f90
r2318 r2417 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Particle loops adapted for sub-box structure, i.e. for each sub-box the 28 ! particle loop runs from start_index up to end_index instead from 1 to 29 ! number_of_particles. This way, it is possible to skip unnecessary 30 ! computations for particles that already completed the LES timestep. 31 ! 32 ! 2318 2017-07-20 17:27:44Z suehring 27 33 ! Get topography top index via Function call 28 34 ! … … 267 273 268 274 DO nb = 0, 7 269 275 270 276 i = ip 271 277 j = jp + block_offset(nb)%j_off … … 287 293 ! 288 294 !-- Determine vertical index of topography top 289 k_wall = get_topography_top_index( jlog, ilog, 's' )295 k_wall = get_topography_top_index( jlog, ilog, 's' ) 290 296 291 297 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 648 654 ELSE ! non-flat topography, e.g., buildings 649 655 650 DO n = 1, number_of_particles 651 i = particles(n)%x * ddx 652 j = particles(n)%y * ddy 653 k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz & 654 + offset_ocean_nzt ! only exact if eq.dist 655 ! 656 !-- In case that there are buildings it has to be determined 657 !-- how many of the gridpoints defining the particle box are 658 !-- situated within a building 659 !-- gp_outside_of_building(1): i,j,k 660 !-- gp_outside_of_building(2): i,j+1,k 661 !-- gp_outside_of_building(3): i,j,k+1 662 !-- gp_outside_of_building(4): i,j+1,k+1 663 !-- gp_outside_of_building(5): i+1,j,k 664 !-- gp_outside_of_building(6): i+1,j+1,k 665 !-- gp_outside_of_building(7): i+1,j,k+1 666 !-- gp_outside_of_building(8): i+1,j+1,k+1 667 668 gp_outside_of_building = 0 669 location = 0.0_wp 670 num_gp = 0 671 672 ! 673 !-- Determine vertical index of topography top at (j,i) 674 k_wall = get_topography_top_index( j, i, 's' ) 675 ! 676 !-- To do: Reconsider order of computations in order to avoid 677 !-- unnecessary index calculations. 678 IF ( k > k_wall .OR. k_wall == 0 ) THEN 679 num_gp = num_gp + 1 680 gp_outside_of_building(1) = 1 681 location(num_gp,1) = i * dx 682 location(num_gp,2) = j * dy 683 location(num_gp,3) = k * dz - 0.5_wp * dz 684 ei(num_gp) = e(k,j,i) 685 dissi(num_gp) = diss(k,j,i) 686 de_dxi(num_gp) = de_dx(k,j,i) 687 de_dyi(num_gp) = de_dy(k,j,i) 688 de_dzi(num_gp) = de_dz(k,j,i) 689 ENDIF 690 691 ! 692 !-- Determine vertical index of topography top at (j+1,i) 693 k_wall = get_topography_top_index( j+1, i, 's' ) 694 695 IF ( k > k_wall .OR. k_wall == 0 ) THEN 696 num_gp = num_gp + 1 697 gp_outside_of_building(2) = 1 698 location(num_gp,1) = i * dx 699 location(num_gp,2) = (j+1) * dy 700 location(num_gp,3) = k * dz - 0.5_wp * dz 701 ei(num_gp) = e(k,j+1,i) 702 dissi(num_gp) = diss(k,j+1,i) 703 de_dxi(num_gp) = de_dx(k,j+1,i) 704 de_dyi(num_gp) = de_dy(k,j+1,i) 705 de_dzi(num_gp) = de_dz(k,j+1,i) 706 ENDIF 707 708 ! 709 !-- Determine vertical index of topography top at (j,i) 710 k_wall = get_topography_top_index( j, i, 's' ) 711 712 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 713 num_gp = num_gp + 1 714 gp_outside_of_building(3) = 1 715 location(num_gp,1) = i * dx 716 location(num_gp,2) = j * dy 717 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 718 ei(num_gp) = e(k+1,j,i) 719 dissi(num_gp) = diss(k+1,j,i) 720 de_dxi(num_gp) = de_dx(k+1,j,i) 721 de_dyi(num_gp) = de_dy(k+1,j,i) 722 de_dzi(num_gp) = de_dz(k+1,j,i) 723 ENDIF 724 725 ! 726 !-- Determine vertical index of topography top at (j+1,i) 727 k_wall = get_topography_top_index( j+1, i, 's' ) 728 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 729 num_gp = num_gp + 1 730 gp_outside_of_building(4) = 1 731 location(num_gp,1) = i * dx 732 location(num_gp,2) = (j+1) * dy 733 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 734 ei(num_gp) = e(k+1,j+1,i) 735 dissi(num_gp) = diss(k+1,j+1,i) 736 de_dxi(num_gp) = de_dx(k+1,j+1,i) 737 de_dyi(num_gp) = de_dy(k+1,j+1,i) 738 de_dzi(num_gp) = de_dz(k+1,j+1,i) 739 ENDIF 740 741 ! 742 !-- Determine vertical index of topography top at (j,i+1) 743 k_wall = get_topography_top_index( j, i+1, 's' ) 744 IF ( k > k_wall .OR. k_wall == 0 ) THEN 745 num_gp = num_gp + 1 746 gp_outside_of_building(5) = 1 747 location(num_gp,1) = (i+1) * dx 748 location(num_gp,2) = j * dy 749 location(num_gp,3) = k * dz - 0.5_wp * dz 750 ei(num_gp) = e(k,j,i+1) 751 dissi(num_gp) = diss(k,j,i+1) 752 de_dxi(num_gp) = de_dx(k,j,i+1) 753 de_dyi(num_gp) = de_dy(k,j,i+1) 754 de_dzi(num_gp) = de_dz(k,j,i+1) 755 ENDIF 756 757 ! 758 !-- Determine vertical index of topography top at (j+1,i+1) 759 k_wall = get_topography_top_index( j+1, i+1, 's' ) 760 761 IF ( k > k_wall .OR. k_wall == 0 ) THEN 762 num_gp = num_gp + 1 763 gp_outside_of_building(6) = 1 764 location(num_gp,1) = (i+1) * dx 765 location(num_gp,2) = (j+1) * dy 766 location(num_gp,3) = k * dz - 0.5_wp * dz 767 ei(num_gp) = e(k,j+1,i+1) 768 dissi(num_gp) = diss(k,j+1,i+1) 769 de_dxi(num_gp) = de_dx(k,j+1,i+1) 770 de_dyi(num_gp) = de_dy(k,j+1,i+1) 771 de_dzi(num_gp) = de_dz(k,j+1,i+1) 772 ENDIF 773 774 ! 775 !-- Determine vertical index of topography top at (j,i+1) 776 k_wall = get_topography_top_index( j, i+1, 's' ) 777 778 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 779 num_gp = num_gp + 1 780 gp_outside_of_building(7) = 1 781 location(num_gp,1) = (i+1) * dx 782 location(num_gp,2) = j * dy 783 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 784 ei(num_gp) = e(k+1,j,i+1) 785 dissi(num_gp) = diss(k+1,j,i+1) 786 de_dxi(num_gp) = de_dx(k+1,j,i+1) 787 de_dyi(num_gp) = de_dy(k+1,j,i+1) 788 de_dzi(num_gp) = de_dz(k+1,j,i+1) 789 ENDIF 790 791 ! 792 !-- Determine vertical index of topography top at (j+1,i+1) 793 k_wall = get_topography_top_index( j+1, i+1, 's' ) 794 795 IF ( k+1 > k_wall .OR. k_wall == 0) THEN 796 num_gp = num_gp + 1 797 gp_outside_of_building(8) = 1 798 location(num_gp,1) = (i+1) * dx 799 location(num_gp,2) = (j+1) * dy 800 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 801 ei(num_gp) = e(k+1,j+1,i+1) 802 dissi(num_gp) = diss(k+1,j+1,i+1) 803 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 804 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 805 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 806 ENDIF 807 ! 808 !-- If all gridpoints are situated outside of a building, then the 809 !-- ordinary interpolation scheme can be used. 810 IF ( num_gp == 8 ) THEN 811 812 x = particles(n)%x - i * dx 813 y = particles(n)%y - j * dy 814 aa = x**2 + y**2 815 bb = ( dx - x )**2 + y**2 816 cc = x**2 + ( dy - y )**2 817 dd = ( dx - x )**2 + ( dy - y )**2 818 gg = aa + bb + cc + dd 819 820 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 821 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 822 ) / ( 3.0_wp * gg ) 823 824 IF ( k == nzt ) THEN 825 e_int(n) = e_int_l 826 ELSE 827 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 828 ( gg - bb ) * e(k+1,j,i+1) + & 829 ( gg - cc ) * e(k+1,j+1,i) + & 830 ( gg - dd ) * e(k+1,j+1,i+1) & 831 ) / ( 3.0_wp * gg ) 832 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 833 ( e_int_u - e_int_l ) 834 ENDIF 835 ! 836 !-- Needed to avoid NaN particle velocities (this might not be 837 !-- required any more) 838 IF ( e_int(n) <= 0.0_wp ) THEN 839 e_int(n) = 1.0E-20_wp 840 ENDIF 841 ! 842 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 843 !-- and all position variables from above (TKE)) 844 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 845 ( gg - bb ) * de_dx(k,j,i+1) + & 846 ( gg - cc ) * de_dx(k,j+1,i) + & 847 ( gg - dd ) * de_dx(k,j+1,i+1) & 848 ) / ( 3.0_wp * gg ) 849 850 IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN 851 de_dx_int(n) = de_dx_int_l 852 ELSE 853 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 854 ( gg - bb ) * de_dx(k+1,j,i+1) + & 855 ( gg - cc ) * de_dx(k+1,j+1,i) + & 856 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 857 ) / ( 3.0_wp * gg ) 858 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / & 859 dz * ( de_dx_int_u - de_dx_int_l ) 860 ENDIF 861 862 ! 863 !-- Interpolate the TKE gradient along y 864 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 865 ( gg - bb ) * de_dy(k,j,i+1) + & 866 ( gg - cc ) * de_dy(k,j+1,i) + & 867 ( gg - dd ) * de_dy(k,j+1,i+1) & 868 ) / ( 3.0_wp * gg ) 869 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 870 de_dy_int(n) = de_dy_int_l 871 ELSE 872 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 873 ( gg - bb ) * de_dy(k+1,j,i+1) + & 874 ( gg - cc ) * de_dy(k+1,j+1,i) + & 875 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 876 ) / ( 3.0_wp * gg ) 877 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / & 878 dz * ( de_dy_int_u - de_dy_int_l ) 879 ENDIF 880 881 ! 882 !-- Interpolate the TKE gradient along z 883 IF ( zv(n) < 0.5_wp * dz ) THEN 884 de_dz_int(n) = 0.0_wp 885 ELSE 886 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 887 ( gg - bb ) * de_dz(k,j,i+1) + & 888 ( gg - cc ) * de_dz(k,j+1,i) + & 889 ( gg - dd ) * de_dz(k,j+1,i+1) & 890 ) / ( 3.0_wp * gg ) 891 892 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 893 de_dz_int(n) = de_dz_int_l 894 ELSE 895 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 896 ( gg - bb ) * de_dz(k+1,j,i+1) + & 897 ( gg - cc ) * de_dz(k+1,j+1,i) + & 898 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 899 ) / ( 3.0_wp * gg ) 900 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /& 901 dz * ( de_dz_int_u - de_dz_int_l ) 902 ENDIF 903 ENDIF 904 905 ! 906 !-- Interpolate the dissipation of TKE 907 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 908 ( gg - bb ) * diss(k,j,i+1) + & 909 ( gg - cc ) * diss(k,j+1,i) + & 910 ( gg - dd ) * diss(k,j+1,i+1) & 911 ) / ( 3.0_wp * gg ) 912 913 IF ( k == nzt ) THEN 914 diss_int(n) = diss_int_l 915 ELSE 916 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 917 ( gg - bb ) * diss(k+1,j,i+1) + & 918 ( gg - cc ) * diss(k+1,j+1,i) + & 919 ( gg - dd ) * diss(k+1,j+1,i+1) & 920 ) / ( 3.0_wp * gg ) 921 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *& 922 ( diss_int_u - diss_int_l ) 923 ENDIF 924 ! 925 !-- Set flag for stochastic equation. 926 term_1_2(n) = 1.0_wp 927 928 ELSE 929 930 ! 931 !-- If wall between gridpoint 1 and gridpoint 5, then 932 !-- Neumann boundary condition has to be applied 933 IF ( gp_outside_of_building(1) == 1 .AND. & 934 gp_outside_of_building(5) == 0 ) THEN 656 DO nb = 0, 7 657 658 DO n = start_index(nb), end_index(nb) 659 660 i = particles(n)%x * ddx 661 j = particles(n)%y * ddy 662 k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz & 663 + offset_ocean_nzt ! only exact if eq.dist 664 ! 665 !-- In case that there are buildings it has to be determined 666 !-- how many of the gridpoints defining the particle box are 667 !-- situated within a building 668 !-- gp_outside_of_building(1): i,j,k 669 !-- gp_outside_of_building(2): i,j+1,k 670 !-- gp_outside_of_building(3): i,j,k+1 671 !-- gp_outside_of_building(4): i,j+1,k+1 672 !-- gp_outside_of_building(5): i+1,j,k 673 !-- gp_outside_of_building(6): i+1,j+1,k 674 !-- gp_outside_of_building(7): i+1,j,k+1 675 !-- gp_outside_of_building(8): i+1,j+1,k+1 676 677 gp_outside_of_building = 0 678 location = 0.0_wp 679 num_gp = 0 680 681 ! 682 !-- Determine vertical index of topography top at (j,i) 683 k_wall = get_topography_top_index( j, i, 's' ) 684 ! 685 !-- To do: Reconsider order of computations in order to avoid 686 !-- unnecessary index calculations. 687 IF ( k > k_wall .OR. k_wall == 0 ) THEN 935 688 num_gp = num_gp + 1 936 location(num_gp,1) = i * dx + 0.5_wp * dx 689 gp_outside_of_building(1) = 1 690 location(num_gp,1) = i * dx 937 691 location(num_gp,2) = j * dy 938 692 location(num_gp,3) = k * dz - 0.5_wp * dz 939 693 ei(num_gp) = e(k,j,i) 940 694 dissi(num_gp) = diss(k,j,i) 941 de_dxi(num_gp) = 0.0_wp695 de_dxi(num_gp) = de_dx(k,j,i) 942 696 de_dyi(num_gp) = de_dy(k,j,i) 943 697 de_dzi(num_gp) = de_dz(k,j,i) 944 698 ENDIF 945 699 946 IF ( gp_outside_of_building(5) == 1 .AND. & 947 gp_outside_of_building(1) == 0 ) THEN 700 ! 701 !-- Determine vertical index of topography top at (j+1,i) 702 k_wall = get_topography_top_index( j+1, i, 's' ) 703 704 IF ( k > k_wall .OR. k_wall == 0 ) THEN 948 705 num_gp = num_gp + 1 949 location(num_gp,1) = i * dx + 0.5_wp * dx 706 gp_outside_of_building(2) = 1 707 location(num_gp,1) = i * dx 708 location(num_gp,2) = (j+1) * dy 709 location(num_gp,3) = k * dz - 0.5_wp * dz 710 ei(num_gp) = e(k,j+1,i) 711 dissi(num_gp) = diss(k,j+1,i) 712 de_dxi(num_gp) = de_dx(k,j+1,i) 713 de_dyi(num_gp) = de_dy(k,j+1,i) 714 de_dzi(num_gp) = de_dz(k,j+1,i) 715 ENDIF 716 717 ! 718 !-- Determine vertical index of topography top at (j,i) 719 k_wall = get_topography_top_index( j, i, 's' ) 720 721 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 722 num_gp = num_gp + 1 723 gp_outside_of_building(3) = 1 724 location(num_gp,1) = i * dx 725 location(num_gp,2) = j * dy 726 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 727 ei(num_gp) = e(k+1,j,i) 728 dissi(num_gp) = diss(k+1,j,i) 729 de_dxi(num_gp) = de_dx(k+1,j,i) 730 de_dyi(num_gp) = de_dy(k+1,j,i) 731 de_dzi(num_gp) = de_dz(k+1,j,i) 732 ENDIF 733 734 ! 735 !-- Determine vertical index of topography top at (j+1,i) 736 k_wall = get_topography_top_index( j+1, i, 's' ) 737 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 738 num_gp = num_gp + 1 739 gp_outside_of_building(4) = 1 740 location(num_gp,1) = i * dx 741 location(num_gp,2) = (j+1) * dy 742 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 743 ei(num_gp) = e(k+1,j+1,i) 744 dissi(num_gp) = diss(k+1,j+1,i) 745 de_dxi(num_gp) = de_dx(k+1,j+1,i) 746 de_dyi(num_gp) = de_dy(k+1,j+1,i) 747 de_dzi(num_gp) = de_dz(k+1,j+1,i) 748 ENDIF 749 750 ! 751 !-- Determine vertical index of topography top at (j,i+1) 752 k_wall = get_topography_top_index( j, i+1, 's' ) 753 IF ( k > k_wall .OR. k_wall == 0 ) THEN 754 num_gp = num_gp + 1 755 gp_outside_of_building(5) = 1 756 location(num_gp,1) = (i+1) * dx 950 757 location(num_gp,2) = j * dy 951 758 location(num_gp,3) = k * dz - 0.5_wp * dz 952 759 ei(num_gp) = e(k,j,i+1) 953 760 dissi(num_gp) = diss(k,j,i+1) 954 de_dxi(num_gp) = 0.0_wp761 de_dxi(num_gp) = de_dx(k,j,i+1) 955 762 de_dyi(num_gp) = de_dy(k,j,i+1) 956 763 de_dzi(num_gp) = de_dz(k,j,i+1) … … 958 765 959 766 ! 960 !-- If wall between gridpoint 5 and gridpoint 6, then961 !-- then Neumann boundary condition has to be applied 962 IF ( gp_outside_of_building(5) == 1 .AND. & 963 gp_outside_of_building(6)== 0 ) THEN767 !-- Determine vertical index of topography top at (j+1,i+1) 768 k_wall = get_topography_top_index( j+1, i+1, 's' ) 769 770 IF ( k > k_wall .OR. k_wall == 0 ) THEN 964 771 num_gp = num_gp + 1 772 gp_outside_of_building(6) = 1 965 773 location(num_gp,1) = (i+1) * dx 966 location(num_gp,2) = j * dy + 0.5_wp * dy 967 location(num_gp,3) = k * dz - 0.5_wp * dz 968 ei(num_gp) = e(k,j,i+1) 969 dissi(num_gp) = diss(k,j,i+1) 970 de_dxi(num_gp) = de_dx(k,j,i+1) 971 de_dyi(num_gp) = 0.0_wp 972 de_dzi(num_gp) = de_dz(k,j,i+1) 973 ENDIF 974 975 IF ( gp_outside_of_building(6) == 1 .AND. & 976 gp_outside_of_building(5) == 0 ) THEN 977 num_gp = num_gp + 1 978 location(num_gp,1) = (i+1) * dx 979 location(num_gp,2) = j * dy + 0.5_wp * dy 774 location(num_gp,2) = (j+1) * dy 980 775 location(num_gp,3) = k * dz - 0.5_wp * dz 981 776 ei(num_gp) = e(k,j+1,i+1) 982 777 dissi(num_gp) = diss(k,j+1,i+1) 983 778 de_dxi(num_gp) = de_dx(k,j+1,i+1) 984 de_dyi(num_gp) = 0.0_wp985 de_dzi(num_gp) = de_dz(k,j+1,i+1)986 ENDIF987 988 !989 !-- If wall between gridpoint 2 and gridpoint 6, then990 !-- Neumann boundary condition has to be applied991 IF ( gp_outside_of_building(2) == 1 .AND. &992 gp_outside_of_building(6) == 0 ) THEN993 num_gp = num_gp + 1994 location(num_gp,1) = i * dx + 0.5_wp * dx995 location(num_gp,2) = (j+1) * dy996 location(num_gp,3) = k * dz - 0.5_wp * dz997 ei(num_gp) = e(k,j+1,i)998 dissi(num_gp) = diss(k,j+1,i)999 de_dxi(num_gp) = 0.0_wp1000 de_dyi(num_gp) = de_dy(k,j+1,i)1001 de_dzi(num_gp) = de_dz(k,j+1,i)1002 ENDIF1003 1004 IF ( gp_outside_of_building(6) == 1 .AND. &1005 gp_outside_of_building(2) == 0 ) THEN1006 num_gp = num_gp + 11007 location(num_gp,1) = i * dx + 0.5_wp * dx1008 location(num_gp,2) = (j+1) * dy1009 location(num_gp,3) = k * dz - 0.5_wp * dz1010 ei(num_gp) = e(k,j+1,i+1)1011 dissi(num_gp) = diss(k,j+1,i+1)1012 de_dxi(num_gp) = 0.0_wp1013 779 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1014 780 de_dzi(num_gp) = de_dz(k,j+1,i+1) … … 1016 782 1017 783 ! 1018 !-- If wall between gridpoint 1 and gridpoint 2, then1019 !-- Neumann boundary condition has to be applied 1020 IF ( gp_outside_of_building(1) == 1 .AND. & 1021 gp_outside_of_building(2)== 0 ) THEN784 !-- Determine vertical index of topography top at (j,i+1) 785 k_wall = get_topography_top_index( j, i+1, 's' ) 786 787 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 1022 788 num_gp = num_gp + 1 1023 location(num_gp,1) = i * dx 1024 location(num_gp,2) = j * dy + 0.5_wp * dy 1025 location(num_gp,3) = k * dz - 0.5_wp * dz 1026 ei(num_gp) = e(k,j,i) 1027 dissi(num_gp) = diss(k,j,i) 1028 de_dxi(num_gp) = de_dx(k,j,i) 1029 de_dyi(num_gp) = 0.0_wp 1030 de_dzi(num_gp) = de_dz(k,j,i) 1031 ENDIF 1032 1033 IF ( gp_outside_of_building(2) == 1 .AND. & 1034 gp_outside_of_building(1) == 0 ) THEN 1035 num_gp = num_gp + 1 1036 location(num_gp,1) = i * dx 1037 location(num_gp,2) = j * dy + 0.5_wp * dy 1038 location(num_gp,3) = k * dz - 0.5_wp * dz 1039 ei(num_gp) = e(k,j+1,i) 1040 dissi(num_gp) = diss(k,j+1,i) 1041 de_dxi(num_gp) = de_dx(k,j+1,i) 1042 de_dyi(num_gp) = 0.0_wp 1043 de_dzi(num_gp) = de_dz(k,j+1,i) 1044 ENDIF 1045 1046 ! 1047 !-- If wall between gridpoint 3 and gridpoint 7, then 1048 !-- Neumann boundary condition has to be applied 1049 IF ( gp_outside_of_building(3) == 1 .AND. & 1050 gp_outside_of_building(7) == 0 ) THEN 1051 num_gp = num_gp + 1 1052 location(num_gp,1) = i * dx + 0.5_wp * dx 1053 location(num_gp,2) = j * dy 1054 location(num_gp,3) = k * dz + 0.5_wp * dz 1055 ei(num_gp) = e(k+1,j,i) 1056 dissi(num_gp) = diss(k+1,j,i) 1057 de_dxi(num_gp) = 0.0_wp 1058 de_dyi(num_gp) = de_dy(k+1,j,i) 1059 de_dzi(num_gp) = de_dz(k+1,j,i) 1060 ENDIF 1061 1062 IF ( gp_outside_of_building(7) == 1 .AND. & 1063 gp_outside_of_building(3) == 0 ) THEN 1064 num_gp = num_gp + 1 1065 location(num_gp,1) = i * dx + 0.5_wp * dx 1066 location(num_gp,2) = j * dy 1067 location(num_gp,3) = k * dz + 0.5_wp * dz 1068 ei(num_gp) = e(k+1,j,i+1) 1069 dissi(num_gp) = diss(k+1,j,i+1) 1070 de_dxi(num_gp) = 0.0_wp 1071 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1072 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1073 ENDIF 1074 1075 ! 1076 !-- If wall between gridpoint 7 and gridpoint 8, then 1077 !-- Neumann boundary condition has to be applied 1078 IF ( gp_outside_of_building(7) == 1 .AND. & 1079 gp_outside_of_building(8) == 0 ) THEN 1080 num_gp = num_gp + 1 1081 location(num_gp,1) = (i+1) * dx 1082 location(num_gp,2) = j * dy + 0.5_wp * dy 1083 location(num_gp,3) = k * dz + 0.5_wp * dz 1084 ei(num_gp) = e(k+1,j,i+1) 1085 dissi(num_gp) = diss(k+1,j,i+1) 1086 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1087 de_dyi(num_gp) = 0.0_wp 1088 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1089 ENDIF 1090 1091 IF ( gp_outside_of_building(8) == 1 .AND. & 1092 gp_outside_of_building(7) == 0 ) THEN 1093 num_gp = num_gp + 1 1094 location(num_gp,1) = (i+1) * dx 1095 location(num_gp,2) = j * dy + 0.5_wp * dy 1096 location(num_gp,3) = k * dz + 0.5_wp * dz 1097 ei(num_gp) = e(k+1,j+1,i+1) 1098 dissi(num_gp) = diss(k+1,j+1,i+1) 1099 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1100 de_dyi(num_gp) = 0.0_wp 1101 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1102 ENDIF 1103 1104 ! 1105 !-- If wall between gridpoint 4 and gridpoint 8, then 1106 !-- Neumann boundary condition has to be applied 1107 IF ( gp_outside_of_building(4) == 1 .AND. & 1108 gp_outside_of_building(8) == 0 ) THEN 1109 num_gp = num_gp + 1 1110 location(num_gp,1) = i * dx + 0.5_wp * dx 1111 location(num_gp,2) = (j+1) * dy 1112 location(num_gp,3) = k * dz + 0.5_wp * dz 1113 ei(num_gp) = e(k+1,j+1,i) 1114 dissi(num_gp) = diss(k+1,j+1,i) 1115 de_dxi(num_gp) = 0.0_wp 1116 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1117 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1118 ENDIF 1119 1120 IF ( gp_outside_of_building(8) == 1 .AND. & 1121 gp_outside_of_building(4) == 0 ) THEN 1122 num_gp = num_gp + 1 1123 location(num_gp,1) = i * dx + 0.5_wp * dx 1124 location(num_gp,2) = (j+1) * dy 1125 location(num_gp,3) = k * dz + 0.5_wp * dz 1126 ei(num_gp) = e(k+1,j+1,i+1) 1127 dissi(num_gp) = diss(k+1,j+1,i+1) 1128 de_dxi(num_gp) = 0.0_wp 1129 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1130 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1131 ENDIF 1132 1133 ! 1134 !-- If wall between gridpoint 3 and gridpoint 4, then 1135 !-- Neumann boundary condition has to be applied 1136 IF ( gp_outside_of_building(3) == 1 .AND. & 1137 gp_outside_of_building(4) == 0 ) THEN 1138 num_gp = num_gp + 1 1139 location(num_gp,1) = i * dx 1140 location(num_gp,2) = j * dy + 0.5_wp * dy 1141 location(num_gp,3) = k * dz + 0.5_wp * dz 1142 ei(num_gp) = e(k+1,j,i) 1143 dissi(num_gp) = diss(k+1,j,i) 1144 de_dxi(num_gp) = de_dx(k+1,j,i) 1145 de_dyi(num_gp) = 0.0_wp 1146 de_dzi(num_gp) = de_dz(k+1,j,i) 1147 ENDIF 1148 1149 IF ( gp_outside_of_building(4) == 1 .AND. & 1150 gp_outside_of_building(3) == 0 ) THEN 1151 num_gp = num_gp + 1 1152 location(num_gp,1) = i * dx 1153 location(num_gp,2) = j * dy + 0.5_wp * dy 1154 location(num_gp,3) = k * dz + 0.5_wp * dz 1155 ei(num_gp) = e(k+1,j+1,i) 1156 dissi(num_gp) = diss(k+1,j+1,i) 1157 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1158 de_dyi(num_gp) = 0.0_wp 1159 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1160 ENDIF 1161 1162 ! 1163 !-- If wall between gridpoint 1 and gridpoint 3, then 1164 !-- Neumann boundary condition has to be applied 1165 !-- (only one case as only building beneath is possible) 1166 IF ( gp_outside_of_building(1) == 0 .AND. & 1167 gp_outside_of_building(3) == 1 ) THEN 1168 num_gp = num_gp + 1 1169 location(num_gp,1) = i * dx 1170 location(num_gp,2) = j * dy 1171 location(num_gp,3) = k * dz 1172 ei(num_gp) = e(k+1,j,i) 1173 dissi(num_gp) = diss(k+1,j,i) 1174 de_dxi(num_gp) = de_dx(k+1,j,i) 1175 de_dyi(num_gp) = de_dy(k+1,j,i) 1176 de_dzi(num_gp) = 0.0_wp 1177 ENDIF 1178 1179 ! 1180 !-- If wall between gridpoint 5 and gridpoint 7, then 1181 !-- Neumann boundary condition has to be applied 1182 !-- (only one case as only building beneath is possible) 1183 IF ( gp_outside_of_building(5) == 0 .AND. & 1184 gp_outside_of_building(7) == 1 ) THEN 1185 num_gp = num_gp + 1 789 gp_outside_of_building(7) = 1 1186 790 location(num_gp,1) = (i+1) * dx 1187 791 location(num_gp,2) = j * dy 1188 location(num_gp,3) = k* dz792 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 1189 793 ei(num_gp) = e(k+1,j,i+1) 1190 794 dissi(num_gp) = diss(k+1,j,i+1) 1191 795 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1192 796 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1193 de_dzi(num_gp) = 0.0_wp 1194 ENDIF 1195 1196 ! 1197 !-- If wall between gridpoint 2 and gridpoint 4, then 1198 !-- Neumann boundary condition has to be applied 1199 !-- (only one case as only building beneath is possible) 1200 IF ( gp_outside_of_building(2) == 0 .AND. & 1201 gp_outside_of_building(4) == 1 ) THEN 797 de_dzi(num_gp) = de_dz(k+1,j,i+1) 798 ENDIF 799 800 ! 801 !-- Determine vertical index of topography top at (j+1,i+1) 802 k_wall = get_topography_top_index( j+1, i+1, 's' ) 803 804 IF ( k+1 > k_wall .OR. k_wall == 0) THEN 1202 805 num_gp = num_gp + 1 1203 location(num_gp,1) = i * dx 1204 location(num_gp,2) = (j+1) * dy 1205 location(num_gp,3) = k * dz 1206 ei(num_gp) = e(k+1,j+1,i) 1207 dissi(num_gp) = diss(k+1,j+1,i) 1208 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1209 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1210 de_dzi(num_gp) = 0.0_wp 1211 ENDIF 1212 1213 ! 1214 !-- If wall between gridpoint 6 and gridpoint 8, then 1215 !-- Neumann boundary condition has to be applied 1216 !-- (only one case as only building beneath is possible) 1217 IF ( gp_outside_of_building(6) == 0 .AND. & 1218 gp_outside_of_building(8) == 1 ) THEN 1219 num_gp = num_gp + 1 806 gp_outside_of_building(8) = 1 1220 807 location(num_gp,1) = (i+1) * dx 1221 808 location(num_gp,2) = (j+1) * dy 1222 location(num_gp,3) = k* dz809 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 1223 810 ei(num_gp) = e(k+1,j+1,i+1) 1224 811 dissi(num_gp) = diss(k+1,j+1,i+1) 1225 812 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1226 813 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1227 de_dzi(num_gp) = 0.0_wp 1228 ENDIF 1229 1230 ! 1231 !-- Carry out the interpolation 1232 IF ( num_gp == 1 ) THEN 814 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 815 ENDIF 816 ! 817 !-- If all gridpoints are situated outside of a building, then the 818 !-- ordinary interpolation scheme can be used. 819 IF ( num_gp == 8 ) THEN 820 821 x = particles(n)%x - i * dx 822 y = particles(n)%y - j * dy 823 aa = x**2 + y**2 824 bb = ( dx - x )**2 + y**2 825 cc = x**2 + ( dy - y )**2 826 dd = ( dx - x )**2 + ( dy - y )**2 827 gg = aa + bb + cc + dd 828 829 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 830 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 831 ) / ( 3.0_wp * gg ) 832 833 IF ( k == nzt ) THEN 834 e_int(n) = e_int_l 835 ELSE 836 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 837 ( gg - bb ) * e(k+1,j,i+1) + & 838 ( gg - cc ) * e(k+1,j+1,i) + & 839 ( gg - dd ) * e(k+1,j+1,i+1) & 840 ) / ( 3.0_wp * gg ) 841 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 842 ( e_int_u - e_int_l ) 843 ENDIF 1233 844 ! 1234 !-- If only one of the gridpoints is situated outside of the 1235 !-- building, it follows that the values at the particle 1236 !-- location are the same as the gridpoint values 1237 e_int(n) = ei(num_gp) 1238 diss_int(n) = dissi(num_gp) 1239 de_dx_int(n) = de_dxi(num_gp) 1240 de_dy_int(n) = de_dyi(num_gp) 1241 de_dz_int(n) = de_dzi(num_gp) 845 !-- Needed to avoid NaN particle velocities (this might not be 846 !-- required any more) 847 IF ( e_int(n) <= 0.0_wp ) THEN 848 e_int(n) = 1.0E-20_wp 849 ENDIF 850 ! 851 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 852 !-- and all position variables from above (TKE)) 853 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 854 ( gg - bb ) * de_dx(k,j,i+1) + & 855 ( gg - cc ) * de_dx(k,j+1,i) + & 856 ( gg - dd ) * de_dx(k,j+1,i+1) & 857 ) / ( 3.0_wp * gg ) 858 859 IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN 860 de_dx_int(n) = de_dx_int_l 861 ELSE 862 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 863 ( gg - bb ) * de_dx(k+1,j,i+1) + & 864 ( gg - cc ) * de_dx(k+1,j+1,i) + & 865 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 866 ) / ( 3.0_wp * gg ) 867 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / & 868 dz * ( de_dx_int_u - de_dx_int_l ) 869 ENDIF 870 871 ! 872 !-- Interpolate the TKE gradient along y 873 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 874 ( gg - bb ) * de_dy(k,j,i+1) + & 875 ( gg - cc ) * de_dy(k,j+1,i) + & 876 ( gg - dd ) * de_dy(k,j+1,i+1) & 877 ) / ( 3.0_wp * gg ) 878 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 879 de_dy_int(n) = de_dy_int_l 880 ELSE 881 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 882 ( gg - bb ) * de_dy(k+1,j,i+1) + & 883 ( gg - cc ) * de_dy(k+1,j+1,i) + & 884 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 885 ) / ( 3.0_wp * gg ) 886 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / & 887 dz * ( de_dy_int_u - de_dy_int_l ) 888 ENDIF 889 890 ! 891 !-- Interpolate the TKE gradient along z 892 IF ( zv(n) < 0.5_wp * dz ) THEN 893 de_dz_int(n) = 0.0_wp 894 ELSE 895 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 896 ( gg - bb ) * de_dz(k,j,i+1) + & 897 ( gg - cc ) * de_dz(k,j+1,i) + & 898 ( gg - dd ) * de_dz(k,j+1,i+1) & 899 ) / ( 3.0_wp * gg ) 900 901 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 902 de_dz_int(n) = de_dz_int_l 903 ELSE 904 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 905 ( gg - bb ) * de_dz(k+1,j,i+1) + & 906 ( gg - cc ) * de_dz(k+1,j+1,i) + & 907 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 908 ) / ( 3.0_wp * gg ) 909 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /& 910 dz * ( de_dz_int_u - de_dz_int_l ) 911 ENDIF 912 ENDIF 913 914 ! 915 !-- Interpolate the dissipation of TKE 916 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 917 ( gg - bb ) * diss(k,j,i+1) + & 918 ( gg - cc ) * diss(k,j+1,i) + & 919 ( gg - dd ) * diss(k,j+1,i+1) & 920 ) / ( 3.0_wp * gg ) 921 922 IF ( k == nzt ) THEN 923 diss_int(n) = diss_int_l 924 ELSE 925 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 926 ( gg - bb ) * diss(k+1,j,i+1) + & 927 ( gg - cc ) * diss(k+1,j+1,i) + & 928 ( gg - dd ) * diss(k+1,j+1,i+1) & 929 ) / ( 3.0_wp * gg ) 930 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *& 931 ( diss_int_u - diss_int_l ) 932 ENDIF 933 ! 934 !-- Set flag for stochastic equation. 935 term_1_2(n) = 1.0_wp 936 937 ELSE 938 939 ! 940 !-- If wall between gridpoint 1 and gridpoint 5, then 941 !-- Neumann boundary condition has to be applied 942 IF ( gp_outside_of_building(1) == 1 .AND. & 943 gp_outside_of_building(5) == 0 ) THEN 944 num_gp = num_gp + 1 945 location(num_gp,1) = i * dx + 0.5_wp * dx 946 location(num_gp,2) = j * dy 947 location(num_gp,3) = k * dz - 0.5_wp * dz 948 ei(num_gp) = e(k,j,i) 949 dissi(num_gp) = diss(k,j,i) 950 de_dxi(num_gp) = 0.0_wp 951 de_dyi(num_gp) = de_dy(k,j,i) 952 de_dzi(num_gp) = de_dz(k,j,i) 953 ENDIF 954 955 IF ( gp_outside_of_building(5) == 1 .AND. & 956 gp_outside_of_building(1) == 0 ) THEN 957 num_gp = num_gp + 1 958 location(num_gp,1) = i * dx + 0.5_wp * dx 959 location(num_gp,2) = j * dy 960 location(num_gp,3) = k * dz - 0.5_wp * dz 961 ei(num_gp) = e(k,j,i+1) 962 dissi(num_gp) = diss(k,j,i+1) 963 de_dxi(num_gp) = 0.0_wp 964 de_dyi(num_gp) = de_dy(k,j,i+1) 965 de_dzi(num_gp) = de_dz(k,j,i+1) 966 ENDIF 967 968 ! 969 !-- If wall between gridpoint 5 and gridpoint 6, then 970 !-- then Neumann boundary condition has to be applied 971 IF ( gp_outside_of_building(5) == 1 .AND. & 972 gp_outside_of_building(6) == 0 ) THEN 973 num_gp = num_gp + 1 974 location(num_gp,1) = (i+1) * dx 975 location(num_gp,2) = j * dy + 0.5_wp * dy 976 location(num_gp,3) = k * dz - 0.5_wp * dz 977 ei(num_gp) = e(k,j,i+1) 978 dissi(num_gp) = diss(k,j,i+1) 979 de_dxi(num_gp) = de_dx(k,j,i+1) 980 de_dyi(num_gp) = 0.0_wp 981 de_dzi(num_gp) = de_dz(k,j,i+1) 982 ENDIF 983 984 IF ( gp_outside_of_building(6) == 1 .AND. & 985 gp_outside_of_building(5) == 0 ) THEN 986 num_gp = num_gp + 1 987 location(num_gp,1) = (i+1) * dx 988 location(num_gp,2) = j * dy + 0.5_wp * dy 989 location(num_gp,3) = k * dz - 0.5_wp * dz 990 ei(num_gp) = e(k,j+1,i+1) 991 dissi(num_gp) = diss(k,j+1,i+1) 992 de_dxi(num_gp) = de_dx(k,j+1,i+1) 993 de_dyi(num_gp) = 0.0_wp 994 de_dzi(num_gp) = de_dz(k,j+1,i+1) 995 ENDIF 996 997 ! 998 !-- If wall between gridpoint 2 and gridpoint 6, then 999 !-- Neumann boundary condition has to be applied 1000 IF ( gp_outside_of_building(2) == 1 .AND. & 1001 gp_outside_of_building(6) == 0 ) THEN 1002 num_gp = num_gp + 1 1003 location(num_gp,1) = i * dx + 0.5_wp * dx 1004 location(num_gp,2) = (j+1) * dy 1005 location(num_gp,3) = k * dz - 0.5_wp * dz 1006 ei(num_gp) = e(k,j+1,i) 1007 dissi(num_gp) = diss(k,j+1,i) 1008 de_dxi(num_gp) = 0.0_wp 1009 de_dyi(num_gp) = de_dy(k,j+1,i) 1010 de_dzi(num_gp) = de_dz(k,j+1,i) 1011 ENDIF 1012 1013 IF ( gp_outside_of_building(6) == 1 .AND. & 1014 gp_outside_of_building(2) == 0 ) THEN 1015 num_gp = num_gp + 1 1016 location(num_gp,1) = i * dx + 0.5_wp * dx 1017 location(num_gp,2) = (j+1) * dy 1018 location(num_gp,3) = k * dz - 0.5_wp * dz 1019 ei(num_gp) = e(k,j+1,i+1) 1020 dissi(num_gp) = diss(k,j+1,i+1) 1021 de_dxi(num_gp) = 0.0_wp 1022 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1023 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1024 ENDIF 1025 1026 ! 1027 !-- If wall between gridpoint 1 and gridpoint 2, then 1028 !-- Neumann boundary condition has to be applied 1029 IF ( gp_outside_of_building(1) == 1 .AND. & 1030 gp_outside_of_building(2) == 0 ) THEN 1031 num_gp = num_gp + 1 1032 location(num_gp,1) = i * dx 1033 location(num_gp,2) = j * dy + 0.5_wp * dy 1034 location(num_gp,3) = k * dz - 0.5_wp * dz 1035 ei(num_gp) = e(k,j,i) 1036 dissi(num_gp) = diss(k,j,i) 1037 de_dxi(num_gp) = de_dx(k,j,i) 1038 de_dyi(num_gp) = 0.0_wp 1039 de_dzi(num_gp) = de_dz(k,j,i) 1040 ENDIF 1041 1042 IF ( gp_outside_of_building(2) == 1 .AND. & 1043 gp_outside_of_building(1) == 0 ) THEN 1044 num_gp = num_gp + 1 1045 location(num_gp,1) = i * dx 1046 location(num_gp,2) = j * dy + 0.5_wp * dy 1047 location(num_gp,3) = k * dz - 0.5_wp * dz 1048 ei(num_gp) = e(k,j+1,i) 1049 dissi(num_gp) = diss(k,j+1,i) 1050 de_dxi(num_gp) = de_dx(k,j+1,i) 1051 de_dyi(num_gp) = 0.0_wp 1052 de_dzi(num_gp) = de_dz(k,j+1,i) 1053 ENDIF 1054 1055 ! 1056 !-- If wall between gridpoint 3 and gridpoint 7, then 1057 !-- Neumann boundary condition has to be applied 1058 IF ( gp_outside_of_building(3) == 1 .AND. & 1059 gp_outside_of_building(7) == 0 ) THEN 1060 num_gp = num_gp + 1 1061 location(num_gp,1) = i * dx + 0.5_wp * dx 1062 location(num_gp,2) = j * dy 1063 location(num_gp,3) = k * dz + 0.5_wp * dz 1064 ei(num_gp) = e(k+1,j,i) 1065 dissi(num_gp) = diss(k+1,j,i) 1066 de_dxi(num_gp) = 0.0_wp 1067 de_dyi(num_gp) = de_dy(k+1,j,i) 1068 de_dzi(num_gp) = de_dz(k+1,j,i) 1069 ENDIF 1070 1071 IF ( gp_outside_of_building(7) == 1 .AND. & 1072 gp_outside_of_building(3) == 0 ) THEN 1073 num_gp = num_gp + 1 1074 location(num_gp,1) = i * dx + 0.5_wp * dx 1075 location(num_gp,2) = j * dy 1076 location(num_gp,3) = k * dz + 0.5_wp * dz 1077 ei(num_gp) = e(k+1,j,i+1) 1078 dissi(num_gp) = diss(k+1,j,i+1) 1079 de_dxi(num_gp) = 0.0_wp 1080 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1081 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1082 ENDIF 1083 1084 ! 1085 !-- If wall between gridpoint 7 and gridpoint 8, then 1086 !-- Neumann boundary condition has to be applied 1087 IF ( gp_outside_of_building(7) == 1 .AND. & 1088 gp_outside_of_building(8) == 0 ) THEN 1089 num_gp = num_gp + 1 1090 location(num_gp,1) = (i+1) * dx 1091 location(num_gp,2) = j * dy + 0.5_wp * dy 1092 location(num_gp,3) = k * dz + 0.5_wp * dz 1093 ei(num_gp) = e(k+1,j,i+1) 1094 dissi(num_gp) = diss(k+1,j,i+1) 1095 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1096 de_dyi(num_gp) = 0.0_wp 1097 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1098 ENDIF 1099 1100 IF ( gp_outside_of_building(8) == 1 .AND. & 1101 gp_outside_of_building(7) == 0 ) THEN 1102 num_gp = num_gp + 1 1103 location(num_gp,1) = (i+1) * dx 1104 location(num_gp,2) = j * dy + 0.5_wp * dy 1105 location(num_gp,3) = k * dz + 0.5_wp * dz 1106 ei(num_gp) = e(k+1,j+1,i+1) 1107 dissi(num_gp) = diss(k+1,j+1,i+1) 1108 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1109 de_dyi(num_gp) = 0.0_wp 1110 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1111 ENDIF 1112 1113 ! 1114 !-- If wall between gridpoint 4 and gridpoint 8, then 1115 !-- Neumann boundary condition has to be applied 1116 IF ( gp_outside_of_building(4) == 1 .AND. & 1117 gp_outside_of_building(8) == 0 ) THEN 1118 num_gp = num_gp + 1 1119 location(num_gp,1) = i * dx + 0.5_wp * dx 1120 location(num_gp,2) = (j+1) * dy 1121 location(num_gp,3) = k * dz + 0.5_wp * dz 1122 ei(num_gp) = e(k+1,j+1,i) 1123 dissi(num_gp) = diss(k+1,j+1,i) 1124 de_dxi(num_gp) = 0.0_wp 1125 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1126 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1127 ENDIF 1128 1129 IF ( gp_outside_of_building(8) == 1 .AND. & 1130 gp_outside_of_building(4) == 0 ) THEN 1131 num_gp = num_gp + 1 1132 location(num_gp,1) = i * dx + 0.5_wp * dx 1133 location(num_gp,2) = (j+1) * dy 1134 location(num_gp,3) = k * dz + 0.5_wp * dz 1135 ei(num_gp) = e(k+1,j+1,i+1) 1136 dissi(num_gp) = diss(k+1,j+1,i+1) 1137 de_dxi(num_gp) = 0.0_wp 1138 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1139 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1140 ENDIF 1141 1142 ! 1143 !-- If wall between gridpoint 3 and gridpoint 4, then 1144 !-- Neumann boundary condition has to be applied 1145 IF ( gp_outside_of_building(3) == 1 .AND. & 1146 gp_outside_of_building(4) == 0 ) THEN 1147 num_gp = num_gp + 1 1148 location(num_gp,1) = i * dx 1149 location(num_gp,2) = j * dy + 0.5_wp * dy 1150 location(num_gp,3) = k * dz + 0.5_wp * dz 1151 ei(num_gp) = e(k+1,j,i) 1152 dissi(num_gp) = diss(k+1,j,i) 1153 de_dxi(num_gp) = de_dx(k+1,j,i) 1154 de_dyi(num_gp) = 0.0_wp 1155 de_dzi(num_gp) = de_dz(k+1,j,i) 1156 ENDIF 1157 1158 IF ( gp_outside_of_building(4) == 1 .AND. & 1159 gp_outside_of_building(3) == 0 ) THEN 1160 num_gp = num_gp + 1 1161 location(num_gp,1) = i * dx 1162 location(num_gp,2) = j * dy + 0.5_wp * dy 1163 location(num_gp,3) = k * dz + 0.5_wp * dz 1164 ei(num_gp) = e(k+1,j+1,i) 1165 dissi(num_gp) = diss(k+1,j+1,i) 1166 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1167 de_dyi(num_gp) = 0.0_wp 1168 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1169 ENDIF 1170 1171 ! 1172 !-- If wall between gridpoint 1 and gridpoint 3, then 1173 !-- Neumann boundary condition has to be applied 1174 !-- (only one case as only building beneath is possible) 1175 IF ( gp_outside_of_building(1) == 0 .AND. & 1176 gp_outside_of_building(3) == 1 ) THEN 1177 num_gp = num_gp + 1 1178 location(num_gp,1) = i * dx 1179 location(num_gp,2) = j * dy 1180 location(num_gp,3) = k * dz 1181 ei(num_gp) = e(k+1,j,i) 1182 dissi(num_gp) = diss(k+1,j,i) 1183 de_dxi(num_gp) = de_dx(k+1,j,i) 1184 de_dyi(num_gp) = de_dy(k+1,j,i) 1185 de_dzi(num_gp) = 0.0_wp 1186 ENDIF 1187 1188 ! 1189 !-- If wall between gridpoint 5 and gridpoint 7, then 1190 !-- Neumann boundary condition has to be applied 1191 !-- (only one case as only building beneath is possible) 1192 IF ( gp_outside_of_building(5) == 0 .AND. & 1193 gp_outside_of_building(7) == 1 ) THEN 1194 num_gp = num_gp + 1 1195 location(num_gp,1) = (i+1) * dx 1196 location(num_gp,2) = j * dy 1197 location(num_gp,3) = k * dz 1198 ei(num_gp) = e(k+1,j,i+1) 1199 dissi(num_gp) = diss(k+1,j,i+1) 1200 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1201 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1202 de_dzi(num_gp) = 0.0_wp 1203 ENDIF 1204 1205 ! 1206 !-- If wall between gridpoint 2 and gridpoint 4, then 1207 !-- Neumann boundary condition has to be applied 1208 !-- (only one case as only building beneath is possible) 1209 IF ( gp_outside_of_building(2) == 0 .AND. & 1210 gp_outside_of_building(4) == 1 ) THEN 1211 num_gp = num_gp + 1 1212 location(num_gp,1) = i * dx 1213 location(num_gp,2) = (j+1) * dy 1214 location(num_gp,3) = k * dz 1215 ei(num_gp) = e(k+1,j+1,i) 1216 dissi(num_gp) = diss(k+1,j+1,i) 1217 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1218 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1219 de_dzi(num_gp) = 0.0_wp 1220 ENDIF 1221 1222 ! 1223 !-- If wall between gridpoint 6 and gridpoint 8, then 1224 !-- Neumann boundary condition has to be applied 1225 !-- (only one case as only building beneath is possible) 1226 IF ( gp_outside_of_building(6) == 0 .AND. & 1227 gp_outside_of_building(8) == 1 ) THEN 1228 num_gp = num_gp + 1 1229 location(num_gp,1) = (i+1) * dx 1230 location(num_gp,2) = (j+1) * dy 1231 location(num_gp,3) = k * dz 1232 ei(num_gp) = e(k+1,j+1,i+1) 1233 dissi(num_gp) = diss(k+1,j+1,i+1) 1234 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1235 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1236 de_dzi(num_gp) = 0.0_wp 1237 ENDIF 1238 1239 ! 1240 !-- Carry out the interpolation 1241 IF ( num_gp == 1 ) THEN 1242 ! 1243 !-- If only one of the gridpoints is situated outside of the 1244 !-- building, it follows that the values at the particle 1245 !-- location are the same as the gridpoint values 1246 e_int(n) = ei(num_gp) 1247 diss_int(n) = dissi(num_gp) 1248 de_dx_int(n) = de_dxi(num_gp) 1249 de_dy_int(n) = de_dyi(num_gp) 1250 de_dz_int(n) = de_dzi(num_gp) 1251 ! 1252 !-- Set flag for stochastic equation which disables calculation 1253 !-- of drift and memory term near topography. 1254 term_1_2(n) = 0.0_wp 1255 ELSE IF ( num_gp > 1 ) THEN 1256 1257 d_sum = 0.0_wp 1258 ! 1259 !-- Evaluation of the distances between the gridpoints 1260 !-- contributing to the interpolated values, and the particle 1261 !-- location 1262 DO agp = 1, num_gp 1263 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1264 + ( particles(n)%y-location(agp,2) )**2 & 1265 + ( zv(n)-location(agp,3) )**2 1266 d_sum = d_sum + d_gp_pl(agp) 1267 ENDDO 1268 1269 ! 1270 !-- Finally the interpolation can be carried out 1271 e_int(n) = 0.0_wp 1272 diss_int(n) = 0.0_wp 1273 de_dx_int(n) = 0.0_wp 1274 de_dy_int(n) = 0.0_wp 1275 de_dz_int(n) = 0.0_wp 1276 DO agp = 1, num_gp 1277 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1278 ei(agp) / ( (num_gp-1) * d_sum ) 1279 diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1280 dissi(agp) / ( (num_gp-1) * d_sum ) 1281 de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1282 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1283 de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1284 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1285 de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1286 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1287 ENDDO 1288 1289 ENDIF 1290 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1291 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1292 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1293 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1294 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1242 1295 ! 1243 1296 !-- Set flag for stochastic equation which disables calculation 1244 1297 !-- of drift and memory term near topography. 1245 1298 term_1_2(n) = 0.0_wp 1246 ELSE IF ( num_gp > 1 ) THEN 1247 1248 d_sum = 0.0_wp 1249 ! 1250 !-- Evaluation of the distances between the gridpoints 1251 !-- contributing to the interpolated values, and the particle 1252 !-- location 1253 DO agp = 1, num_gp 1254 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1255 + ( particles(n)%y-location(agp,2) )**2 & 1256 + ( zv(n)-location(agp,3) )**2 1257 d_sum = d_sum + d_gp_pl(agp) 1258 ENDDO 1259 1260 ! 1261 !-- Finally the interpolation can be carried out 1262 e_int(n) = 0.0_wp 1263 diss_int(n) = 0.0_wp 1264 de_dx_int(n) = 0.0_wp 1265 de_dy_int(n) = 0.0_wp 1266 de_dz_int(n) = 0.0_wp 1267 DO agp = 1, num_gp 1268 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1269 ei(agp) / ( (num_gp-1) * d_sum ) 1270 diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1271 dissi(agp) / ( (num_gp-1) * d_sum ) 1272 de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1273 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1274 de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1275 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1276 de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1277 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1278 ENDDO 1279 1280 ENDIF 1281 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1282 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1283 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1284 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1285 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1286 ! 1287 !-- Set flag for stochastic equation which disables calculation 1288 !-- of drift and memory term near topography. 1289 term_1_2(n) = 0.0_wp 1290 ENDIF 1299 ENDIF 1300 ENDDO 1291 1301 ENDDO 1292 1302 ENDIF … … 1345 1355 ENDDO 1346 1356 1347 DO n = 1, number_of_particles1348 1349 rg(n,1) = random_gauss( iran_part, 5.0_wp )1350 rg(n,2) = random_gauss( iran_part, 5.0_wp )1351 rg(n,3) = random_gauss( iran_part, 5.0_wp )1352 1357 DO nb = 0, 7 1358 DO n = start_index(nb), end_index(nb) 1359 rg(n,1) = random_gauss( iran_part, 5.0_wp ) 1360 rg(n,2) = random_gauss( iran_part, 5.0_wp ) 1361 rg(n,3) = random_gauss( iran_part, 5.0_wp ) 1362 ENDDO 1353 1363 ENDDO 1354 1364 1355 DO n = 1, number_of_particles 1356 ! 1357 !-- Calculate the Lagrangian timescale according to Weil et al. (2004). 1358 lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / & 1359 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp ) 1360 1361 ! 1362 !-- Calculate the next particle timestep. dt_gap is the time needed to 1363 !-- complete the current LES timestep. 1364 dt_gap = dt_3d - particles(n)%dt_sum 1365 dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap ) 1366 1367 ! 1368 !-- The particle timestep should not be too small in order to prevent 1369 !-- the number of particle timesteps of getting too large 1370 IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap ) THEN 1371 dt_particle(n) = dt_min_part 1372 ENDIF 1373 1374 ! 1375 !-- Calculate the SGS velocity components 1376 IF ( particles(n)%age == 0.0_wp ) THEN 1377 ! 1378 !-- For new particles the SGS components are derived from the SGS 1379 !-- TKE. Limit the Gaussian random number to the interval 1380 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 1381 !-- from becoming unrealistically large. 1382 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1383 ( rg(n,1) - 1.0_wp ) 1384 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1385 ( rg(n,2) - 1.0_wp ) 1386 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1387 ( rg(n,3) - 1.0_wp ) 1388 1389 ELSE 1390 ! 1391 !-- Restriction of the size of the new timestep: compared to the 1392 !-- previous timestep the increase must not exceed 200% 1393 1394 dt_particle_m = particles(n)%age - particles(n)%age_m 1395 IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN 1396 dt_particle(n) = 2.0_wp * dt_particle_m 1365 DO nb = 0, 7 1366 DO n = start_index(nb), end_index(nb) 1367 1368 ! 1369 !-- Calculate the Lagrangian timescale according to Weil et al. (2004). 1370 lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / & 1371 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp ) 1372 1373 ! 1374 !-- Calculate the next particle timestep. dt_gap is the time needed to 1375 !-- complete the current LES timestep. 1376 dt_gap = dt_3d - particles(n)%dt_sum 1377 dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap ) 1378 1379 ! 1380 !-- The particle timestep should not be too small in order to prevent 1381 !-- the number of particle timesteps of getting too large 1382 IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap ) THEN 1383 dt_particle(n) = dt_min_part 1397 1384 ENDIF 1398 1385 1399 1386 ! 1400 !-- For old particles the SGS components are correlated with the 1401 !-- values from the previous timestep. Random numbers have also to 1402 !-- be limited (see above). 1403 !-- As negative values for the subgrid TKE are not allowed, the 1404 !-- change of the subgrid TKE with time cannot be smaller than 1405 !-- -e_int(n)/dt_particle. This value is used as a lower boundary 1406 !-- value for the change of TKE 1407 1408 de_dt_min = - e_int(n) / dt_particle(n) 1409 1410 de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m 1411 1412 IF ( de_dt < de_dt_min ) THEN 1413 de_dt = de_dt_min 1387 !-- Calculate the SGS velocity components 1388 IF ( particles(n)%age == 0.0_wp ) THEN 1389 ! 1390 !-- For new particles the SGS components are derived from the SGS 1391 !-- TKE. Limit the Gaussian random number to the interval 1392 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 1393 !-- from becoming unrealistically large. 1394 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 1395 + 1E-20_wp ) * & 1396 ( rg(n,1) - 1.0_wp ) 1397 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 1398 + 1E-20_wp ) * & 1399 ( rg(n,2) - 1.0_wp ) 1400 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 1401 + 1E-20_wp ) * & 1402 ( rg(n,3) - 1.0_wp ) 1403 1404 ELSE 1405 ! 1406 !-- Restriction of the size of the new timestep: compared to the 1407 !-- previous timestep the increase must not exceed 200%. First, 1408 !-- check if age > age_m, in order to prevent that particles get zero 1409 !-- timestep. 1410 dt_particle_m = MERGE( dt_particle(n), & 1411 particles(n)%age - particles(n)%age_m, & 1412 particles(n)%age - particles(n)%age_m < & 1413 1E-8_wp ) 1414 IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN 1415 dt_particle(n) = 2.0_wp * dt_particle_m 1416 ENDIF 1417 1418 ! 1419 !-- For old particles the SGS components are correlated with the 1420 !-- values from the previous timestep. Random numbers have also to 1421 !-- be limited (see above). 1422 !-- As negative values for the subgrid TKE are not allowed, the 1423 !-- change of the subgrid TKE with time cannot be smaller than 1424 !-- -e_int(n)/dt_particle. This value is used as a lower boundary 1425 !-- value for the change of TKE 1426 de_dt_min = - e_int(n) / dt_particle(n) 1427 1428 de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m 1429 1430 IF ( de_dt < de_dt_min ) THEN 1431 de_dt = de_dt_min 1432 ENDIF 1433 1434 CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),& 1435 de_dx_int(n), de_dt, diss_int(n), & 1436 dt_particle(n), rg(n,1), term_1_2(n) ) 1437 1438 CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),& 1439 de_dy_int(n), de_dt, diss_int(n), & 1440 dt_particle(n), rg(n,2), term_1_2(n) ) 1441 1442 CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),& 1443 de_dz_int(n), de_dt, diss_int(n), & 1444 dt_particle(n), rg(n,3), term_1_2(n) ) 1445 1414 1446 ENDIF 1415 1447 1416 CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n), & 1417 de_dx_int(n), de_dt, diss_int(n), & 1418 dt_particle(n), rg(n,1), term_1_2(n) ) 1419 1420 CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n), & 1421 de_dy_int(n), de_dt, diss_int(n), & 1422 dt_particle(n), rg(n,2), term_1_2(n) ) 1423 1424 CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n), & 1425 de_dz_int(n), de_dt, diss_int(n), & 1426 dt_particle(n), rg(n,3), term_1_2(n) ) 1427 1428 ENDIF 1429 1430 u_int(n) = u_int(n) + particles(n)%rvar1 1431 v_int(n) = v_int(n) + particles(n)%rvar2 1432 w_int(n) = w_int(n) + particles(n)%rvar3 1433 ! 1434 !-- Store the SGS TKE of the current timelevel which is needed for 1435 !-- for calculating the SGS particle velocities at the next timestep 1436 particles(n)%e_m = e_int(n) 1448 u_int(n) = u_int(n) + particles(n)%rvar1 1449 v_int(n) = v_int(n) + particles(n)%rvar2 1450 w_int(n) = w_int(n) + particles(n)%rvar3 1451 ! 1452 !-- Store the SGS TKE of the current timelevel which is needed for 1453 !-- for calculating the SGS particle velocities at the next timestep 1454 particles(n)%e_m = e_int(n) 1455 ENDDO 1437 1456 ENDDO 1438 1457 … … 1444 1463 1445 1464 ENDIF 1446 !1447 !-- Store the old age of the particle ( needed to prevent that a1448 !-- particle crosses several PEs during one timestep, and for the1449 !-- evaluation of the subgrid particle velocity fluctuations )1450 particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age1451 1465 1452 1466 dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio 1453 1467 1454 1468 IF ( ANY( dens_ratio == 0.0_wp ) ) THEN 1455 DO n = 1, number_of_particles 1456 1457 ! 1458 !-- Particle advection 1459 IF ( dens_ratio(n) == 0.0_wp ) THEN 1460 ! 1461 !-- Pure passive transport (without particle inertia) 1462 particles(n)%x = xv(n) + u_int(n) * dt_particle(n) 1463 particles(n)%y = yv(n) + v_int(n) * dt_particle(n) 1464 particles(n)%z = zv(n) + w_int(n) * dt_particle(n) 1465 1466 particles(n)%speed_x = u_int(n) 1467 particles(n)%speed_y = v_int(n) 1468 particles(n)%speed_z = w_int(n) 1469 1470 ELSE 1469 DO nb = 0, 7 1470 DO n = start_index(nb), end_index(nb) 1471 1472 ! 1473 !-- Particle advection 1474 IF ( dens_ratio(n) == 0.0_wp ) THEN 1475 ! 1476 !-- Pure passive transport (without particle inertia) 1477 particles(n)%x = xv(n) + u_int(n) * dt_particle(n) 1478 particles(n)%y = yv(n) + v_int(n) * dt_particle(n) 1479 particles(n)%z = zv(n) + w_int(n) * dt_particle(n) 1480 1481 particles(n)%speed_x = u_int(n) 1482 particles(n)%speed_y = v_int(n) 1483 particles(n)%speed_z = w_int(n) 1484 1485 ELSE 1486 ! 1487 !-- Transport of particles with inertia 1488 particles(n)%x = particles(n)%x + particles(n)%speed_x * & 1489 dt_particle(n) 1490 particles(n)%y = particles(n)%y + particles(n)%speed_y * & 1491 dt_particle(n) 1492 particles(n)%z = particles(n)%z + particles(n)%speed_z * & 1493 dt_particle(n) 1494 1495 ! 1496 !-- Update of the particle velocity 1497 IF ( cloud_droplets ) THEN 1498 ! 1499 !-- Terminal velocity is computed for vertical direction (Rogers et 1500 !-- al., 1993, J. Appl. Meteorol.) 1501 diameter = particles(n)%radius * 2000.0_wp !diameter in mm 1502 IF ( diameter <= d0_rog ) THEN 1503 w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) ) 1504 ELSE 1505 w_s = a_rog - b_rog * EXP( -c_rog * diameter ) 1506 ENDIF 1507 1508 ! 1509 !-- If selected, add random velocities following Soelch and Kaercher 1510 !-- (2010, Q. J. R. Meteorol. Soc.) 1511 IF ( use_sgs_for_particles ) THEN 1512 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 1513 RL = EXP( -1.0_wp * dt_3d / lagr_timescale ) 1514 sigma = SQRT( e(kp,jp,ip) ) 1515 1516 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1517 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1518 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1519 1520 particles(n)%rvar1 = RL * particles(n)%rvar1 + & 1521 SQRT( 1.0_wp - RL**2 ) * sigma * rg1 1522 particles(n)%rvar2 = RL * particles(n)%rvar2 + & 1523 SQRT( 1.0_wp - RL**2 ) * sigma * rg2 1524 particles(n)%rvar3 = RL * particles(n)%rvar3 + & 1525 SQRT( 1.0_wp - RL**2 ) * sigma * rg3 1526 1527 particles(n)%speed_x = u_int(n) + particles(n)%rvar1 1528 particles(n)%speed_y = v_int(n) + particles(n)%rvar2 1529 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s 1530 ELSE 1531 particles(n)%speed_x = u_int(n) 1532 particles(n)%speed_y = v_int(n) 1533 particles(n)%speed_z = w_int(n) - w_s 1534 ENDIF 1535 1536 ELSE 1537 1538 IF ( use_sgs_for_particles ) THEN 1539 exp_arg = particle_groups(particles(n)%group)%exp_arg 1540 exp_term = EXP( -exp_arg * dt_particle(n) ) 1541 ELSE 1542 exp_arg = particle_groups(particles(n)%group)%exp_arg 1543 exp_term = particle_groups(particles(n)%group)%exp_term 1544 ENDIF 1545 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 1546 u_int(n) * ( 1.0_wp - exp_term ) 1547 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 1548 v_int(n) * ( 1.0_wp - exp_term ) 1549 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 1550 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * & 1551 g / exp_arg ) * ( 1.0_wp - exp_term ) 1552 ENDIF 1553 1554 ENDIF 1555 ENDDO 1556 ENDDO 1557 1558 ELSE 1559 1560 DO nb = 0, 7 1561 DO n = start_index(nb), end_index(nb) 1471 1562 ! 1472 1563 !-- Transport of particles with inertia 1473 particles(n)%x = particles(n)%x + particles(n)%speed_x * & 1474 dt_particle(n) 1475 particles(n)%y = particles(n)%y + particles(n)%speed_y * & 1476 dt_particle(n) 1477 particles(n)%z = particles(n)%z + particles(n)%speed_z * & 1478 dt_particle(n) 1479 1564 particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n) 1565 particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n) 1566 particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n) 1480 1567 ! 1481 1568 !-- Update of the particle velocity 1482 1569 IF ( cloud_droplets ) THEN 1483 1570 ! 1484 !-- Terminal velocity is computed for vertical direction (Rogers et 1485 !-- al.,1993, J. Appl. Meteorol.)1571 !-- Terminal velocity is computed for vertical direction (Rogers et al., 1572 !-- 1993, J. Appl. Meteorol.) 1486 1573 diameter = particles(n)%radius * 2000.0_wp !diameter in mm 1487 1574 IF ( diameter <= d0_rog ) THEN … … 1495 1582 !-- (2010, Q. J. R. Meteorol. Soc.) 1496 1583 IF ( use_sgs_for_particles ) THEN 1497 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )1498 RL = EXP( -1.0_wp * dt_3d / lagr_timescale )1499 sigma = SQRT( e(kp,jp,ip) )1500 1501 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1502 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1503 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1504 1505 particles(n)%rvar1 = RL * particles(n)%rvar1 +&1506 SQRT( 1.0_wp - RL**2 ) * sigma * rg11507 particles(n)%rvar2 = RL * particles(n)%rvar2 +&1508 SQRT( 1.0_wp - RL**2 ) * sigma * rg21509 particles(n)%rvar3 = RL * particles(n)%rvar3 +&1510 SQRT( 1.0_wp - RL**2 ) * sigma * rg31511 1512 particles(n)%speed_x = u_int(n) + particles(n)%rvar11513 particles(n)%speed_y = v_int(n) + particles(n)%rvar21514 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s1584 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 1585 RL = EXP( -1.0_wp * dt_3d / lagr_timescale ) 1586 sigma = SQRT( e(kp,jp,ip) ) 1587 1588 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1589 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1590 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1591 1592 particles(n)%rvar1 = RL * particles(n)%rvar1 + & 1593 SQRT( 1.0_wp - RL**2 ) * sigma * rg1 1594 particles(n)%rvar2 = RL * particles(n)%rvar2 + & 1595 SQRT( 1.0_wp - RL**2 ) * sigma * rg2 1596 particles(n)%rvar3 = RL * particles(n)%rvar3 + & 1597 SQRT( 1.0_wp - RL**2 ) * sigma * rg3 1598 1599 particles(n)%speed_x = u_int(n) + particles(n)%rvar1 1600 particles(n)%speed_y = v_int(n) + particles(n)%rvar2 1601 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s 1515 1602 ELSE 1516 particles(n)%speed_x = u_int(n)1517 particles(n)%speed_y = v_int(n)1518 particles(n)%speed_z = w_int(n) - w_s1603 particles(n)%speed_x = u_int(n) 1604 particles(n)%speed_y = v_int(n) 1605 particles(n)%speed_z = w_int(n) - w_s 1519 1606 ENDIF 1520 1607 … … 1528 1615 exp_term = particle_groups(particles(n)%group)%exp_term 1529 1616 ENDIF 1530 particles(n)%speed_x = particles(n)%speed_x * exp_term + &1617 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 1531 1618 u_int(n) * ( 1.0_wp - exp_term ) 1532 particles(n)%speed_y = particles(n)%speed_y * exp_term + &1619 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 1533 1620 v_int(n) * ( 1.0_wp - exp_term ) 1534 particles(n)%speed_z = particles(n)%speed_z * exp_term + &1535 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &1536 g /exp_arg ) * ( 1.0_wp - exp_term )1621 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 1622 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / & 1623 exp_arg ) * ( 1.0_wp - exp_term ) 1537 1624 ENDIF 1538 1625 ENDDO 1626 ENDDO 1627 1628 ENDIF 1629 1630 ! 1631 !-- Store the old age of the particle ( needed to prevent that a 1632 !-- particle crosses several PEs during one timestep, and for the 1633 !-- evaluation of the subgrid particle velocity fluctuations ) 1634 particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age 1635 1636 DO nb = 0, 7 1637 DO n = start_index(nb), end_index(nb) 1638 ! 1639 !-- Increment the particle age and the total time that the particle 1640 !-- has advanced within the particle timestep procedure 1641 particles(n)%age = particles(n)%age + dt_particle(n) 1642 particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n) 1643 1644 ! 1645 !-- Check whether there is still a particle that has not yet completed 1646 !-- the total LES timestep 1647 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN 1648 dt_3d_reached_l = .FALSE. 1539 1649 ENDIF 1540 1650 1541 1651 ENDDO 1542 1543 ELSE1544 1545 DO n = 1, number_of_particles1546 1547 !-- Transport of particles with inertia1548 particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)1549 particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)1550 particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)1551 !1552 !-- Update of the particle velocity1553 IF ( cloud_droplets ) THEN1554 !1555 !-- Terminal velocity is computed for vertical direction (Rogers et al.,1556 !-- 1993, J. Appl. Meteorol.)1557 diameter = particles(n)%radius * 2000.0_wp !diameter in mm1558 IF ( diameter <= d0_rog ) THEN1559 w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )1560 ELSE1561 w_s = a_rog - b_rog * EXP( -c_rog * diameter )1562 ENDIF1563 1564 !1565 !-- If selected, add random velocities following Soelch and Kaercher1566 !-- (2010, Q. J. R. Meteorol. Soc.)1567 IF ( use_sgs_for_particles ) THEN1568 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )1569 RL = EXP( -1.0_wp * dt_3d / lagr_timescale )1570 sigma = SQRT( e(kp,jp,ip) )1571 1572 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1573 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1574 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp1575 1576 particles(n)%rvar1 = RL * particles(n)%rvar1 + &1577 SQRT( 1.0_wp - RL**2 ) * sigma * rg11578 particles(n)%rvar2 = RL * particles(n)%rvar2 + &1579 SQRT( 1.0_wp - RL**2 ) * sigma * rg21580 particles(n)%rvar3 = RL * particles(n)%rvar3 + &1581 SQRT( 1.0_wp - RL**2 ) * sigma * rg31582 1583 particles(n)%speed_x = u_int(n) + particles(n)%rvar11584 particles(n)%speed_y = v_int(n) + particles(n)%rvar21585 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s1586 ELSE1587 particles(n)%speed_x = u_int(n)1588 particles(n)%speed_y = v_int(n)1589 particles(n)%speed_z = w_int(n) - w_s1590 ENDIF1591 1592 ELSE1593 1594 IF ( use_sgs_for_particles ) THEN1595 exp_arg = particle_groups(particles(n)%group)%exp_arg1596 exp_term = EXP( -exp_arg * dt_particle(n) )1597 ELSE1598 exp_arg = particle_groups(particles(n)%group)%exp_arg1599 exp_term = particle_groups(particles(n)%group)%exp_term1600 ENDIF1601 particles(n)%speed_x = particles(n)%speed_x * exp_term + &1602 u_int(n) * ( 1.0_wp - exp_term )1603 particles(n)%speed_y = particles(n)%speed_y * exp_term + &1604 v_int(n) * ( 1.0_wp - exp_term )1605 particles(n)%speed_z = particles(n)%speed_z * exp_term + &1606 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &1607 exp_arg ) * ( 1.0_wp - exp_term )1608 ENDIF1609 1610 ENDDO1611 1612 ENDIF1613 1614 DO n = 1, number_of_particles1615 !1616 !-- Increment the particle age and the total time that the particle1617 !-- has advanced within the particle timestep procedure1618 particles(n)%age = particles(n)%age + dt_particle(n)1619 particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)1620 1621 !1622 !-- Check whether there is still a particle that has not yet completed1623 !-- the total LES timestep1624 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN1625 dt_3d_reached_l = .FALSE.1626 ENDIF1627 1628 1652 ENDDO 1629 1653 … … 1693 1717 !-- subgrid-scale velocity component is set to zero, in order to prevent a 1694 1718 !-- velocity build-up. 1695 1696 1719 !-- This case, set also previous subgrid-scale component to zero. 1697 1720 v_sgs = v_sgs * fac + term1 + term2 + term3 -
palm/trunk/SOURCE/lpm_pack_arrays.f90
r2101 r2417 25 25 ! ----------------- 26 26 ! $Id$ 27 ! New routine which sorts particles into particles that completed and not 28 ! completed the LES timestep. 29 ! 30 ! 2101 2017-01-05 16:42:31Z suehring 27 31 ! 28 32 ! 2000 2016-08-20 18:09:15Z knoop … … 61 65 !> deletion and move data with higher index values to these free indices. 62 66 !> Determine the new number of particles. 67 !> Moreover, particles are also sorted into groups finished and not finished 68 !> its timestep. 63 69 !------------------------------------------------------------------------------! 64 70 MODULE lpm_pack_arrays_mod … … 70 76 71 77 PRIVATE 72 PUBLIC lpm_pack_all_arrays, lpm_pack_arrays 78 PUBLIC lpm_pack_all_arrays, lpm_pack_arrays, lpm_sort 73 79 74 80 INTERFACE lpm_pack_all_arrays … … 80 86 END INTERFACE lpm_pack_arrays 81 87 82 CONTAINS 88 INTERFACE lpm_sort 89 MODULE PROCEDURE lpm_sort 90 END INTERFACE lpm_sort 91 92 93 CONTAINS 83 94 84 95 !------------------------------------------------------------------------------! … … 236 247 END SUBROUTINE lpm_pack_and_sort 237 248 249 !------------------------------------------------------------------------------! 250 ! Description: 251 ! ------------ 252 !> Sort particles in each sub-grid box into two groups: particles that already 253 !> completed the LES timestep, and particles that need further timestepping to 254 !> complete the LES timestep. 255 !------------------------------------------------------------------------------! 256 SUBROUTINE lpm_sort 257 258 USE control_parameters, & 259 ONLY: dt_3d 260 261 USE indices, & 262 ONLY: nxl, nxr, nys, nyn, nzb, nzt 263 264 USE kinds 265 266 IMPLICIT NONE 267 268 INTEGER(iwp) :: end_index !< particle end index for each sub-box 269 INTEGER(iwp) :: i !< index of particle grid box in x-direction 270 INTEGER(iwp) :: j !< index of particle grid box in y-direction 271 INTEGER(iwp) :: k !< index of particle grid box in z-direction 272 INTEGER(iwp) :: n !< running index for number of particles 273 INTEGER(iwp) :: nb !< index of subgrid boux 274 INTEGER(iwp) :: nf !< indices for particles in each sub-box that already finalized their substeps 275 INTEGER(iwp) :: nnf !< indices for particles in each sub-box that need further treatment 276 INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps 277 INTEGER(iwp) :: start_index !< particle start index for each sub-box 278 279 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles !< temporary particle array 280 281 DO i = nxl, nxr 282 DO j = nys, nyn 283 DO k = nzb+1, nzt 284 285 number_of_particles = prt_count(k,j,i) 286 IF ( number_of_particles <= 0 ) CYCLE 287 288 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 289 290 DO nb = 0, 7 291 ! 292 !-- Obtain start and end index for each subgrid box 293 start_index = grid_particles(k,j,i)%start_index(nb) 294 end_index = grid_particles(k,j,i)%end_index(nb) 295 ! 296 !-- Allocate temporary array used for sorting. 297 ALLOCATE( sort_particles(start_index:end_index) ) 298 ! 299 !-- Determine number of particles already completed the LES 300 !-- timestep, and write them into a temporary array. 301 nf = start_index 302 num_finalized = 0 303 DO n = start_index, end_index 304 IF ( dt_3d - particles(n)%dt_sum < 1E-8_wp ) THEN 305 sort_particles(nf) = particles(n) 306 nf = nf + 1 307 num_finalized = num_finalized + 1 308 ENDIF 309 ENDDO 310 ! 311 !-- Determine number of particles that not completed the LES 312 !-- timestep, and write them into a temporary array. 313 nnf = nf 314 DO n = start_index, end_index 315 IF ( dt_3d - particles(n)%dt_sum > 1E-8_wp ) THEN 316 sort_particles(nnf) = particles(n) 317 nnf = nnf + 1 318 ENDIF 319 ENDDO 320 ! 321 !-- Write back sorted particles 322 particles(start_index:end_index) = & 323 sort_particles(start_index:end_index) 324 ! 325 !-- Determine updated start_index, used to masked already 326 !-- completed particles. 327 grid_particles(k,j,i)%start_index(nb) = & 328 grid_particles(k,j,i)%start_index(nb) & 329 + num_finalized 330 ! 331 !-- Deallocate dummy array 332 DEALLOCATE ( sort_particles ) 333 ! 334 !-- Finally, if number of non-completed particles is non zero 335 !-- in any of the sub-boxes, set control flag appropriately. 336 IF ( nnf > nf ) & 337 grid_particles(k,j,i)%time_loop_done = .FALSE. 338 339 ENDDO 340 ENDDO 341 ENDDO 342 ENDDO 343 344 END SUBROUTINE lpm_sort 345 238 346 239 347 END module lpm_pack_arrays_mod -
palm/trunk/SOURCE/user_lpm_advec.f90
r2101 r2417 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Particle loops adapted for sub-box structure, i.e. for each sub-box the 28 ! particle loop runs from start_index up to end_index instead from 1 to 29 ! number_of_particles. 30 ! 31 ! 2101 2017-01-05 16:42:31Z suehring 27 32 ! 28 33 ! 2000 2016-08-20 18:09:15Z knoop … … 70 75 IMPLICIT NONE 71 76 72 INTEGER(iwp) :: ip !< 73 INTEGER(iwp) :: jp !< 74 INTEGER(iwp) :: kp !< 75 INTEGER(iwp) :: n !< 77 INTEGER(iwp) :: ip !< index of particle grid box, x-direction 78 INTEGER(iwp) :: jp !< index of particle grid box, y-direction 79 INTEGER(iwp) :: kp !< index of particle grid box, z-direction 80 INTEGER(iwp) :: n !< particle index 81 INTEGER(iwp) :: nb !< index of sub-box particles are sorted in 82 83 INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current sub-box 84 INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current sub-box 76 85 77 86 ! … … 79 88 ! number_of_particles = prt_count(kp,jp,ip) 80 89 ! particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 90 ! 91 ! start_index = grid_particles(kp,jp,ip)%start_index 92 ! end_index = grid_particles(kp,jp,ip)%end_index 93 ! 81 94 ! IF ( number_of_particles <= 0 ) CYCLE 82 95 ! DO n = 1, number_of_particles 83 ! 96 ! DO nb = 0, 7 97 ! DO n = start_index(nb), end_index(nb) 98 ! particles(n)%xxx = 99 ! ENDDO 84 100 ! ENDDO 85 101
Note: See TracChangeset
for help on using the changeset viewer.
