Changeset 2698 for palm/trunk/SOURCE/lpm_advec.f90
- Timestamp:
- Dec 14, 2017 6:46:24 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_advec.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r2696 r2698 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Particle interpolations at walls in case of SGS velocities revised and not 23 ! required parts are removed. (responsible Philipp Thiele) 24 ! 25 ! Bugfix in get_topography_top_index 23 26 ! 24 27 ! Former revisions: … … 149 152 150 153 USE indices, & 151 ONLY: nzb, nzt 154 ONLY: nzb, nzt, wall_flags_0 152 155 153 156 USE kinds … … 163 166 164 167 USE surface_mod, & 165 ONLY: get_topography_top_index , surf_def_h, surf_lsm_h, surf_usm_h168 ONLY: get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h 166 169 167 170 IMPLICIT NONE 171 172 LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall 168 173 169 174 INTEGER(iwp) :: agp !< loop variable … … 302 307 !-- above topography (Prandtl-layer height) 303 308 !-- Determine vertical index of topography top 304 k_wall = get_topography_top_index ( jp, ip, 's' )309 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 305 310 306 311 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 397 402 ! 398 403 !-- Determine vertical index of topography top 399 k_wall = get_topography_top_index ( jp,ip, 's' )404 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 400 405 401 406 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN … … 528 533 !-- velocities 529 534 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 530 531 IF ( topography == 'flat' ) THEN 532 533 DO nb = 0,7 534 535 536 DO nb = 0,7 537 538 subbox_at_wall = .FALSE. 539 ! 540 !-- In case of topography check if subbox is adjacent to a wall 541 IF ( .NOT. topography == 'flat' ) THEN 542 i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) ) 543 j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) ) 544 k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) ) 545 IF ( .NOT. BTEST(wall_flags_0(k, jp, ip), 0) .OR. & 546 .NOT. BTEST(wall_flags_0(kp, j, ip), 0) .OR. & 547 .NOT. BTEST(wall_flags_0(kp, jp, i ), 0) ) & 548 THEN 549 subbox_at_wall = .TRUE. 550 ENDIF 551 ENDIF 552 IF ( subbox_at_wall ) THEN 553 e_int(start_index(nb):end_index(nb)) = e(kp,jp,ip) 554 diss_int(start_index(nb):end_index(nb)) = diss(kp,jp,ip) 555 de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip) 556 de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip) 557 de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip) 558 ! 559 !-- Set flag for stochastic equation. 560 term_1_2(start_index(nb):end_index(nb)) = 0.0_wp 561 ELSE 535 562 i = ip + block_offset(nb)%i_off 536 563 j = jp + block_offset(nb)%j_off … … 600 627 ELSE 601 628 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 602 ( gg - bb ) * de_dy(k+1,j,i+1) + &603 ( gg - cc ) * de_dy(k+1,j+1,i) + &604 ( gg - dd ) * de_dy(k+1,j+1,i+1) &629 ( gg - bb ) * de_dy(k+1,j,i+1) + & 630 ( gg - cc ) * de_dy(k+1,j+1,i) + & 631 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 605 632 ) / ( 3.0_wp * gg ) 606 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &607 ( de_dy_int_u - de_dy_int_l )633 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * & 634 ( de_dy_int_u - de_dy_int_l ) 608 635 ENDIF 609 636 … … 638 665 ( gg - cc ) * diss(k,j+1,i) + & 639 666 ( gg - dd ) * diss(k,j+1,i+1) & 640 ) / ( 3.0_wp * gg )667 ) / ( 3.0_wp * gg ) 641 668 642 669 IF ( k == nzt ) THEN … … 649 676 ) / ( 3.0_wp * gg ) 650 677 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * & 651 ( diss_int_u - diss_int_l )678 ( diss_int_u - diss_int_l ) 652 679 ENDIF 653 680 … … 655 682 !-- Set flag for stochastic equation. 656 683 term_1_2(n) = 1.0_wp 657 658 684 ENDDO 659 ENDDO 660 661 ELSE ! non-flat topography, e.g., buildings 662 663 DO nb = 0, 7 664 665 DO n = start_index(nb), end_index(nb) 666 667 i = particles(n)%x * ddx 668 j = particles(n)%y * ddy 669 k = ( zv(n) + dz * atmos_ocean_sign ) / dz & 670 + offset_ocean_nzt ! only exact if eq.dist 671 ! 672 !-- In case that there are buildings it has to be determined 673 !-- how many of the gridpoints defining the particle box are 674 !-- situated within a building 675 !-- gp_outside_of_building(1): i,j,k 676 !-- gp_outside_of_building(2): i,j+1,k 677 !-- gp_outside_of_building(3): i,j,k+1 678 !-- gp_outside_of_building(4): i,j+1,k+1 679 !-- gp_outside_of_building(5): i+1,j,k 680 !-- gp_outside_of_building(6): i+1,j+1,k 681 !-- gp_outside_of_building(7): i+1,j,k+1 682 !-- gp_outside_of_building(8): i+1,j+1,k+1 683 684 gp_outside_of_building = 0 685 location = 0.0_wp 686 num_gp = 0 687 688 ! 689 !-- Determine vertical index of topography top at (j,i) 690 k_wall = get_topography_top_index( j, i, 's' ) 691 ! 692 !-- To do: Reconsider order of computations in order to avoid 693 !-- unnecessary index calculations. 694 IF ( k > k_wall .OR. k_wall == 0 ) THEN 695 num_gp = num_gp + 1 696 gp_outside_of_building(1) = 1 697 location(num_gp,1) = i * dx 698 location(num_gp,2) = j * dy 699 location(num_gp,3) = k * dz - 0.5_wp * dz 700 ei(num_gp) = e(k,j,i) 701 dissi(num_gp) = diss(k,j,i) 702 de_dxi(num_gp) = de_dx(k,j,i) 703 de_dyi(num_gp) = de_dy(k,j,i) 704 de_dzi(num_gp) = de_dz(k,j,i) 705 ENDIF 706 707 ! 708 !-- Determine vertical index of topography top at (j+1,i) 709 k_wall = get_topography_top_index( j+1, i, 's' ) 710 711 IF ( k > k_wall .OR. k_wall == 0 ) THEN 712 num_gp = num_gp + 1 713 gp_outside_of_building(2) = 1 714 location(num_gp,1) = i * dx 715 location(num_gp,2) = (j+1) * dy 716 location(num_gp,3) = k * dz - 0.5_wp * dz 717 ei(num_gp) = e(k,j+1,i) 718 dissi(num_gp) = diss(k,j+1,i) 719 de_dxi(num_gp) = de_dx(k,j+1,i) 720 de_dyi(num_gp) = de_dy(k,j+1,i) 721 de_dzi(num_gp) = de_dz(k,j+1,i) 722 ENDIF 723 724 ! 725 !-- Determine vertical index of topography top at (j,i) 726 k_wall = get_topography_top_index( j, i, 's' ) 727 728 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 729 num_gp = num_gp + 1 730 gp_outside_of_building(3) = 1 731 location(num_gp,1) = i * dx 732 location(num_gp,2) = j * dy 733 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 734 ei(num_gp) = e(k+1,j,i) 735 dissi(num_gp) = diss(k+1,j,i) 736 de_dxi(num_gp) = de_dx(k+1,j,i) 737 de_dyi(num_gp) = de_dy(k+1,j,i) 738 de_dzi(num_gp) = de_dz(k+1,j,i) 739 ENDIF 740 741 ! 742 !-- Determine vertical index of topography top at (j+1,i) 743 k_wall = get_topography_top_index( j+1, i, 's' ) 744 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 745 num_gp = num_gp + 1 746 gp_outside_of_building(4) = 1 747 location(num_gp,1) = i * dx 748 location(num_gp,2) = (j+1) * dy 749 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 750 ei(num_gp) = e(k+1,j+1,i) 751 dissi(num_gp) = diss(k+1,j+1,i) 752 de_dxi(num_gp) = de_dx(k+1,j+1,i) 753 de_dyi(num_gp) = de_dy(k+1,j+1,i) 754 de_dzi(num_gp) = de_dz(k+1,j+1,i) 755 ENDIF 756 757 ! 758 !-- Determine vertical index of topography top at (j,i+1) 759 k_wall = get_topography_top_index( j, i+1, 's' ) 760 IF ( k > k_wall .OR. k_wall == 0 ) THEN 761 num_gp = num_gp + 1 762 gp_outside_of_building(5) = 1 763 location(num_gp,1) = (i+1) * dx 764 location(num_gp,2) = j * dy 765 location(num_gp,3) = k * dz - 0.5_wp * dz 766 ei(num_gp) = e(k,j,i+1) 767 dissi(num_gp) = diss(k,j,i+1) 768 de_dxi(num_gp) = de_dx(k,j,i+1) 769 de_dyi(num_gp) = de_dy(k,j,i+1) 770 de_dzi(num_gp) = de_dz(k,j,i+1) 771 ENDIF 772 773 ! 774 !-- Determine vertical index of topography top at (j+1,i+1) 775 k_wall = get_topography_top_index( j+1, i+1, 's' ) 776 777 IF ( k > k_wall .OR. k_wall == 0 ) THEN 778 num_gp = num_gp + 1 779 gp_outside_of_building(6) = 1 780 location(num_gp,1) = (i+1) * dx 781 location(num_gp,2) = (j+1) * dy 782 location(num_gp,3) = k * dz - 0.5_wp * dz 783 ei(num_gp) = e(k,j+1,i+1) 784 dissi(num_gp) = diss(k,j+1,i+1) 785 de_dxi(num_gp) = de_dx(k,j+1,i+1) 786 de_dyi(num_gp) = de_dy(k,j+1,i+1) 787 de_dzi(num_gp) = de_dz(k,j+1,i+1) 788 ENDIF 789 790 ! 791 !-- Determine vertical index of topography top at (j,i+1) 792 k_wall = get_topography_top_index( j, i+1, 's' ) 793 794 IF ( k+1 > k_wall .OR. k_wall == 0 ) THEN 795 num_gp = num_gp + 1 796 gp_outside_of_building(7) = 1 797 location(num_gp,1) = (i+1) * dx 798 location(num_gp,2) = j * dy 799 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 800 ei(num_gp) = e(k+1,j,i+1) 801 dissi(num_gp) = diss(k+1,j,i+1) 802 de_dxi(num_gp) = de_dx(k+1,j,i+1) 803 de_dyi(num_gp) = de_dy(k+1,j,i+1) 804 de_dzi(num_gp) = de_dz(k+1,j,i+1) 805 ENDIF 806 807 ! 808 !-- Determine vertical index of topography top at (j+1,i+1) 809 k_wall = get_topography_top_index( j+1, i+1, 's' ) 810 811 IF ( k+1 > k_wall .OR. k_wall == 0) THEN 812 num_gp = num_gp + 1 813 gp_outside_of_building(8) = 1 814 location(num_gp,1) = (i+1) * dx 815 location(num_gp,2) = (j+1) * dy 816 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz 817 ei(num_gp) = e(k+1,j+1,i+1) 818 dissi(num_gp) = diss(k+1,j+1,i+1) 819 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 820 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 821 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 822 ENDIF 823 ! 824 !-- If all gridpoints are situated outside of a building, then the 825 !-- ordinary interpolation scheme can be used. 826 IF ( num_gp == 8 ) THEN 827 828 x = particles(n)%x - i * dx 829 y = particles(n)%y - j * dy 830 aa = x**2 + y**2 831 bb = ( dx - x )**2 + y**2 832 cc = x**2 + ( dy - y )**2 833 dd = ( dx - x )**2 + ( dy - y )**2 834 gg = aa + bb + cc + dd 835 836 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 837 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 838 ) / ( 3.0_wp * gg ) 839 840 IF ( k == nzt ) THEN 841 e_int(n) = e_int_l 842 ELSE 843 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 844 ( gg - bb ) * e(k+1,j,i+1) + & 845 ( gg - cc ) * e(k+1,j+1,i) + & 846 ( gg - dd ) * e(k+1,j+1,i+1) & 847 ) / ( 3.0_wp * gg ) 848 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 849 ( e_int_u - e_int_l ) 850 ENDIF 851 ! 852 !-- Needed to avoid NaN particle velocities (this might not be 853 !-- required any more) 854 IF ( e_int(n) <= 0.0_wp ) THEN 855 e_int(n) = 1.0E-20_wp 856 ENDIF 857 ! 858 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 859 !-- and all position variables from above (TKE)) 860 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 861 ( gg - bb ) * de_dx(k,j,i+1) + & 862 ( gg - cc ) * de_dx(k,j+1,i) + & 863 ( gg - dd ) * de_dx(k,j+1,i+1) & 864 ) / ( 3.0_wp * gg ) 865 866 IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN 867 de_dx_int(n) = de_dx_int_l 868 ELSE 869 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 870 ( gg - bb ) * de_dx(k+1,j,i+1) + & 871 ( gg - cc ) * de_dx(k+1,j+1,i) + & 872 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 873 ) / ( 3.0_wp * gg ) 874 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / & 875 dz * ( de_dx_int_u - de_dx_int_l ) 876 ENDIF 877 878 ! 879 !-- Interpolate the TKE gradient along y 880 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 881 ( gg - bb ) * de_dy(k,j,i+1) + & 882 ( gg - cc ) * de_dy(k,j+1,i) + & 883 ( gg - dd ) * de_dy(k,j+1,i+1) & 884 ) / ( 3.0_wp * gg ) 885 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 886 de_dy_int(n) = de_dy_int_l 887 ELSE 888 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 889 ( gg - bb ) * de_dy(k+1,j,i+1) + & 890 ( gg - cc ) * de_dy(k+1,j+1,i) + & 891 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 892 ) / ( 3.0_wp * gg ) 893 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / & 894 dz * ( de_dy_int_u - de_dy_int_l ) 895 ENDIF 896 897 ! 898 !-- Interpolate the TKE gradient along z 899 IF ( zv(n) < 0.5_wp * dz ) THEN 900 de_dz_int(n) = 0.0_wp 901 ELSE 902 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 903 ( gg - bb ) * de_dz(k,j,i+1) + & 904 ( gg - cc ) * de_dz(k,j+1,i) + & 905 ( gg - dd ) * de_dz(k,j+1,i+1) & 906 ) / ( 3.0_wp * gg ) 907 908 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 909 de_dz_int(n) = de_dz_int_l 910 ELSE 911 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 912 ( gg - bb ) * de_dz(k+1,j,i+1) + & 913 ( gg - cc ) * de_dz(k+1,j+1,i) + & 914 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 915 ) / ( 3.0_wp * gg ) 916 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /& 917 dz * ( de_dz_int_u - de_dz_int_l ) 918 ENDIF 919 ENDIF 920 921 ! 922 !-- Interpolate the dissipation of TKE 923 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 924 ( gg - bb ) * diss(k,j,i+1) + & 925 ( gg - cc ) * diss(k,j+1,i) + & 926 ( gg - dd ) * diss(k,j+1,i+1) & 927 ) / ( 3.0_wp * gg ) 928 929 IF ( k == nzt ) THEN 930 diss_int(n) = diss_int_l 931 ELSE 932 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 933 ( gg - bb ) * diss(k+1,j,i+1) + & 934 ( gg - cc ) * diss(k+1,j+1,i) + & 935 ( gg - dd ) * diss(k+1,j+1,i+1) & 936 ) / ( 3.0_wp * gg ) 937 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *& 938 ( diss_int_u - diss_int_l ) 939 ENDIF 940 ! 941 !-- Set flag for stochastic equation. 942 term_1_2(n) = 1.0_wp 943 944 ELSE 945 946 ! 947 !-- If wall between gridpoint 1 and gridpoint 5, then 948 !-- Neumann boundary condition has to be applied 949 IF ( gp_outside_of_building(1) == 1 .AND. & 950 gp_outside_of_building(5) == 0 ) THEN 951 num_gp = num_gp + 1 952 location(num_gp,1) = i * dx + 0.5_wp * dx 953 location(num_gp,2) = j * dy 954 location(num_gp,3) = k * dz - 0.5_wp * dz 955 ei(num_gp) = e(k,j,i) 956 dissi(num_gp) = diss(k,j,i) 957 de_dxi(num_gp) = 0.0_wp 958 de_dyi(num_gp) = de_dy(k,j,i) 959 de_dzi(num_gp) = de_dz(k,j,i) 960 ENDIF 961 962 IF ( gp_outside_of_building(5) == 1 .AND. & 963 gp_outside_of_building(1) == 0 ) THEN 964 num_gp = num_gp + 1 965 location(num_gp,1) = i * dx + 0.5_wp * dx 966 location(num_gp,2) = j * 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) = 0.0_wp 971 de_dyi(num_gp) = de_dy(k,j,i+1) 972 de_dzi(num_gp) = de_dz(k,j,i+1) 973 ENDIF 974 975 ! 976 !-- If wall between gridpoint 5 and gridpoint 6, then 977 !-- then Neumann boundary condition has to be applied 978 IF ( gp_outside_of_building(5) == 1 .AND. & 979 gp_outside_of_building(6) == 0 ) THEN 980 num_gp = num_gp + 1 981 location(num_gp,1) = (i+1) * dx 982 location(num_gp,2) = j * dy + 0.5_wp * dy 983 location(num_gp,3) = k * dz - 0.5_wp * dz 984 ei(num_gp) = e(k,j,i+1) 985 dissi(num_gp) = diss(k,j,i+1) 986 de_dxi(num_gp) = de_dx(k,j,i+1) 987 de_dyi(num_gp) = 0.0_wp 988 de_dzi(num_gp) = de_dz(k,j,i+1) 989 ENDIF 990 991 IF ( gp_outside_of_building(6) == 1 .AND. & 992 gp_outside_of_building(5) == 0 ) THEN 993 num_gp = num_gp + 1 994 location(num_gp,1) = (i+1) * dx 995 location(num_gp,2) = j * dy + 0.5_wp * dy 996 location(num_gp,3) = k * dz - 0.5_wp * dz 997 ei(num_gp) = e(k,j+1,i+1) 998 dissi(num_gp) = diss(k,j+1,i+1) 999 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1000 de_dyi(num_gp) = 0.0_wp 1001 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1002 ENDIF 1003 1004 ! 1005 !-- If wall between gridpoint 2 and gridpoint 6, then 1006 !-- Neumann boundary condition has to be applied 1007 IF ( gp_outside_of_building(2) == 1 .AND. & 1008 gp_outside_of_building(6) == 0 ) THEN 1009 num_gp = num_gp + 1 1010 location(num_gp,1) = i * dx + 0.5_wp * dx 1011 location(num_gp,2) = (j+1) * dy 1012 location(num_gp,3) = k * dz - 0.5_wp * dz 1013 ei(num_gp) = e(k,j+1,i) 1014 dissi(num_gp) = diss(k,j+1,i) 1015 de_dxi(num_gp) = 0.0_wp 1016 de_dyi(num_gp) = de_dy(k,j+1,i) 1017 de_dzi(num_gp) = de_dz(k,j+1,i) 1018 ENDIF 1019 1020 IF ( gp_outside_of_building(6) == 1 .AND. & 1021 gp_outside_of_building(2) == 0 ) THEN 1022 num_gp = num_gp + 1 1023 location(num_gp,1) = i * dx + 0.5_wp * dx 1024 location(num_gp,2) = (j+1) * dy 1025 location(num_gp,3) = k * dz - 0.5_wp * dz 1026 ei(num_gp) = e(k,j+1,i+1) 1027 dissi(num_gp) = diss(k,j+1,i+1) 1028 de_dxi(num_gp) = 0.0_wp 1029 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1030 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1031 ENDIF 1032 1033 ! 1034 !-- If wall between gridpoint 1 and gridpoint 2, then 1035 !-- Neumann boundary condition has to be applied 1036 IF ( gp_outside_of_building(1) == 1 .AND. & 1037 gp_outside_of_building(2) == 0 ) THEN 1038 num_gp = num_gp + 1 1039 location(num_gp,1) = i * dx 1040 location(num_gp,2) = j * dy + 0.5_wp * dy 1041 location(num_gp,3) = k * dz - 0.5_wp * dz 1042 ei(num_gp) = e(k,j,i) 1043 dissi(num_gp) = diss(k,j,i) 1044 de_dxi(num_gp) = de_dx(k,j,i) 1045 de_dyi(num_gp) = 0.0_wp 1046 de_dzi(num_gp) = de_dz(k,j,i) 1047 ENDIF 1048 1049 IF ( gp_outside_of_building(2) == 1 .AND. & 1050 gp_outside_of_building(1) == 0 ) THEN 1051 num_gp = num_gp + 1 1052 location(num_gp,1) = i * dx 1053 location(num_gp,2) = j * dy + 0.5_wp * dy 1054 location(num_gp,3) = k * dz - 0.5_wp * dz 1055 ei(num_gp) = e(k,j+1,i) 1056 dissi(num_gp) = diss(k,j+1,i) 1057 de_dxi(num_gp) = de_dx(k,j+1,i) 1058 de_dyi(num_gp) = 0.0_wp 1059 de_dzi(num_gp) = de_dz(k,j+1,i) 1060 ENDIF 1061 1062 ! 1063 !-- If wall between gridpoint 3 and gridpoint 7, then 1064 !-- Neumann boundary condition has to be applied 1065 IF ( gp_outside_of_building(3) == 1 .AND. & 1066 gp_outside_of_building(7) == 0 ) THEN 1067 num_gp = num_gp + 1 1068 location(num_gp,1) = i * dx + 0.5_wp * dx 1069 location(num_gp,2) = j * dy 1070 location(num_gp,3) = k * dz + 0.5_wp * dz 1071 ei(num_gp) = e(k+1,j,i) 1072 dissi(num_gp) = diss(k+1,j,i) 1073 de_dxi(num_gp) = 0.0_wp 1074 de_dyi(num_gp) = de_dy(k+1,j,i) 1075 de_dzi(num_gp) = de_dz(k+1,j,i) 1076 ENDIF 1077 1078 IF ( gp_outside_of_building(7) == 1 .AND. & 1079 gp_outside_of_building(3) == 0 ) THEN 1080 num_gp = num_gp + 1 1081 location(num_gp,1) = i * dx + 0.5_wp * dx 1082 location(num_gp,2) = j * 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) = 0.0_wp 1087 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1088 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1089 ENDIF 1090 1091 ! 1092 !-- If wall between gridpoint 7 and gridpoint 8, then 1093 !-- Neumann boundary condition has to be applied 1094 IF ( gp_outside_of_building(7) == 1 .AND. & 1095 gp_outside_of_building(8) == 0 ) THEN 1096 num_gp = num_gp + 1 1097 location(num_gp,1) = (i+1) * dx 1098 location(num_gp,2) = j * dy + 0.5_wp * dy 1099 location(num_gp,3) = k * dz + 0.5_wp * dz 1100 ei(num_gp) = e(k+1,j,i+1) 1101 dissi(num_gp) = diss(k+1,j,i+1) 1102 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1103 de_dyi(num_gp) = 0.0_wp 1104 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1105 ENDIF 1106 1107 IF ( gp_outside_of_building(8) == 1 .AND. & 1108 gp_outside_of_building(7) == 0 ) THEN 1109 num_gp = num_gp + 1 1110 location(num_gp,1) = (i+1) * dx 1111 location(num_gp,2) = j * dy + 0.5_wp * dy 1112 location(num_gp,3) = k * dz + 0.5_wp * dz 1113 ei(num_gp) = e(k+1,j+1,i+1) 1114 dissi(num_gp) = diss(k+1,j+1,i+1) 1115 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1116 de_dyi(num_gp) = 0.0_wp 1117 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1118 ENDIF 1119 1120 ! 1121 !-- If wall between gridpoint 4 and gridpoint 8, then 1122 !-- Neumann boundary condition has to be applied 1123 IF ( gp_outside_of_building(4) == 1 .AND. & 1124 gp_outside_of_building(8) == 0 ) THEN 1125 num_gp = num_gp + 1 1126 location(num_gp,1) = i * dx + 0.5_wp * dx 1127 location(num_gp,2) = (j+1) * dy 1128 location(num_gp,3) = k * dz + 0.5_wp * dz 1129 ei(num_gp) = e(k+1,j+1,i) 1130 dissi(num_gp) = diss(k+1,j+1,i) 1131 de_dxi(num_gp) = 0.0_wp 1132 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1133 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1134 ENDIF 1135 1136 IF ( gp_outside_of_building(8) == 1 .AND. & 1137 gp_outside_of_building(4) == 0 ) THEN 1138 num_gp = num_gp + 1 1139 location(num_gp,1) = i * dx + 0.5_wp * dx 1140 location(num_gp,2) = (j+1) * dy 1141 location(num_gp,3) = k * dz + 0.5_wp * dz 1142 ei(num_gp) = e(k+1,j+1,i+1) 1143 dissi(num_gp) = diss(k+1,j+1,i+1) 1144 de_dxi(num_gp) = 0.0_wp 1145 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1146 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1147 ENDIF 1148 1149 ! 1150 !-- If wall between gridpoint 3 and gridpoint 4, then 1151 !-- Neumann boundary condition has to be applied 1152 IF ( gp_outside_of_building(3) == 1 .AND. & 1153 gp_outside_of_building(4) == 0 ) THEN 1154 num_gp = num_gp + 1 1155 location(num_gp,1) = i * dx 1156 location(num_gp,2) = j * dy + 0.5_wp * dy 1157 location(num_gp,3) = k * dz + 0.5_wp * dz 1158 ei(num_gp) = e(k+1,j,i) 1159 dissi(num_gp) = diss(k+1,j,i) 1160 de_dxi(num_gp) = de_dx(k+1,j,i) 1161 de_dyi(num_gp) = 0.0_wp 1162 de_dzi(num_gp) = de_dz(k+1,j,i) 1163 ENDIF 1164 1165 IF ( gp_outside_of_building(4) == 1 .AND. & 1166 gp_outside_of_building(3) == 0 ) THEN 1167 num_gp = num_gp + 1 1168 location(num_gp,1) = i * dx 1169 location(num_gp,2) = j * dy + 0.5_wp * dy 1170 location(num_gp,3) = k * dz + 0.5_wp * dz 1171 ei(num_gp) = e(k+1,j+1,i) 1172 dissi(num_gp) = diss(k+1,j+1,i) 1173 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1174 de_dyi(num_gp) = 0.0_wp 1175 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1176 ENDIF 1177 1178 ! 1179 !-- If wall between gridpoint 1 and gridpoint 3, then 1180 !-- Neumann boundary condition has to be applied 1181 !-- (only one case as only building beneath is possible) 1182 IF ( gp_outside_of_building(1) == 0 .AND. & 1183 gp_outside_of_building(3) == 1 ) THEN 1184 num_gp = num_gp + 1 1185 location(num_gp,1) = i * dx 1186 location(num_gp,2) = j * dy 1187 location(num_gp,3) = k * dz 1188 ei(num_gp) = e(k+1,j,i) 1189 dissi(num_gp) = diss(k+1,j,i) 1190 de_dxi(num_gp) = de_dx(k+1,j,i) 1191 de_dyi(num_gp) = de_dy(k+1,j,i) 1192 de_dzi(num_gp) = 0.0_wp 1193 ENDIF 1194 1195 ! 1196 !-- If wall between gridpoint 5 and gridpoint 7, then 1197 !-- Neumann boundary condition has to be applied 1198 !-- (only one case as only building beneath is possible) 1199 IF ( gp_outside_of_building(5) == 0 .AND. & 1200 gp_outside_of_building(7) == 1 ) THEN 1201 num_gp = num_gp + 1 1202 location(num_gp,1) = (i+1) * dx 1203 location(num_gp,2) = j * dy 1204 location(num_gp,3) = k * dz 1205 ei(num_gp) = e(k+1,j,i+1) 1206 dissi(num_gp) = diss(k+1,j,i+1) 1207 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1208 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1209 de_dzi(num_gp) = 0.0_wp 1210 ENDIF 1211 1212 ! 1213 !-- If wall between gridpoint 2 and gridpoint 4, then 1214 !-- Neumann boundary condition has to be applied 1215 !-- (only one case as only building beneath is possible) 1216 IF ( gp_outside_of_building(2) == 0 .AND. & 1217 gp_outside_of_building(4) == 1 ) THEN 1218 num_gp = num_gp + 1 1219 location(num_gp,1) = i * dx 1220 location(num_gp,2) = (j+1) * dy 1221 location(num_gp,3) = k * dz 1222 ei(num_gp) = e(k+1,j+1,i) 1223 dissi(num_gp) = diss(k+1,j+1,i) 1224 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1225 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1226 de_dzi(num_gp) = 0.0_wp 1227 ENDIF 1228 1229 ! 1230 !-- If wall between gridpoint 6 and gridpoint 8, then 1231 !-- Neumann boundary condition has to be applied 1232 !-- (only one case as only building beneath is possible) 1233 IF ( gp_outside_of_building(6) == 0 .AND. & 1234 gp_outside_of_building(8) == 1 ) THEN 1235 num_gp = num_gp + 1 1236 location(num_gp,1) = (i+1) * dx 1237 location(num_gp,2) = (j+1) * dy 1238 location(num_gp,3) = k * dz 1239 ei(num_gp) = e(k+1,j+1,i+1) 1240 dissi(num_gp) = diss(k+1,j+1,i+1) 1241 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1242 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1243 de_dzi(num_gp) = 0.0_wp 1244 ENDIF 1245 1246 ! 1247 !-- Carry out the interpolation 1248 IF ( num_gp == 1 ) THEN 1249 ! 1250 !-- If only one of the gridpoints is situated outside of the 1251 !-- building, it follows that the values at the particle 1252 !-- location are the same as the gridpoint values 1253 e_int(n) = ei(num_gp) 1254 diss_int(n) = dissi(num_gp) 1255 de_dx_int(n) = de_dxi(num_gp) 1256 de_dy_int(n) = de_dyi(num_gp) 1257 de_dz_int(n) = de_dzi(num_gp) 1258 ! 1259 !-- Set flag for stochastic equation which disables calculation 1260 !-- of drift and memory term near topography. 1261 term_1_2(n) = 0.0_wp 1262 ELSE IF ( num_gp > 1 ) THEN 1263 1264 d_sum = 0.0_wp 1265 ! 1266 !-- Evaluation of the distances between the gridpoints 1267 !-- contributing to the interpolated values, and the particle 1268 !-- location 1269 DO agp = 1, num_gp 1270 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1271 + ( particles(n)%y-location(agp,2) )**2 & 1272 + ( zv(n)-location(agp,3) )**2 1273 d_sum = d_sum + d_gp_pl(agp) 1274 ENDDO 1275 1276 ! 1277 !-- Finally the interpolation can be carried out 1278 e_int(n) = 0.0_wp 1279 diss_int(n) = 0.0_wp 1280 de_dx_int(n) = 0.0_wp 1281 de_dy_int(n) = 0.0_wp 1282 de_dz_int(n) = 0.0_wp 1283 DO agp = 1, num_gp 1284 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1285 ei(agp) / ( (num_gp-1) * d_sum ) 1286 diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1287 dissi(agp) / ( (num_gp-1) * d_sum ) 1288 de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1289 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1290 de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1291 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1292 de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * & 1293 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1294 ENDDO 1295 1296 ENDIF 1297 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1298 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1299 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1300 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1301 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1302 ! 1303 !-- Set flag for stochastic equation which disables calculation 1304 !-- of drift and memory term near topography. 1305 term_1_2(n) = 0.0_wp 1306 ENDIF 1307 ENDDO 1308 ENDDO 1309 ENDIF 685 ENDIF 686 ENDDO 1310 687 1311 688 DO nb = 0,7
Note: See TracChangeset
for help on using the changeset viewer.
