Changeset 57 for palm/trunk/SOURCE/advec_particles.f90
- Timestamp:
- Mar 9, 2007 12:05:41 PM (17 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_particles.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_particles.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Particle reflection at vertical walls (by Jin Zhang), regarding vertical 7 ! walls in the SGS model, 8 ! + user interface user_advec_particles 6 9 ! TEST: PRINT statements on unit 9 (commented out) 7 10 ! … … 37 40 ! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z, 38 41 ! izuf renamed iran, output of particle time series 39 !40 ! Revision 1.28 2006/02/23 09:41:21 raasch41 ! Only a fraction of the particles may have tails, therefore output for unit42 ! 85 and 90 has changed,43 ! sorting of particles in gridbox order after each timestep, optimization of44 ! the droplet growth by collision,45 ! 0.5 replaced by 0.4999999999 in setting the random fluctuations,46 ! psl, psr, pdx, etc. are now 1d arrays, i.e. they depend on the particle group,47 ! improve particle release at PE boundaries, nt_anz renamed48 ! current_timestep_number, idum (particle_type) renamed tail_id,49 ! error number argument for handle_netcdf_error,50 !51 ! Revision 1.27 2005/10/20 14:04:15 raasch52 ! droplet growth by collision completed (still has to be optimized for speed),53 ! new routine collision_efficiency (see end of this file),54 ! 2*r replaced by r in the exponential terms (momentum equation for particles),55 ! number of particles really used is additionally output on the netcdf particle56 ! data file,57 ! test print statement removed58 !59 ! Revision 1.26 2005/06/26 19:50:14 raasch60 ! Cloud droplet features added, radius is used instead of diameter,61 ! particle type initial data assignment extended,62 ! output on unit 80 for non-parallel runs63 !64 ! Revision 1.25 2005/05/18 14:53:11 raasch65 ! Extensions for NetCDF output66 !67 ! Revision 1.24 2005/04/23 08:24:58 raasch68 ! Revised calculation of output time counters regarding a possible decrease of69 ! the output time interval in case of restart runs70 !71 ! Revision 1.23 2005/03/26 14:20:05 raasch72 ! calculation of vertical particle velocity (with inertia) corrected, exp_arg73 ! had a wrong sign74 !75 ! Revision 1.22 2004/04/30 07:48:57 raasch76 ! Particle type initial data assignment extended77 !78 ! Revision 1.21 2004/01/28 15:00:48 raasch79 ! Output of particle infos in subroutine allocate_prt_memory on demand only80 !81 ! Revision 1.20 2003/10/29 08:38:52 raasch82 ! Module random_function_mod is used, modifications for new particle group83 ! feature, version number is written on binary file84 !85 ! Revision 1.19 2003/03/16 09:19:32 raasch86 ! Two underscores (_) are placed in front of all define-strings87 !88 ! Revision 1.18 2003/03/14 13:36:31 raasch89 ! Error in reflection boundary condition removed: particle speed must also90 ! be inverted91 !92 ! Revision 1.17 2003/03/04 11:23:20 raasch93 ! Error in particle inertia part removed (exp_arg must not contain the timestep)94 ! Particle velocities are also stored in particles in case of zero density ratio95 !96 ! Revision 1.16 2002/09/12 12:55:29 raasch97 ! Transport of particles with inertia implemented, write density_ratio to98 ! restart file99 !100 ! Revision 1.15 2002/06/11 12:11:56 raasch101 ! Declaration of u_int_l, u_int_u, v_int_l, v_int_u added102 !103 ! Revision 1.14 2002/05/02 18:46:34 18:46:34 raasch (Siegfried Raasch)104 ! Horizontal velocity components for particle advection are now interpolated105 ! between horizontal grid levels106 !107 ! Revision 1.13 2002/04/24 18:46:11 18:46:11 raasch (Siegfried Raasch)108 ! Temporary particle-arrays trlp, trnp, trrp and trsp are now allocatable arrays109 !110 ! Revision 1.12 2002/04/16 07:44:07 raasch111 ! Additional boundary conditions added, particle data for later analysis can be112 ! written on file113 !114 ! Revision 1.11 2001/11/09 13:38:20 raasch115 ! Colour information for particle tails is stored in array116 ! particle_tail_coordinates, adding information to the particle tails117 ! moved after call of user_particle_attributes118 !119 ! Revision 1.10 2001/08/21 08:19:48 raasch120 ! Storage of particle tails, precision of vertical advection improved,121 ! particle attributes (e.g. size and color) can be defined by the user122 !123 ! Revision 1.9 2001/07/12 12:01:27 raasch124 ! Random start positions of initial particles possible125 !126 ! Revision 1.8 2001/03/29 17:27:10 raasch127 ! Translation of remaining German identifiers (variables, subroutines, etc.)128 !129 ! Revision 1.7 2001/02/13 10:24:00 raasch130 ! Particle advection possible in case of galilei transformation131 !132 ! Revision 1.6 2001/01/25 06:51:04 raasch133 ! Module test_variables removed, writing of particle informations is optional134 ! now135 !136 ! Revision 1.5 2001/01/02 17:15:29 raasch137 ! Unit 80 is immediately closed after usage every time. Particle positions138 ! are written on unit 90 instead of 81 (subroutine write_particles).139 !140 ! Revision 1.4 2000/12/28 12:43:54 raasch141 ! Packing of particles-array is done explicitly (not by PACK-function, which142 ! seems to have problems when a large number of PEs is used)143 ! Missing translations included,144 ! code is used only optionally (cpp-directives are added)145 !146 ! Revision 1.3 2000/01/20 08:51:37 letzel147 ! All comments translated into English148 !149 ! Revision 1.2 1999/11/25 16:41:48 raasch150 ! Indexfehler im nichtparallelen Teil korrigiert151 42 ! 152 43 ! Revision 1.1 1999/11/25 16:16:06 raasch … … 693 584 DO j = nys, nyn 694 585 DO k = nzb, nzt+1 695 de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 696 de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 586 587 IF ( k <= nzb_s_inner(j,i-1) .AND. & 588 k > nzb_s_inner(j,i) .AND. & 589 k > nzb_s_inner(j,i+1) ) THEN 590 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 591 ( e(k,j,i+1) - e(k,j,i) ) * ddx 592 ELSEIF ( k > nzb_s_inner(j,i-1) .AND. & 593 k > nzb_s_inner(j,i) .AND. & 594 k <= nzb_s_inner(j,i+1) ) THEN 595 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 596 ( e(k,j,i) - e(k,j,i-1) ) * ddx 597 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) & 598 THEN 599 de_dx(k,j,i) = 0.0 600 ELSEIF ( k < nzb_s_inner(j,i-1) .AND. k < nzb_s_inner(j,i) ) & 601 THEN 602 de_dx(k,j,i) = 0.0 603 ELSE 604 de_dx(k,j,i) = sgs_wfu_part * & 605 ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 606 ENDIF 607 608 IF ( k <= nzb_s_inner(j-1,i) .AND. & 609 k > nzb_s_inner(j,i) .AND. & 610 k > nzb_s_inner(j+1,i) ) THEN 611 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 612 ( e(k,j+1,i) - e(k,j,i) ) * ddy 613 ELSEIF ( k > nzb_s_inner(j-1,i) .AND. & 614 k > nzb_s_inner(j,i) .AND. & 615 k <= nzb_s_inner(j+1,i) ) THEN 616 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 617 ( e(k,j,i) - e(k,j-1,i) ) * ddy 618 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) & 619 THEN 620 de_dy(k,j,i) = 0.0 621 ELSEIF ( k < nzb_s_inner(j-1,i) .AND. k < nzb_s_inner(j,i) ) & 622 THEN 623 de_dy(k,j,i) = 0.0 624 ELSE 625 de_dy(k,j,i) = sgs_wfv_part * & 626 ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 627 ENDIF 628 697 629 ENDDO 698 630 ENDDO … … 703 635 DO i = nxl, nxr 704 636 DO j = nys, nyn 705 DO k = nzb+2, nzt-1 637 638 DO k = nzb_s_inner(j,i)+2, nzt-1 706 639 de_dz(k,j,i) = 2.0 * sgs_wfw_part * & 707 640 ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) ) 708 641 ENDDO 709 de_dz(nzb,j,i) = 0.0 710 de_dz(nzb+1,j,i) = 2.0 * sgs_wfw_part * & 711 ( e(nzb+2,j,i) - e(nzb+1,j,i) ) / & 712 ( zu(nzb+2) - zu(nzb+1) ) 642 643 k = nzb_s_inner(j,i) 644 de_dz(nzb:k,j,i) = 0.0 645 de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) & 646 / ( zu(k+2) - zu(k+1) ) 713 647 de_dz(nzt,j,i) = 0.0 714 648 de_dz(nzt+1,j,i) = 0.0 715 649 ENDDO 716 ENDDO 650 ENDDO 717 651 718 652 ! … … 967 901 k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if eq.dist 968 902 969 x = particles(n)%x - i * dx 970 y = particles(n)%y - j * dy 971 aa = x**2 + y**2 972 bb = ( dx - x )**2 + y**2 973 cc = x**2 + ( dy - y )**2 974 dd = ( dx - x )**2 + ( dy - y )**2 975 gg = aa + bb + cc + dd 976 977 e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 978 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) & 979 ) / ( 3.0 * gg ) 980 981 IF ( k+1 == nzt+1 ) THEN 982 e_int = e_int_l 903 IF ( topography == 'flat' ) THEN 904 905 x = particles(n)%x - i * dx 906 y = particles(n)%y - j * dy 907 aa = x**2 + y**2 908 bb = ( dx - x )**2 + y**2 909 cc = x**2 + ( dy - y )**2 910 dd = ( dx - x )**2 + ( dy - y )**2 911 gg = aa + bb + cc + dd 912 913 e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 914 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) & 915 ) / ( 3.0 * gg ) 916 917 IF ( k+1 == nzt+1 ) THEN 918 e_int = e_int_l 919 ELSE 920 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 921 ( gg - bb ) * e(k+1,j,i+1) + & 922 ( gg - cc ) * e(k+1,j+1,i) + & 923 ( gg - dd ) * e(k+1,j+1,i+1) & 924 ) / ( 3.0 * gg ) 925 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 926 ( e_int_u - e_int_l ) 927 ENDIF 928 929 ! 930 !-- Interpolate the TKE gradient along x (adopt incides i,j,k and 931 !-- all position variables from above (TKE)) 932 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 933 ( gg - bb ) * de_dx(k,j,i+1) + & 934 ( gg - cc ) * de_dx(k,j+1,i) + & 935 ( gg - dd ) * de_dx(k,j+1,i+1) & 936 ) / ( 3.0 * gg ) 937 938 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 939 de_dx_int = de_dx_int_l 940 ELSE 941 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 942 ( gg - bb ) * de_dx(k+1,j,i+1) + & 943 ( gg - cc ) * de_dx(k+1,j+1,i) + & 944 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 945 ) / ( 3.0 * gg ) 946 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * & 947 ( de_dx_int_u - de_dx_int_l ) 948 ENDIF 949 950 ! 951 !-- Interpolate the TKE gradient along y 952 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 953 ( gg - bb ) * de_dy(k,j,i+1) + & 954 ( gg - cc ) * de_dy(k,j+1,i) + & 955 ( gg - dd ) * de_dy(k,j+1,i+1) & 956 ) / ( 3.0 * gg ) 957 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 958 de_dy_int = de_dy_int_l 959 ELSE 960 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 961 ( gg - bb ) * de_dy(k+1,j,i+1) + & 962 ( gg - cc ) * de_dy(k+1,j+1,i) + & 963 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 964 ) / ( 3.0 * gg ) 965 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * & 966 ( de_dy_int_u - de_dy_int_l ) 967 ENDIF 968 969 ! 970 !-- Interpolate the TKE gradient along z 971 IF ( particles(n)%z < 0.5 * dz ) THEN 972 de_dz_int = 0.0 973 ELSE 974 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 975 ( gg - bb ) * de_dz(k,j,i+1) + & 976 ( gg - cc ) * de_dz(k,j+1,i) + & 977 ( gg - dd ) * de_dz(k,j+1,i+1) & 978 ) / ( 3.0 * gg ) 979 980 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 981 de_dz_int = de_dz_int_l 982 ELSE 983 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 984 ( gg - bb ) * de_dz(k+1,j,i+1) + & 985 ( gg - cc ) * de_dz(k+1,j+1,i) + & 986 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 987 ) / ( 3.0 * gg ) 988 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * & 989 ( de_dz_int_u - de_dz_int_l ) 990 ENDIF 991 ENDIF 992 993 ! 994 !-- Interpolate the dissipation of TKE 995 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 996 ( gg - bb ) * diss(k,j,i+1) + & 997 ( gg - cc ) * diss(k,j+1,i) + & 998 ( gg - dd ) * diss(k,j+1,i+1) & 999 ) / ( 3.0 * gg ) 1000 1001 IF ( k+1 == nzt+1 ) THEN 1002 diss_int = diss_int_l 1003 ELSE 1004 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1005 ( gg - bb ) * diss(k+1,j,i+1) + & 1006 ( gg - cc ) * diss(k+1,j+1,i) + & 1007 ( gg - dd ) * diss(k+1,j+1,i+1) & 1008 ) / ( 3.0 * gg ) 1009 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * & 1010 ( diss_int_u - diss_int_l ) 1011 ENDIF 1012 983 1013 ELSE 984 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 985 ( gg - bb ) * e(k+1,j,i+1) + & 986 ( gg - cc ) * e(k+1,j+1,i) + & 987 ( gg - dd ) * e(k+1,j+1,i+1) & 988 ) / ( 3.0 * gg ) 989 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 990 ( e_int_u - e_int_l ) 991 ENDIF 992 993 ! 994 !-- Interpolate the TKE gradient along x (adopt incides i,j,k and all 995 !-- position variables from above (TKE)) 996 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 997 ( gg - bb ) * de_dx(k,j,i+1) + & 998 ( gg - cc ) * de_dx(k,j+1,i) + & 999 ( gg - dd ) * de_dx(k,j+1,i+1) & 1000 ) / ( 3.0 * gg ) 1001 1002 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1003 de_dx_int = de_dx_int_l 1004 ELSE 1005 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1006 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1007 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1008 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1009 ) / ( 3.0 * gg ) 1010 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * & 1011 ( de_dx_int_u - de_dx_int_l ) 1012 ENDIF 1013 1014 ! 1015 !-- Interpolate the TKE gradient along y 1016 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1017 ( gg - bb ) * de_dy(k,j,i+1) + & 1018 ( gg - cc ) * de_dy(k,j+1,i) + & 1019 ( gg - dd ) * de_dy(k,j+1,i+1) & 1020 ) / ( 3.0 * gg ) 1021 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1022 de_dy_int = de_dy_int_l 1023 ELSE 1024 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1025 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1026 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1027 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1028 ) / ( 3.0 * gg ) 1029 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * & 1030 ( de_dy_int_u - de_dy_int_l ) 1031 ENDIF 1032 1033 ! 1034 !-- Interpolate the TKE gradient along z 1035 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1036 ( gg - bb ) * de_dz(k,j,i+1) + & 1037 ( gg - cc ) * de_dz(k,j+1,i) + & 1038 ( gg - dd ) * de_dz(k,j+1,i+1) & 1039 ) / ( 3.0 * gg ) 1040 1041 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1042 de_dz_int = de_dz_int_l 1043 ELSE 1044 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1045 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1046 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1047 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1048 ) / ( 3.0 * gg ) 1049 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * & 1050 ( de_dz_int_u - de_dz_int_l ) 1051 ENDIF 1052 1053 ! 1054 !-- Interpolate the dissipation of TKE 1055 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1056 ( gg - bb ) * diss(k,j,i+1) + & 1057 ( gg - cc ) * diss(k,j+1,i) + & 1058 ( gg - dd ) * diss(k,j+1,i+1) & 1059 ) / ( 3.0 * gg ) 1060 1061 IF ( k+1 == nzt+1 ) THEN 1062 diss_int = diss_int_l 1063 ELSE 1064 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1065 ( gg - bb ) * diss(k+1,j,i+1) + & 1066 ( gg - cc ) * diss(k+1,j+1,i) + & 1067 ( gg - dd ) * diss(k+1,j+1,i+1) & 1068 ) / ( 3.0 * gg ) 1069 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * & 1070 ( diss_int_u - diss_int_l ) 1071 ENDIF 1014 1015 ! 1016 !-- In case that there are buildings it has to be determined 1017 !-- how many of the gridpoints defining the particle box are 1018 !-- situated within a building 1019 !-- gp_outside_of_building(1): i,j,k 1020 !-- gp_outside_of_building(2): i,j+1,k 1021 !-- gp_outside_of_building(3): i,j,k+1 1022 !-- gp_outside_of_building(4): i,j+1,k+1 1023 !-- gp_outside_of_building(5): i+1,j,k 1024 !-- gp_outside_of_building(6): i+1,j+1,k 1025 !-- gp_outside_of_building(7): i+1,j,k+1 1026 !-- gp_outside_of_building(8): i+1,j+1,k+1 1027 1028 gp_outside_of_building = 0 1029 location = 0.0 1030 num_gp = 0 1031 1032 IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1033 num_gp = num_gp + 1 1034 gp_outside_of_building(1) = 1 1035 location(num_gp,1) = i * dx 1036 location(num_gp,2) = j * dx 1037 location(num_gp,3) = k * dz - 0.5 * dz 1038 ei(num_gp) = e(k,j,i) 1039 dissi(num_gp) = diss(k,j,i) 1040 de_dxi(num_gp) = de_dx(k,j,i) 1041 de_dyi(num_gp) = de_dy(k,j,i) 1042 de_dzi(num_gp) = de_dz(k,j,i) 1043 ENDIF 1044 1045 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1046 THEN 1047 num_gp = num_gp + 1 1048 gp_outside_of_building(2) = 1 1049 location(num_gp,1) = i * dx 1050 location(num_gp,2) = (j+1) * dx 1051 location(num_gp,3) = k * dz - 0.5 * dz 1052 ei(num_gp) = e(k,j+1,i) 1053 dissi(num_gp) = diss(k,j+1,i) 1054 de_dxi(num_gp) = de_dx(k,j+1,i) 1055 de_dyi(num_gp) = de_dy(k,j+1,i) 1056 de_dzi(num_gp) = de_dz(k,j+1,i) 1057 ENDIF 1058 1059 IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1060 num_gp = num_gp + 1 1061 gp_outside_of_building(3) = 1 1062 location(num_gp,1) = i * dx 1063 location(num_gp,2) = j * dy 1064 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1065 ei(num_gp) = e(k+1,j,i) 1066 dissi(num_gp) = diss(k+1,j,i) 1067 de_dxi(num_gp) = de_dx(k+1,j,i) 1068 de_dyi(num_gp) = de_dy(k+1,j,i) 1069 de_dzi(num_gp) = de_dz(k+1,j,i) 1070 ENDIF 1071 1072 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1073 THEN 1074 num_gp = num_gp + 1 1075 gp_outside_of_building(4) = 1 1076 location(num_gp,1) = i * dx 1077 location(num_gp,2) = (j+1) * dy 1078 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1079 ei(num_gp) = e(k+1,j+1,i) 1080 dissi(num_gp) = diss(k+1,j+1,i) 1081 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1082 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1083 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1084 ENDIF 1085 1086 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1087 THEN 1088 num_gp = num_gp + 1 1089 gp_outside_of_building(5) = 1 1090 location(num_gp,1) = (i+1) * dx 1091 location(num_gp,2) = j * dy 1092 location(num_gp,3) = k * dz - 0.5 * dz 1093 ei(num_gp) = e(k,j,i+1) 1094 dissi(num_gp) = diss(k,j,i+1) 1095 de_dxi(num_gp) = de_dx(k,j,i+1) 1096 de_dyi(num_gp) = de_dy(k,j,i+1) 1097 de_dzi(num_gp) = de_dz(k,j,i+1) 1098 ENDIF 1099 1100 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) & 1101 THEN 1102 num_gp = num_gp + 1 1103 gp_outside_of_building(6) = 1 1104 location(num_gp,1) = (i+1) * dx 1105 location(num_gp,2) = (j+1) * dy 1106 location(num_gp,3) = k * dz - 0.5 * dz 1107 ei(num_gp) = e(k,j+1,i+1) 1108 dissi(num_gp) = diss(k,j+1,i+1) 1109 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1110 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1111 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1112 ENDIF 1113 1114 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1115 THEN 1116 num_gp = num_gp + 1 1117 gp_outside_of_building(7) = 1 1118 location(num_gp,1) = (i+1) * dx 1119 location(num_gp,2) = j * dy 1120 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1121 ei(num_gp) = e(k+1,j,i+1) 1122 dissi(num_gp) = diss(k+1,j,i+1) 1123 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1124 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1125 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1126 ENDIF 1127 1128 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)& 1129 THEN 1130 num_gp = num_gp + 1 1131 gp_outside_of_building(8) = 1 1132 location(num_gp,1) = (i+1) * dx 1133 location(num_gp,2) = (j+1) * dy 1134 location(num_gp,3) = (k+1) * dz - 0.5 * 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) = de_dx(k+1,j+1,i+1) 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 all gridpoints are situated outside of a building, then the 1144 !-- ordinary interpolation scheme can be used. 1145 IF ( num_gp == 8 ) THEN 1146 1147 x = particles(n)%x - i * dx 1148 y = particles(n)%y - j * dy 1149 aa = x**2 + y**2 1150 bb = ( dx - x )**2 + y**2 1151 cc = x**2 + ( dy - y )**2 1152 dd = ( dx - x )**2 + ( dy - y )**2 1153 gg = aa + bb + cc + dd 1154 1155 e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 1156 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)& 1157 ) / ( 3.0 * gg ) 1158 1159 IF ( k+1 == nzt+1 ) THEN 1160 e_int = e_int_l 1161 ELSE 1162 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 1163 ( gg - bb ) * e(k+1,j,i+1) + & 1164 ( gg - cc ) * e(k+1,j+1,i) + & 1165 ( gg - dd ) * e(k+1,j+1,i+1) & 1166 ) / ( 3.0 * gg ) 1167 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 1168 ( e_int_u - e_int_l ) 1169 ENDIF 1170 1171 ! 1172 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 1173 !-- and all position variables from above (TKE)) 1174 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 1175 ( gg - bb ) * de_dx(k,j,i+1) + & 1176 ( gg - cc ) * de_dx(k,j+1,i) + & 1177 ( gg - dd ) * de_dx(k,j+1,i+1) & 1178 ) / ( 3.0 * gg ) 1179 1180 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1181 de_dx_int = de_dx_int_l 1182 ELSE 1183 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1184 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1185 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1186 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1187 ) / ( 3.0 * gg ) 1188 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / & 1189 dz * ( de_dx_int_u - de_dx_int_l ) 1190 ENDIF 1191 1192 ! 1193 !-- Interpolate the TKE gradient along y 1194 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1195 ( gg - bb ) * de_dy(k,j,i+1) + & 1196 ( gg - cc ) * de_dy(k,j+1,i) + & 1197 ( gg - dd ) * de_dy(k,j+1,i+1) & 1198 ) / ( 3.0 * gg ) 1199 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1200 de_dy_int = de_dy_int_l 1201 ELSE 1202 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1203 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1204 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1205 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1206 ) / ( 3.0 * gg ) 1207 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / & 1208 dz * ( de_dy_int_u - de_dy_int_l ) 1209 ENDIF 1210 1211 ! 1212 !-- Interpolate the TKE gradient along z 1213 IF ( particles(n)%z < 0.5 * dz ) THEN 1214 de_dz_int = 0.0 1215 ELSE 1216 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1217 ( gg - bb ) * de_dz(k,j,i+1) + & 1218 ( gg - cc ) * de_dz(k,j+1,i) + & 1219 ( gg - dd ) * de_dz(k,j+1,i+1) & 1220 ) / ( 3.0 * gg ) 1221 1222 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1223 de_dz_int = de_dz_int_l 1224 ELSE 1225 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1226 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1227 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1228 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1229 ) / ( 3.0 * gg ) 1230 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /& 1231 dz * ( de_dz_int_u - de_dz_int_l ) 1232 ENDIF 1233 ENDIF 1234 1235 ! 1236 !-- Interpolate the dissipation of TKE 1237 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1238 ( gg - bb ) * diss(k,j,i+1) + & 1239 ( gg - cc ) * diss(k,j+1,i) + & 1240 ( gg - dd ) * diss(k,j+1,i+1) & 1241 ) / ( 3.0 * gg ) 1242 1243 IF ( k+1 == nzt+1 ) THEN 1244 diss_int = diss_int_l 1245 ELSE 1246 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1247 ( gg - bb ) * diss(k+1,j,i+1) + & 1248 ( gg - cc ) * diss(k+1,j+1,i) + & 1249 ( gg - dd ) * diss(k+1,j+1,i+1) & 1250 ) / ( 3.0 * gg ) 1251 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *& 1252 ( diss_int_u - diss_int_l ) 1253 ENDIF 1254 1255 ELSE 1256 1257 ! 1258 !-- If wall between gridpoint 1 and gridpoint 5, then 1259 !-- Neumann boundary condition has to be applied 1260 IF ( gp_outside_of_building(1) == 1 .AND. & 1261 gp_outside_of_building(5) == 0 ) THEN 1262 num_gp = num_gp + 1 1263 location(num_gp,1) = i * dx + 0.5 * dx 1264 location(num_gp,2) = j * dy 1265 location(num_gp,3) = k * dz - 0.5 * dz 1266 ei(num_gp) = e(k,j,i) 1267 dissi(num_gp) = diss(k,j,i) 1268 de_dxi(num_gp) = 0.0 1269 de_dyi(num_gp) = de_dy(k,j,i) 1270 de_dzi(num_gp) = de_dz(k,j,i) 1271 ENDIF 1272 1273 IF ( gp_outside_of_building(5) == 1 .AND. & 1274 gp_outside_of_building(1) == 0 ) THEN 1275 num_gp = num_gp + 1 1276 location(num_gp,1) = i * dx + 0.5 * dx 1277 location(num_gp,2) = j * dy 1278 location(num_gp,3) = k * dz - 0.5 * dz 1279 ei(num_gp) = e(k,j,i+1) 1280 dissi(num_gp) = diss(k,j,i+1) 1281 de_dxi(num_gp) = 0.0 1282 de_dyi(num_gp) = de_dy(k,j,i+1) 1283 de_dzi(num_gp) = de_dz(k,j,i+1) 1284 ENDIF 1285 1286 ! 1287 !-- If wall between gridpoint 5 and gridpoint 6, then 1288 !-- then Neumann boundary condition has to be applied 1289 IF ( gp_outside_of_building(5) == 1 .AND. & 1290 gp_outside_of_building(6) == 0 ) THEN 1291 num_gp = num_gp + 1 1292 location(num_gp,1) = (i+1) * dx 1293 location(num_gp,2) = j * dy + 0.5 * dy 1294 location(num_gp,3) = k * dz - 0.5 * dz 1295 ei(num_gp) = e(k,j,i+1) 1296 dissi(num_gp) = diss(k,j,i+1) 1297 de_dxi(num_gp) = de_dx(k,j,i+1) 1298 de_dyi(num_gp) = 0.0 1299 de_dzi(num_gp) = de_dz(k,j,i+1) 1300 ENDIF 1301 1302 IF ( gp_outside_of_building(6) == 1 .AND. & 1303 gp_outside_of_building(5) == 0 ) THEN 1304 num_gp = num_gp + 1 1305 location(num_gp,1) = (i+1) * dx 1306 location(num_gp,2) = j * dy + 0.5 * dy 1307 location(num_gp,3) = k * dz - 0.5 * dz 1308 ei(num_gp) = e(k,j+1,i+1) 1309 dissi(num_gp) = diss(k,j+1,i+1) 1310 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1311 de_dyi(num_gp) = 0.0 1312 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1313 ENDIF 1314 1315 ! 1316 !-- If wall between gridpoint 2 and gridpoint 6, then 1317 !-- Neumann boundary condition has to be applied 1318 IF ( gp_outside_of_building(2) == 1 .AND. & 1319 gp_outside_of_building(6) == 0 ) THEN 1320 num_gp = num_gp + 1 1321 location(num_gp,1) = i * dx + 0.5 * dx 1322 location(num_gp,2) = (j+1) * dy 1323 location(num_gp,3) = k * dz - 0.5 * dz 1324 ei(num_gp) = e(k,j+1,i) 1325 dissi(num_gp) = diss(k,j+1,i) 1326 de_dxi(num_gp) = 0.0 1327 de_dyi(num_gp) = de_dy(k,j+1,i) 1328 de_dzi(num_gp) = de_dz(k,j+1,i) 1329 ENDIF 1330 1331 IF ( gp_outside_of_building(6) == 1 .AND. & 1332 gp_outside_of_building(2) == 0 ) THEN 1333 num_gp = num_gp + 1 1334 location(num_gp,1) = i * dx + 0.5 * dx 1335 location(num_gp,2) = (j+1) * dy 1336 location(num_gp,3) = k * dz - 0.5 * dz 1337 ei(num_gp) = e(k,j+1,i+1) 1338 dissi(num_gp) = diss(k,j+1,i+1) 1339 de_dxi(num_gp) = 0.0 1340 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1341 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1342 ENDIF 1343 1344 ! 1345 !-- If wall between gridpoint 1 and gridpoint 2, then 1346 !-- Neumann boundary condition has to be applied 1347 IF ( gp_outside_of_building(1) == 1 .AND. & 1348 gp_outside_of_building(2) == 0 ) THEN 1349 num_gp = num_gp + 1 1350 location(num_gp,1) = i * dx 1351 location(num_gp,2) = j * dy + 0.5 * dy 1352 location(num_gp,3) = k * dz - 0.5 * dz 1353 ei(num_gp) = e(k,j,i) 1354 dissi(num_gp) = diss(k,j,i) 1355 de_dxi(num_gp) = de_dx(k,j,i) 1356 de_dyi(num_gp) = 0.0 1357 de_dzi(num_gp) = de_dz(k,j,i) 1358 ENDIF 1359 1360 IF ( gp_outside_of_building(2) == 1 .AND. & 1361 gp_outside_of_building(1) == 0 ) THEN 1362 num_gp = num_gp + 1 1363 location(num_gp,1) = i * dx 1364 location(num_gp,2) = j * dy + 0.5 * dy 1365 location(num_gp,3) = k * dz - 0.5 * dz 1366 ei(num_gp) = e(k,j+1,i) 1367 dissi(num_gp) = diss(k,j+1,i) 1368 de_dxi(num_gp) = de_dx(k,j+1,i) 1369 de_dyi(num_gp) = 0.0 1370 de_dzi(num_gp) = de_dz(k,j+1,i) 1371 ENDIF 1372 1373 ! 1374 !-- If wall between gridpoint 3 and gridpoint 7, then 1375 !-- Neumann boundary condition has to be applied 1376 IF ( gp_outside_of_building(3) == 1 .AND. & 1377 gp_outside_of_building(7) == 0 ) THEN 1378 num_gp = num_gp + 1 1379 location(num_gp,1) = i * dx + 0.5 * dx 1380 location(num_gp,2) = j * dy 1381 location(num_gp,3) = k * dz + 0.5 * dz 1382 ei(num_gp) = e(k+1,j,i) 1383 dissi(num_gp) = diss(k+1,j,i) 1384 de_dxi(num_gp) = 0.0 1385 de_dyi(num_gp) = de_dy(k+1,j,i) 1386 de_dzi(num_gp) = de_dz(k+1,j,i) 1387 ENDIF 1388 1389 IF ( gp_outside_of_building(7) == 1 .AND. & 1390 gp_outside_of_building(3) == 0 ) THEN 1391 num_gp = num_gp + 1 1392 location(num_gp,1) = i * dx + 0.5 * dx 1393 location(num_gp,2) = j * dy 1394 location(num_gp,3) = k * dz + 0.5 * dz 1395 ei(num_gp) = e(k+1,j,i+1) 1396 dissi(num_gp) = diss(k+1,j,i+1) 1397 de_dxi(num_gp) = 0.0 1398 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1399 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1400 ENDIF 1401 1402 ! 1403 !-- If wall between gridpoint 7 and gridpoint 8, then 1404 !-- Neumann boundary condition has to be applied 1405 IF ( gp_outside_of_building(7) == 1 .AND. & 1406 gp_outside_of_building(8) == 0 ) THEN 1407 num_gp = num_gp + 1 1408 location(num_gp,1) = (i+1) * dx 1409 location(num_gp,2) = j * dy + 0.5 * dy 1410 location(num_gp,3) = k * dz + 0.5 * dz 1411 ei(num_gp) = e(k+1,j,i+1) 1412 dissi(num_gp) = diss(k+1,j,i+1) 1413 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1414 de_dyi(num_gp) = 0.0 1415 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1416 ENDIF 1417 1418 IF ( gp_outside_of_building(8) == 1 .AND. & 1419 gp_outside_of_building(7) == 0 ) THEN 1420 num_gp = num_gp + 1 1421 location(num_gp,1) = (i+1) * dx 1422 location(num_gp,2) = j * dy + 0.5 * dy 1423 location(num_gp,3) = k * dz + 0.5 * dz 1424 ei(num_gp) = e(k+1,j+1,i+1) 1425 dissi(num_gp) = diss(k+1,j+1,i+1) 1426 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1427 de_dyi(num_gp) = 0.0 1428 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1429 ENDIF 1430 1431 ! 1432 !-- If wall between gridpoint 4 and gridpoint 8, then 1433 !-- Neumann boundary condition has to be applied 1434 IF ( gp_outside_of_building(4) == 1 .AND. & 1435 gp_outside_of_building(8) == 0 ) THEN 1436 num_gp = num_gp + 1 1437 location(num_gp,1) = i * dx + 0.5 * dx 1438 location(num_gp,2) = (j+1) * dy 1439 location(num_gp,3) = k * dz + 0.5 * dz 1440 ei(num_gp) = e(k+1,j+1,i) 1441 dissi(num_gp) = diss(k+1,j+1,i) 1442 de_dxi(num_gp) = 0.0 1443 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1444 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1445 ENDIF 1446 1447 IF ( gp_outside_of_building(8) == 1 .AND. & 1448 gp_outside_of_building(4) == 0 ) THEN 1449 num_gp = num_gp + 1 1450 location(num_gp,1) = i * dx + 0.5 * dx 1451 location(num_gp,2) = (j+1) * dy 1452 location(num_gp,3) = k * dz + 0.5 * dz 1453 ei(num_gp) = e(k+1,j+1,i+1) 1454 dissi(num_gp) = diss(k+1,j+1,i+1) 1455 de_dxi(num_gp) = 0.0 1456 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1457 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1458 ENDIF 1459 1460 ! 1461 !-- If wall between gridpoint 3 and gridpoint 4, then 1462 !-- Neumann boundary condition has to be applied 1463 IF ( gp_outside_of_building(3) == 1 .AND. & 1464 gp_outside_of_building(4) == 0 ) THEN 1465 num_gp = num_gp + 1 1466 location(num_gp,1) = i * dx 1467 location(num_gp,2) = j * dy + 0.5 * dy 1468 location(num_gp,3) = k * dz + 0.5 * dz 1469 ei(num_gp) = e(k+1,j,i) 1470 dissi(num_gp) = diss(k+1,j,i) 1471 de_dxi(num_gp) = de_dx(k+1,j,i) 1472 de_dyi(num_gp) = 0.0 1473 de_dzi(num_gp) = de_dz(k+1,j,i) 1474 ENDIF 1475 1476 IF ( gp_outside_of_building(4) == 1 .AND. & 1477 gp_outside_of_building(3) == 0 ) THEN 1478 num_gp = num_gp + 1 1479 location(num_gp,1) = i * dx 1480 location(num_gp,2) = j * dy + 0.5 * dy 1481 location(num_gp,3) = k * dz + 0.5 * dz 1482 ei(num_gp) = e(k+1,j+1,i) 1483 dissi(num_gp) = diss(k+1,j+1,i) 1484 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1485 de_dyi(num_gp) = 0.0 1486 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1487 ENDIF 1488 1489 ! 1490 !-- If wall between gridpoint 1 and gridpoint 3, then 1491 !-- Neumann boundary condition has to be applied 1492 !-- (only one case as only building beneath is possible) 1493 IF ( gp_outside_of_building(1) == 1 .AND. & 1494 gp_outside_of_building(3) == 0 ) THEN 1495 num_gp = num_gp + 1 1496 location(num_gp,1) = i * dx 1497 location(num_gp,2) = j * dy 1498 location(num_gp,3) = k * dz 1499 ei(num_gp) = e(k+1,j,i) 1500 dissi(num_gp) = diss(k+1,j,i) 1501 de_dxi(num_gp) = de_dx(k+1,j,i) 1502 de_dyi(num_gp) = de_dy(k+1,j,i) 1503 de_dzi(num_gp) = 0.0 1504 ENDIF 1505 1506 ! 1507 !-- If wall between gridpoint 5 and gridpoint 7, then 1508 !-- Neumann boundary condition has to be applied 1509 !-- (only one case as only building beneath is possible) 1510 IF ( gp_outside_of_building(5) == 1 .AND. & 1511 gp_outside_of_building(7) == 0 ) THEN 1512 num_gp = num_gp + 1 1513 location(num_gp,1) = (i+1) * dx 1514 location(num_gp,2) = j * dy 1515 location(num_gp,3) = k * dz 1516 ei(num_gp) = e(k+1,j,i+1) 1517 dissi(num_gp) = diss(k+1,j,i+1) 1518 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1519 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1520 de_dzi(num_gp) = 0.0 1521 ENDIF 1522 1523 ! 1524 !-- If wall between gridpoint 2 and gridpoint 4, then 1525 !-- Neumann boundary condition has to be applied 1526 !-- (only one case as only building beneath is possible) 1527 IF ( gp_outside_of_building(2) == 1 .AND. & 1528 gp_outside_of_building(4) == 0 ) THEN 1529 num_gp = num_gp + 1 1530 location(num_gp,1) = i * dx 1531 location(num_gp,2) = (j+1) * dy 1532 location(num_gp,3) = k * dz 1533 ei(num_gp) = e(k+1,j+1,i) 1534 dissi(num_gp) = diss(k+1,j+1,i) 1535 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1536 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1537 de_dzi(num_gp) = 0.0 1538 ENDIF 1539 1540 ! 1541 !-- If wall between gridpoint 6 and gridpoint 8, then 1542 !-- Neumann boundary condition has to be applied 1543 !-- (only one case as only building beneath is possible) 1544 IF ( gp_outside_of_building(6) == 1 .AND. & 1545 gp_outside_of_building(8) == 0 ) THEN 1546 num_gp = num_gp + 1 1547 location(num_gp,1) = (i+1) * dx 1548 location(num_gp,2) = (j+1) * dy 1549 location(num_gp,3) = k * dz 1550 ei(num_gp) = e(k+1,j+1,i+1) 1551 dissi(num_gp) = diss(k+1,j+1,i+1) 1552 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1553 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1554 de_dzi(num_gp) = 0.0 1555 ENDIF 1556 1557 ! 1558 !-- Carry out the interpolation 1559 IF ( num_gp == 1 ) THEN 1560 ! 1561 !-- If only one of the gridpoints is situated outside of the 1562 !-- building, it follows that the values at the particle 1563 !-- location are the same as the gridpoint values 1564 e_int = ei(num_gp) 1565 1566 ELSE IF ( num_gp > 1 ) THEN 1567 1568 d_sum = 0.0 1569 ! 1570 !-- Evaluation of the distances between the gridpoints 1571 !-- contributing to the interpolated values, and the particle 1572 !-- location 1573 DO agp = 1, num_gp 1574 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1575 + ( particles(n)%y-location(agp,2) )**2 & 1576 + ( particles(n)%z-location(agp,3) )**2 1577 d_sum = d_sum + d_gp_pl(agp) 1578 ENDDO 1579 1580 ! 1581 !-- Finally the interpolation can be carried out 1582 e_int = 0.0 1583 diss_int = 0.0 1584 de_dx_int = 0.0 1585 de_dy_int = 0.0 1586 de_dz_int = 0.0 1587 DO agp = 1, num_gp 1588 e_int = e_int + ( d_sum - d_gp_pl(agp) ) * & 1589 ei(agp) / ( (num_gp-1) * d_sum ) 1590 diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * & 1591 dissi(agp) / ( (num_gp-1) * d_sum ) 1592 de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * & 1593 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1594 de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * & 1595 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1596 de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * & 1597 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1598 ENDDO 1599 1600 ENDIF 1601 1602 ENDIF 1603 1604 ENDIF 1072 1605 1073 1606 ! … … 1148 1681 1149 1682 ELSE 1683 1684 ! 1685 !-- Restriction of the size of the new timestep: compared to the 1686 !-- previous timestep the increase must not exceed 200% 1687 1688 dt_particle_old = particles(n)%age - particles(n)%age_m 1689 IF ( dt_particle > 2.0 * dt_particle_old ) THEN 1690 dt_particle = 2.0 * dt_particle_old 1691 ENDIF 1692 1150 1693 ! 1151 1694 !-- For old particles the SGS components are correlated with the 1152 1695 !-- values from the previous timestep. Random numbers have also to 1153 1696 !-- be limited (see above). 1697 !-- As negative values for the subgrid TKE are not allowed, the 1698 !-- change of the subgrid TKE with time cannot be smaller than 1699 !-- -e_int/dt_particle. This value is used as a lower boundary 1700 !-- value for the change of TKE 1701 1702 de_dt_min = - e_int / dt_particle 1703 1704 de_dt = ( e_int - particles(n)%e_m ) / dt_particle_old 1705 1706 IF ( de_dt < de_dt_min ) THEN 1707 de_dt = de_dt_min 1708 ENDIF 1709 1154 1710 particles(n)%speed_x_sgs = particles(n)%speed_x_sgs - & 1155 1711 fs_int * c_0 * diss_int * particles(n)%speed_x_sgs * & 1156 1712 dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + & 1157 ( ( 2.0 * sgs_wfu_part * ( e_int - particles(n)%e_m ) /&1158 dt_particle ) * particles(n)%speed_x_sgs /&1713 ( 2.0 * sgs_wfu_part * de_dt * & 1714 particles(n)%speed_x_sgs / & 1159 1715 ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int & 1160 1716 ) * dt_particle / 2.0 + & … … 1166 1722 fs_int * c_0 * diss_int * particles(n)%speed_y_sgs * & 1167 1723 dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + & 1168 ( ( 2.0 * sgs_wfv_part * ( e_int - particles(n)%e_m ) /&1169 dt_particle ) * particles(n)%speed_y_sgs /&1724 ( 2.0 * sgs_wfv_part * de_dt * & 1725 particles(n)%speed_y_sgs / & 1170 1726 ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int & 1171 1727 ) * dt_particle / 2.0 + & … … 1177 1733 fs_int * c_0 * diss_int * particles(n)%speed_z_sgs * & 1178 1734 dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + & 1179 ( ( 2.0 * sgs_wfw_part * ( e_int - particles(n)%e_m ) /&1180 dt_particle ) * particles(n)%speed_z_sgs /&1735 ( 2.0 * sgs_wfw_part * de_dt * & 1736 particles(n)%speed_z_sgs / & 1181 1737 ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int & 1182 1738 ) * dt_particle / 2.0 + & … … 1203 1759 1204 1760 ENDIF 1761 1762 ! 1763 !-- Remember the old age of the particle ( needed to prevent that a 1764 !-- particle crosses several PEs during one timestep and for the 1765 !-- evaluation of the subgrid particle velocity fluctuations ) 1766 particles(n)%age_m = particles(n)%age 1767 1205 1768 1206 1769 ! … … 1270 1833 1271 1834 ENDDO ! advection loop 1835 1836 1837 ! 1838 !-- Particle reflection from walls 1839 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' ) 1840 1841 DO n = 1, number_of_particles 1842 1843 !------------------------------------------------------- 1844 i2 = (particles(n)%x + 0.5*dx) * ddx 1845 j2 = (particles(n)%y + 0.5*dy) * ddy 1846 k2 = particles(n)%z / dz + 1 1847 1848 particles_x = particles(n)%x 1849 particles_y = particles(n)%y 1850 particles_z = particles(n)%z 1851 1852 1853 IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN 1854 1855 old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle 1856 old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle 1857 old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle 1858 i1= (old_positions_x +0.5*dx) * ddx 1859 j1 =(old_positions_y +0.5*dy) * ddy 1860 k1 = old_positions_z / dz 1861 1862 1863 !-- CASE 1 1864 1865 IF ( particles(n)%x > old_positions_x .and. & 1866 particles(n)%y > old_positions_y )THEN 1867 1868 1869 1870 1871 t_index = 1 1872 1873 Do i = i1, i2 1874 xline = i*dx+0.5*dx 1875 1876 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 1877 t_index = t_index + 1 1878 ENDDO 1879 1880 Do j = j1, j2 1881 yline = j*dy+0.5*dy 1882 1883 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 1884 t_index = t_index + 1 1885 ENDDO 1886 1887 IF ( particles(n)%z < old_positions_z ) THEN 1888 Do k = k1, k2 , -1 1889 zline = k*dz 1890 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 1891 t_index = t_index + 1 1892 ENDDO 1893 ENDIF 1894 t_index_number = t_index-1 1895 !-- sorting t 1896 inc = 1 1897 jr = 1 1898 DO WHILE ( inc <= t_index_number ) 1899 inc = 3 * inc + 1 1900 ENDDO 1901 1902 DO WHILE ( inc > 1 ) 1903 inc = inc / 3 1904 DO ir = inc+1, t_index_number 1905 tmp_t = t(ir) 1906 jr = ir 1907 1908 DO WHILE ( t(jr-inc) > tmp_t ) 1909 t(jr) = t(jr-inc) 1910 jr = jr - inc 1911 IF ( jr <= inc ) THEN 1912 EXIT 1913 ENDIF 1914 ENDDO 1915 t(jr) = tmp_t 1916 1917 ENDDO 1918 ENDDO 1919 1920 1921 1922 1923 Do t_index = 1, t_index_number 1924 1925 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 1926 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 1927 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 1928 1929 1930 1931 i3 = (positions_x + 0.5*dx) * ddx 1932 j3 = (positions_y + 0.5*dy) * ddy 1933 k3 = positions_z / dz 1934 1935 i5 = positions_x * ddx 1936 j5 = positions_y * ddy 1937 k5 = positions_z / dz 1938 1939 1940 1941 1942 1943 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 1944 1945 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 1946 k3 == nzb_s_inner(j5,i3) ) THEN 1947 1948 particles(n)%z = 2*positions_z - particles_z 1949 particles(n)%speed_z = - particles(n)%speed_z 1950 1951 IF ( use_sgs_for_particles ) THEN 1952 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1953 ENDIF 1954 GOTO 999 1955 ENDIF 1956 ENDIF 1957 1958 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 1959 1960 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 1961 k3 == nzb_s_inner(j3,i5) ) THEN 1962 1963 particles(n)%z = 2*positions_z - particles_z 1964 particles(n)%speed_z = - particles(n)%speed_z 1965 1966 IF ( use_sgs_for_particles ) THEN 1967 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1968 ENDIF 1969 1970 GOTO 999 1971 ENDIF 1972 ENDIF 1973 1974 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 1975 1976 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 1977 k3 == nzb_s_inner(j3,i3) ) THEN 1978 1979 particles(n)%z = 2*positions_z - particles_z 1980 particles(n)%speed_z = - particles(n)%speed_z 1981 1982 IF ( use_sgs_for_particles ) THEN 1983 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1984 ENDIF 1985 1986 GOTO 999 1987 ENDIF 1988 1989 IF ( positions_y == j3*dy-0.5*dy .and. & 1990 positions_z < nzb_s_inner(j3,i3)*dz) THEN 1991 1992 particles(n)%y = 2*positions_y - particles_y 1993 particles(n)%speed_y = - particles(n)%speed_y 1994 1995 IF ( use_sgs_for_particles ) THEN 1996 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 1997 ENDIF 1998 1999 GOTO 999 2000 ENDIF 2001 2002 2003 IF ( positions_x == i3*dx-0.5*dx.and. & 2004 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2005 2006 particles(n)%x = 2*positions_x - particles_x 2007 particles(n)%speed_x = - particles(n)%speed_x 2008 2009 IF ( use_sgs_for_particles ) THEN 2010 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2011 GOTO 999 2012 ENDIF 2013 2014 ENDIF 2015 2016 ENDIF 2017 ENDDO 2018 ENDIF 2019 2020 !-- CASE 2 2021 2022 IF ( particles(n)%x > old_positions_x .and. & 2023 particles(n)%y < old_positions_y )THEN 2024 2025 2026 t_index = 1 2027 2028 Do i = i1, i2 2029 xline = i*dx+0.5*dx 2030 2031 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 2032 t_index = t_index + 1 2033 ENDDO 2034 2035 Do j = j1, j2,-1 2036 yline = j*dy-0.5*dy 2037 2038 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2039 t_index = t_index + 1 2040 ENDDO 2041 2042 IF ( particles(n)%z < old_positions_z ) THEN 2043 Do k = k1, k2 , -1 2044 zline = k*dz 2045 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2046 t_index = t_index + 1 2047 ENDDO 2048 ENDIF 2049 t_index_number = t_index-1 2050 2051 !-- sorting t 2052 inc = 1 2053 jr = 1 2054 DO WHILE ( inc <= t_index_number ) 2055 inc = 3 * inc + 1 2056 ENDDO 2057 2058 DO WHILE ( inc > 1 ) 2059 inc = inc / 3 2060 DO ir = inc+1, t_index_number 2061 tmp_t = t(ir) 2062 jr = ir 2063 2064 DO WHILE ( t(jr-inc) > tmp_t ) 2065 t(jr) = t(jr-inc) 2066 jr = jr - inc 2067 IF ( jr <= inc ) THEN 2068 EXIT 2069 ENDIF 2070 ENDDO 2071 t(jr) = tmp_t 2072 2073 ENDDO 2074 ENDDO 2075 2076 2077 2078 2079 2080 2081 Do t_index = 1, t_index_number 2082 2083 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2084 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2085 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2086 2087 2088 2089 i3 = (positions_x + 0.5*dx) * ddx 2090 j3 = (positions_y + 0.5*dy) * ddy 2091 k3 = positions_z / dz 2092 2093 i5 = positions_x * ddx 2094 j5 = positions_y * ddy 2095 k5 = positions_z / dz 2096 2097 2098 2099 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2100 2101 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2102 k3 == nzb_s_inner(j3,i5) ) THEN 2103 2104 particles(n)%z = 2*positions_z - particles_z 2105 particles(n)%speed_z = - particles(n)%speed_z 2106 2107 IF ( use_sgs_for_particles ) THEN 2108 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2109 ENDIF 2110 2111 GOTO 999 2112 ENDIF 2113 ENDIF 2114 2115 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2116 2117 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2118 k3 == nzb_s_inner(j3,i3) ) THEN 2119 2120 particles(n)%z = 2*positions_z - particles_z 2121 particles(n)%speed_z = - particles(n)%speed_z 2122 2123 IF ( use_sgs_for_particles ) THEN 2124 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2125 ENDIF 2126 2127 GOTO 999 2128 ENDIF 2129 2130 2131 IF ( positions_x == i3*dx-0.5*dx.and.& 2132 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2133 2134 particles(n)%x = 2*positions_x - particles_x 2135 particles(n)%speed_x = - particles(n)%speed_x 2136 2137 IF ( use_sgs_for_particles ) THEN 2138 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2139 ENDIF 2140 2141 GOTO 999 2142 ENDIF 2143 2144 ENDIF 2145 2146 2147 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2148 2149 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2150 k3 == nzb_s_inner(j5,i3) ) THEN 2151 2152 particles(n)%z = 2*positions_z - particles_z 2153 particles(n)%speed_z = - particles(n)%speed_z 2154 IF ( use_sgs_for_particles ) THEN 2155 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2156 ENDIF 2157 2158 GOTO 999 2159 ENDIF 2160 2161 IF ( positions_y == j5 * dy + 0.5*dy.and.& 2162 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2163 2164 particles(n)%y = 2*positions_y - particles_y 2165 particles(n)%speed_y = - particles(n)%speed_y 2166 2167 IF ( use_sgs_for_particles ) THEN 2168 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2169 ENDIF 2170 GOTO 999 2171 ENDIF 2172 2173 ENDIF 2174 ENDDO 2175 ENDIF 2176 2177 2178 !-- CASE 3 2179 2180 2181 2182 IF ( particles(n)%x < old_positions_x .and. & 2183 particles(n)%y > old_positions_y )THEN 2184 2185 2186 t_index = 1 2187 2188 Do i = i1, i2,-1 2189 xline = i*dx-0.5*dx 2190 2191 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2192 t_index = t_index + 1 2193 ENDDO 2194 2195 Do j = j1, j2 2196 yline = j*dy+0.5*dy 2197 2198 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 2199 t_index = t_index + 1 2200 ENDDO 2201 2202 IF ( particles(n)%z < old_positions_z ) THEN 2203 Do k = k1, k2 , -1 2204 zline = k*dz 2205 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2206 t_index = t_index + 1 2207 ENDDO 2208 ENDIF 2209 2210 t_index_number = t_index-1 2211 2212 !-- sorting t 2213 inc = 1 2214 jr = 1 2215 DO WHILE ( inc <= t_index_number ) 2216 inc = 3 * inc + 1 2217 ENDDO 2218 2219 DO WHILE ( inc > 1 ) 2220 inc = inc / 3 2221 DO ir = inc+1, t_index_number 2222 tmp_t = t(ir) 2223 jr = ir 2224 2225 DO WHILE ( t(jr-inc) > tmp_t ) 2226 t(jr) = t(jr-inc) 2227 jr = jr - inc 2228 IF ( jr <= inc ) THEN 2229 EXIT 2230 ENDIF 2231 ENDDO 2232 t(jr) = tmp_t 2233 2234 ENDDO 2235 ENDDO 2236 2237 2238 2239 Do t_index = 1, t_index_number 2240 2241 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2242 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2243 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2244 2245 2246 2247 i3 = (positions_x + 0.5*dx) * ddx 2248 j3 = (positions_y + 0.5*dy) * ddy 2249 k3 = positions_z / dz 2250 2251 i5 = positions_x * ddx 2252 j5 = positions_y * ddy 2253 k5 = positions_z / dz 2254 2255 2256 2257 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2258 2259 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2260 k3 == nzb_s_inner(j5,i3) ) THEN 2261 2262 particles(n)%z = 2*positions_z - particles_z 2263 particles(n)%speed_z = - particles(n)%speed_z 2264 2265 IF ( use_sgs_for_particles ) THEN 2266 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2267 ENDIF 2268 2269 GOTO 999 2270 ENDIF 2271 ENDIF 2272 2273 2274 2275 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2276 2277 2278 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2279 k3 == nzb_s_inner(j3,i3) ) THEN 2280 2281 particles(n)%z = 2*positions_z - particles_z 2282 particles(n)%speed_z = - particles(n)%speed_z 2283 2284 IF ( use_sgs_for_particles ) THEN 2285 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2286 ENDIF 2287 2288 GOTO 999 2289 ENDIF 2290 2291 IF ( positions_y == j3*dy-0.5*dy .and. & 2292 positions_z < nzb_s_inner(j3,i3)*dz) THEN 2293 2294 particles(n)%y = 2*positions_y - particles_y 2295 particles(n)%speed_y = - particles(n)%speed_y 2296 2297 IF ( use_sgs_for_particles ) THEN 2298 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2299 ENDIF 2300 2301 GOTO 999 2302 ENDIF 2303 2304 ENDIF 2305 2306 2307 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2308 2309 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2310 k3 == nzb_s_inner(j3,i5) ) THEN 2311 2312 particles(n)%z = 2*positions_z - particles_z 2313 particles(n)%speed_z = - particles(n)%speed_z 2314 2315 IF ( use_sgs_for_particles ) THEN 2316 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2317 ENDIF 2318 2319 2320 GOTO 999 2321 ENDIF 2322 2323 IF ( positions_x == i5 * dx + 0.5*dx .and. & 2324 positions_z < nzb_s_inner(j3,i5)*dz ) THEN 2325 2326 particles(n)%x= 2*positions_x - particles_x 2327 particles(n)%speed_x = - particles(n)%speed_x 2328 2329 IF ( use_sgs_for_particles ) THEN 2330 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2331 ENDIF 2332 2333 GOTO 999 2334 ENDIF 2335 2336 ENDIF 2337 ENDDO 2338 ENDIF 2339 2340 2341 !-- CASE 4 2342 2343 2344 IF ( particles(n)%x < old_positions_x .and. & 2345 particles(n)%y < old_positions_y )THEN 2346 2347 2348 t_index = 1 2349 2350 Do i = i1, i2,-1 2351 xline = i*dx-0.5*dx 2352 2353 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2354 t_index = t_index + 1 2355 2356 2357 2358 ENDDO 2359 2360 Do j = j1, j2,-1 2361 yline = j*dy-0.5*dy 2362 2363 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2364 t_index = t_index + 1 2365 2366 2367 ENDDO 2368 2369 IF ( particles(n)%z < old_positions_z ) THEN 2370 Do k = k1, k2 , -1 2371 zline = k*dz 2372 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2373 t_index = t_index + 1 2374 2375 2376 ENDDO 2377 ENDIF 2378 2379 t_index_number = t_index-1 2380 2381 !-- sorting t 2382 inc = 1 2383 jr = 1 2384 DO WHILE ( inc <= t_index_number ) 2385 inc = 3 * inc + 1 2386 ENDDO 2387 2388 DO WHILE ( inc > 1 ) 2389 inc = inc / 3 2390 DO ir = inc+1, t_index_number 2391 tmp_t = t(ir) 2392 jr = ir 2393 2394 DO WHILE ( t(jr-inc) > tmp_t ) 2395 t(jr) = t(jr-inc) 2396 jr = jr - inc 2397 IF ( jr <= inc ) THEN 2398 EXIT 2399 ENDIF 2400 ENDDO 2401 t(jr) = tmp_t 2402 2403 ENDDO 2404 ENDDO 2405 2406 2407 Do t_index = 1, t_index_number 2408 2409 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2410 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2411 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2412 2413 2414 2415 i3 = (positions_x + 0.5*dx) * ddx 2416 j3 = (positions_y + 0.5*dy) * ddy 2417 k3 = positions_z / dz 2418 2419 i5 = positions_x * ddx 2420 j5 = positions_y * ddy 2421 k5 = positions_z / dz 2422 2423 2424 2425 2426 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2427 2428 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2429 k3 == nzb_s_inner(j3,i3) ) THEN 2430 2431 particles(n)%z = 2*positions_z - particles_z 2432 particles(n)%speed_z = - particles(n)%speed_z 2433 2434 IF ( use_sgs_for_particles ) THEN 2435 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2436 ENDIF 2437 2438 2439 GOTO 999 2440 ENDIF 2441 ENDIF 2442 2443 2444 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2445 2446 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2447 k3 == nzb_s_inner(j3,i5) ) THEN 2448 2449 particles(n)%z = 2*positions_z - particles_z 2450 particles(n)%speed_z = - particles(n)%speed_z 2451 2452 IF ( use_sgs_for_particles ) THEN 2453 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2454 ENDIF 2455 2456 GOTO 999 2457 ENDIF 2458 2459 IF ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.& 2460 positions_z < nzb_s_inner(j3,i5)*dz) THEN 2461 2462 particles(n)%x= 2*positions_x - particles_x 2463 particles(n)%speed_x = - particles(n)%speed_x 2464 2465 IF ( use_sgs_for_particles ) THEN 2466 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2467 ENDIF 2468 GOTO 999 2469 ENDIF 2470 2471 ENDIF 2472 2473 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0) THEN 2474 2475 2476 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2477 k5 == nzb_s_inner(j5,i3) ) THEN 2478 2479 particles(n)%z = 2*positions_z - particles_z 2480 particles(n)%speed_z = - particles(n)%speed_z 2481 2482 IF ( use_sgs_for_particles ) THEN 2483 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2484 ENDIF 2485 2486 GOTO 999 2487 ENDIF 2488 2489 IF ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.& 2490 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2491 2492 particles(n)%y= 2*positions_y - particles_y 2493 particles(n)%speed_y = - particles(n)%speed_y 2494 2495 IF ( use_sgs_for_particles ) THEN 2496 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2497 ENDIF 2498 2499 GOTO 999 2500 ENDIF 2501 2502 2503 ENDIF 2504 ENDDO 2505 ENDIF 2506 2507 999 CONTINUE 2508 2509 ENDIF 2510 2511 !--------------------------------------------------- 2512 2513 ENDDO ! particle reflection from walls 2514 2515 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' ) 2516 2517 ! 2518 !-- User-defined actions after the evaluation of the new particle position 2519 CALL user_advec_particles 1272 2520 1273 2521 ! … … 1500 2748 trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1501 2749 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1502 0 , 0, 0, 0 )2750 0.0, 0, 0, 0, 0 ) 1503 2751 trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1504 2752 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1505 0 , 0, 0, 0 )2753 0.0, 0, 0, 0, 0 ) 1506 2754 1507 2755 IF ( use_particle_tails ) THEN … … 1933 3181 trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1934 3182 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1935 0 , 0, 0, 0 )3183 0.0, 0, 0, 0, 0 ) 1936 3184 trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1937 3185 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1938 0 , 0, 0, 0 )3186 0.0, 0, 0, 0, 0 ) 1939 3187 1940 3188 IF ( use_particle_tails ) THEN … … 2460 3708 2461 3709 nn = particles(n)%tail_id 3710 3711 ! 3712 !-- Stop if particles have moved further than the length of one 3713 !-- PE subdomain 3714 IF ( ABS(particles(n)%speed_x) > & 3715 ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. & 3716 ABS(particles(n)%speed_y) > & 3717 ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN 3718 3719 PRINT*, '+++ advec_particles: particle too fast. n = ', n 3720 #if defined( __parallel ) 3721 CALL MPI_ABORT( comm2d, 9999, ierr ) 3722 #else 3723 CALL local_stop 3724 #endif 3725 ENDIF 2462 3726 2463 3727 IF ( particles(n)%age > particle_maximum_age .AND. &
Note: See TracChangeset
for help on using the changeset viewer.
