Changeset 4121 for palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
- Timestamp:
- Jul 26, 2019 10:01:22 AM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4114 r4121 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of an simple method for interpolating the velocities to 28 ! particle position 29 ! 30 ! 4114 2019-07-23 14:09:27Z schwenkel 27 31 ! Bugfix: Added working precision for if statement 28 32 ! … … 72 76 ! with scalar grid point of same index. 73 77 ! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod 74 ! lpm_pack_all_arrays -> lpm_sort_ in_subboxes, lpm_pack_arrays -> lpm_pack78 ! lpm_pack_all_arrays -> lpm_sort_and_delete, lpm_pack_arrays -> lpm_pack 75 79 ! lpm_sort -> lpm_sort_timeloop_done 76 80 ! … … 169 173 ! Description: 170 174 ! ------------ 171 !> 175 !> 172 176 !------------------------------------------------------------------------------! 173 177 MODULE lagrangian_particle_model_mod … … 178 182 ONLY: de_dx, de_dy, de_dz, dzw, zu, zw, ql_c, ql_v, ql_vp, hyp, & 179 183 pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p, & 180 d_exner 184 d_exner, u_p, v_p, w_p 181 185 182 186 USE averaging, & … … 263 267 IMPLICIT NONE 264 268 265 CHARACTER(LEN=15) :: aero_species = 'nacl' !< aerosol species 266 CHARACTER(LEN=15) :: aero_type = 'maritime' !< aerosol type 267 CHARACTER(LEN=15) :: bc_par_lr = 'cyclic' !< left/right boundary condition 268 CHARACTER(LEN=15) :: bc_par_ns = 'cyclic' !< north/south boundary condition 269 CHARACTER(LEN=15) :: bc_par_b = 'reflect' !< bottom boundary condition 270 CHARACTER(LEN=15) :: bc_par_t = 'absorb' !< top boundary condition 271 CHARACTER(LEN=15) :: collision_kernel = 'none' !< collision kernel 272 273 CHARACTER(LEN=5) :: splitting_function = 'gamma' !< function for calculation critical weighting factor 274 CHARACTER(LEN=5) :: splitting_mode = 'const' !< splitting mode 269 CHARACTER(LEN=15) :: aero_species = 'nacl' !< aerosol species 270 CHARACTER(LEN=15) :: aero_type = 'maritime' !< aerosol type 271 CHARACTER(LEN=15) :: bc_par_lr = 'cyclic' !< left/right boundary condition 272 CHARACTER(LEN=15) :: bc_par_ns = 'cyclic' !< north/south boundary condition 273 CHARACTER(LEN=15) :: bc_par_b = 'reflect' !< bottom boundary condition 274 CHARACTER(LEN=15) :: bc_par_t = 'absorb' !< top boundary condition 275 CHARACTER(LEN=15) :: collision_kernel = 'none' !< collision kernel 276 277 CHARACTER(LEN=5) :: splitting_function = 'gamma' !< function for calculation critical weighting factor 278 CHARACTER(LEN=5) :: splitting_mode = 'const' !< splitting mode 279 280 CHARACTER(LEN=25) :: particle_interpolation = 'trilinear' !< interpolation method for calculatin the particle 275 281 276 282 INTEGER(iwp) :: deleted_particles = 0 !< number of deleted particles per time step … … 496 502 END INTERFACE dealloc_particles_array 497 503 498 INTERFACE lpm_sort_ in_subboxes499 MODULE PROCEDURE lpm_sort_ in_subboxes500 END INTERFACE lpm_sort_ in_subboxes504 INTERFACE lpm_sort_and_delete 505 MODULE PROCEDURE lpm_sort_and_delete 506 END INTERFACE lpm_sort_and_delete 501 507 502 508 INTERFACE lpm_sort_timeloop_done … … 549 555 particles_per_point, & 550 556 particle_advection_start, & 557 particle_interpolation, & 551 558 particle_maximum_age, & 552 559 pdx, & … … 608 615 particles_per_point, & 609 616 particle_advection_start, & 617 particle_interpolation, & 610 618 particle_maximum_age, & 611 619 pdx, & … … 849 857 IF ( collision_kernel(6:9) == 'fast' ) use_kernel_tables = .TRUE. 850 858 859 ! 860 !-- Subgrid scale velocites with the simple interpolation method for resolved 861 !-- velocites is not implemented for passive particles. However, for cloud 862 !-- it can be combined as the sgs-velocites for active particles are 863 !-- calculated differently, i.e. no subboxes are needed. 864 IF ( .NOT. TRIM(particle_interpolation) == 'trilinear' .AND. & 865 use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 866 message_string = 'subrgrid scale velocities in combination with ' // & 867 'simple interpolation method is not ' // & 868 'implemented' 869 CALL message( 'check_parameters', 'PA0659', 1, 2, 0, 6, 0 ) 870 ENDIF 871 851 872 END SUBROUTINE 852 873 … … 1067 1088 CASE ( 'reflect' ) 1068 1089 ibc_par_t = 2 1069 1090 1070 1091 CASE ( 'nested' ) 1071 1092 ibc_par_t = 3 … … 1087 1108 CASE ( 'reflect' ) 1088 1109 ibc_par_lr = 2 1089 1110 1090 1111 CASE ( 'nested' ) 1091 1112 ibc_par_lr = 3 … … 1107 1128 CASE ( 'reflect' ) 1108 1129 ibc_par_ns = 2 1109 1130 1110 1131 CASE ( 'nested' ) 1111 1132 ibc_par_ns = 3 … … 1546 1567 !-- which is needed for a fast interpolation of the LES fields on the particle 1547 1568 !-- position. 1548 CALL lpm_sort_ in_subboxes1569 CALL lpm_sort_and_delete 1549 1570 ! 1550 1571 !-- Determine the current number of particles … … 1623 1644 !-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions. 1624 1645 IF ( TRIM(aero_type) .EQ. 'polar' ) THEN 1625 na = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6 1626 rm = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6 1646 na = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp 1647 rm = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp 1627 1648 log_sigma = (/ 0.245, 0.300, 0.291 /) 1628 1649 ELSEIF ( TRIM(aero_type) .EQ. 'background' ) THEN 1629 na = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6 1630 rm = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6 1650 na = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp 1651 rm = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp 1631 1652 log_sigma = (/ 0.645, 0.253, 0.425 /) 1632 1653 ELSEIF ( TRIM(aero_type) .EQ. 'maritime' ) THEN 1633 na = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6 1634 rm = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6 1654 na = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp 1655 rm = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp 1635 1656 log_sigma = (/ 0.657, 0.210, 0.396 /) 1636 1657 ELSEIF ( TRIM(aero_type) .EQ. 'continental' ) THEN 1637 na = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6 1638 rm = (/ 0.01, 0.058, 0.9 /) * 1.0E-6 1658 na = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp 1659 rm = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp 1639 1660 log_sigma = (/ 0.161, 0.217, 0.380 /) 1640 1661 ELSEIF ( TRIM(aero_type) .EQ. 'desert' ) THEN 1641 na = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6 1642 rm = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6 1662 na = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp 1663 rm = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp 1643 1664 log_sigma = (/ 0.247, 0.770, 0.438 /) 1644 1665 ELSEIF ( TRIM(aero_type) .EQ. 'rural' ) THEN 1645 na = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6 1646 rm = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6 1666 na = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp 1667 rm = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp 1647 1668 log_sigma = (/ 0.225, 0.557, 0.266 /) 1648 1669 ELSEIF ( TRIM(aero_type) .EQ. 'urban' ) THEN 1649 na = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6 1650 rm = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6 1670 na = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp 1671 rm = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp 1651 1672 log_sigma = (/ 0.245, 0.666, 0.337 /) 1652 1673 ELSEIF ( TRIM(aero_type) .EQ. 'user' ) THEN … … 1679 1700 particles(n)%aux1 = r_mid 1680 1701 particles(n)%weight_factor = & 1681 ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) * &1682 EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) + &1683 na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) * &1684 EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) + &1685 na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) * &1686 EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) ) &1702 ( na(1) / ( SQRT( 2.0_wp * pi ) * log_sigma(1) ) * & 1703 EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0_wp * log_sigma(1)**2 ) ) + & 1704 na(2) / ( SQRT( 2.0_wp * pi ) * log_sigma(2) ) * & 1705 EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0_wp * log_sigma(2)**2 ) ) + & 1706 na(3) / ( SQRT( 2.0_wp * pi ) * log_sigma(3) ) * & 1707 EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0_wp * log_sigma(3)**2 ) ) & 1687 1708 ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dzw(kp) ) 1688 1709 … … 1700 1721 ! 1701 1722 !-- Unnecessary particles will be deleted 1702 IF ( particles(n)%weight_factor .LE. 0.0 ) particles(n)%particle_mask = .FALSE.1723 IF ( particles(n)%weight_factor .LE. 0.0_wp ) particles(n)%particle_mask = .FALSE. 1703 1724 1704 1725 ENDDO … … 1741 1762 1742 1763 END SUBROUTINE lpm_init_aerosols 1743 1744 1764 1765 1745 1766 !------------------------------------------------------------------------------! 1746 1767 ! Description: … … 1752 1773 !------------------------------------------------------------------------------! 1753 1774 SUBROUTINE lpm_init_sgs_tke 1754 1775 1755 1776 USE statistics, & 1756 1777 ONLY: flow_statistics_called, hom, sums, sums_l … … 2013 2034 LOGICAL :: first_loop_stride !< 2014 2035 2015 2036 2016 2037 SELECT CASE ( location ) 2017 2038 … … 2098 2119 CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' ) 2099 2120 CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' ) 2100 2121 2101 2122 ! 2102 2123 !-- If particle advection includes SGS velocity components, calculate the … … 2116 2137 IF ( .NOT. first_loop_stride .AND. use_sgs_for_particles ) & 2117 2138 CALL lpm_sort_timeloop_done 2118 2119 2139 DO i = nxl, nxr 2120 2140 DO j = nys, nyn … … 2164 2184 ! 2165 2185 !-- Particle advection 2166 CALL lpm_advec( i,j,k)2186 CALL lpm_advec(TRIM(particle_interpolation),i,j,k) 2167 2187 ! 2168 2188 !-- Particle reflection from walls. Only applied if the particles … … 2203 2223 ENDDO 2204 2224 ENDDO 2205 2206 2225 steps = steps + 1 2207 2226 dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done) … … 2216 2235 dt_3d_reached = dt_3d_reached_l 2217 2236 #endif 2218 2219 2237 CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' ) 2220 2238 … … 2237 2255 2238 2256 ! 2239 !-- IF .FALSE., lpm_sort_ in_subboxesis done inside pcmp2257 !-- IF .FALSE., lpm_sort_and_delete is done inside pcmp 2240 2258 IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run ) THEN 2241 2259 ! 2242 2260 !-- Pack particles (eliminate those marked for deletion), 2243 2261 !-- determine new number of particles 2244 CALL lpm_sort_ in_subboxes2245 ! 2262 CALL lpm_sort_and_delete 2263 2246 2264 !-- Initialize variables for the next (sub-) timestep, i.e., for marking 2247 2265 !-- those particles to be deleted after the timestep … … 2253 2271 first_loop_stride = .FALSE. 2254 2272 ENDDO ! timestep loop 2255 ! 2273 ! 2256 2274 !-- in case of nested runs do the transfer of particles after every full model time step 2257 2275 IF ( nested_run ) THEN … … 2260 2278 CALL pmcp_p_delete_particles_in_fine_grid_area 2261 2279 2262 CALL lpm_sort_ in_subboxes2280 CALL lpm_sort_and_delete 2263 2281 2264 2282 deleted_particles = 0 … … 2284 2302 2285 2303 CALL cpu_log( log_point(25), 'lpm', 'stop' ) 2286 2304 2287 2305 ! ! 2288 2306 ! !-- Output of particle time series … … 2295 2313 ! ENDIF 2296 2314 ! ENDIF 2297 2315 2298 2316 CASE DEFAULT 2299 2317 CONTINUE … … 2940 2958 !-- Must be called to sort particles into blocks, which is needed for a fast 2941 2959 !-- interpolation of the LES fields on the particle position. 2942 CALL lpm_sort_ in_subboxes2943 2960 CALL lpm_sort_and_delete 2961 2944 2962 2945 2963 END SUBROUTINE lpm_rrd_local_particles … … 2949 2967 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 2950 2968 nysc, nys_on_file, tmp_3d, found ) 2951 2952 2969 2970 2953 2971 USE control_parameters, & 2954 ONLY: length, restart_string 2972 ONLY: length, restart_string 2955 2973 2956 2974 INTEGER(iwp) :: k !< … … 2972 2990 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 2973 2991 2974 2992 2975 2993 found = .TRUE. 2976 2994 2977 2995 SELECT CASE ( restart_string(1:length) ) 2978 2996 2979 2997 CASE ( 'iran' ) ! matching random numbers is still unresolved issue 2980 2998 IF ( k == 1 ) READ ( 13 ) iran, iran_part 2981 2999 2982 3000 CASE ( 'pc_av' ) 2983 3001 IF ( .NOT. ALLOCATED( pc_av ) ) THEN … … 3154 3172 !> grid scale) as well as the boundary conditions of particles. As a next step 3155 3173 !> this submodule should be excluded as an own file. 3156 !------------------------------------------------------------------------------! 3157 SUBROUTINE lpm_advec (ip,jp,kp) 3158 3174 !------------------------------------------------------------------------------! 3175 SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp) 3176 3177 CHARACTER (LEN=*), INTENT(IN) :: interpolation_method !< 3159 3178 LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall 3160 3179 3161 3180 INTEGER(iwp) :: i !< index variable along x 3162 INTEGER(iwp) :: ip !< index variable along x 3163 INTEGER(iwp) :: j !< index variable along y 3164 INTEGER(iwp) :: jp !< index variable along y 3165 INTEGER(iwp) :: k !< index variable along z 3181 INTEGER(iwp) :: i_next !< index variable along x 3182 INTEGER(iwp) :: ip !< index variable along x 3183 INTEGER(iwp) :: iteration_steps = 1 !< amount of iterations steps for corrector step 3184 INTEGER(iwp) :: j !< index variable along y 3185 INTEGER(iwp) :: j_next !< index variable along y 3186 INTEGER(iwp) :: jp !< index variable along y 3187 INTEGER(iwp) :: k !< index variable along z 3166 3188 INTEGER(iwp) :: k_wall !< vertical index of topography top 3167 INTEGER(iwp) :: kp !< index variable along z 3189 INTEGER(iwp) :: kp !< index variable along z 3190 INTEGER(iwp) :: k_next !< index variable along z 3168 3191 INTEGER(iwp) :: kw !< index variable along z 3192 INTEGER(iwp) :: kkw !< index variable along z 3169 3193 INTEGER(iwp) :: n !< loop variable over all particles in a grid box 3170 3194 INTEGER(iwp) :: nb !< block number particles are sorted in 3171 INTEGER(iwp) :: surf_start !< Index on surface data-type for current grid box 3195 INTEGER(iwp) :: particle_end !< end index for partilce loop 3196 INTEGER(iwp) :: particle_start !< start index for particle loop 3197 INTEGER(iwp) :: surf_start !< Index on surface data-type for current grid box 3198 INTEGER(iwp) :: subbox_end !< end index for loop over subboxes in particle advection 3199 INTEGER(iwp) :: subbox_start !< start index for loop over subboxes in particle advection 3200 INTEGER(iwp) :: nn !< loop variable over iterations steps 3172 3201 3173 3202 INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block 3174 3203 INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block 3175 3204 3176 REAL(wp) :: aa !< dummy argument for horizontal particle interpolation 3177 REAL(wp) :: bb !< dummy argument for horizontal particle interpolation 3178 REAL(wp) :: cc !< dummy argument for horizontal particle interpolation 3205 REAL(wp) :: aa !< dummy argument for horizontal particle interpolation 3206 REAL(wp) :: alpha !< interpolation facor for x-direction 3207 3208 REAL(wp) :: bb !< dummy argument for horizontal particle interpolation 3209 REAL(wp) :: beta !< interpolation facor for y-direction 3210 REAL(wp) :: cc !< dummy argument for horizontal particle interpolation 3179 3211 REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation 3180 3212 REAL(wp) :: dd !< dummy argument for horizontal particle interpolation … … 3197 3229 REAL(wp) :: exp_arg !< argument in the exponent - particle radius 3198 3230 REAL(wp) :: exp_term !< exponent term 3231 REAL(wp) :: gamma !< interpolation facor for z-direction 3199 3232 REAL(wp) :: gg !< dummy argument for horizontal particle interpolation 3200 3233 REAL(wp) :: height_p !< dummy argument for logarithmic interpolation … … 3208 3241 REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level 3209 3242 REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level 3243 REAL(wp) :: unext !< calculated particle u-velocity of corrector step 3210 3244 REAL(wp) :: us_int !< friction velocity at particle grid box 3211 3245 REAL(wp) :: usws_int !< surface momentum flux (u component) at particle grid box … … 3213 3247 REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level 3214 3248 REAL(wp) :: vsws_int !< surface momentum flux (u component) at particle grid box 3249 REAL(wp) :: vnext !< calculated particle v-velocity of corrector step 3215 3250 REAL(wp) :: vv_int !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection 3216 3251 REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level 3217 3252 REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level 3253 REAL(wp) :: wnext !< calculated particle w-velocity of corrector step 3218 3254 REAL(wp) :: w_s !< terminal velocity of droplets 3219 3255 REAL(wp) :: x !< dummy argument for horizontal particle interpolation 3220 REAL(wp) :: y !< dummy argument for horizontal particle interpolation 3256 REAL(wp) :: xp !< calculated particle position in x of predictor step 3257 REAL(wp) :: y !< dummy argument for horizontal particle interpolation 3258 REAL(wp) :: yp !< calculated particle position in y of predictor step 3221 3259 REAL(wp) :: z_p !< surface layer height (0.5 dz) 3260 REAL(wp) :: zp !< calculated particle position in z of predictor step 3222 3261 3223 3262 REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity … … 3249 3288 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 3250 3289 3251 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers 3290 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers 3252 3291 3253 3292 CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' ) 3254 3255 3293 ! 3256 3294 !-- Determine height of Prandtl layer and distance between Prandtl-layer 3257 !-- height and horizontal mean roughness height, which are required for 3258 !-- vertical logarithmic interpolation of horizontal particle speeds 3295 !-- height and horizontal mean roughness height, which are required for 3296 !-- vertical logarithmic interpolation of horizontal particle speeds 3259 3297 !-- (for particles below first vertical grid level). 3260 3298 z_p = zu(nzb+1) - zw(nzb) 3261 3299 d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global ) 3262 3300 3263 start_index = grid_particles(kp,jp,ip)%start_index3264 end_index = grid_particles(kp,jp,ip)%end_index3265 3266 3301 xv = particles(1:number_of_particles)%x 3267 3302 yv = particles(1:number_of_particles)%y 3268 3303 zv = particles(1:number_of_particles)%z 3269 3270 DO nb = 0, 7 3271 ! 3272 !-- Interpolate u velocity-component 3273 i = ip 3274 j = jp + block_offset(nb)%j_off 3275 k = kp + block_offset(nb)%k_off 3276 3277 DO n = start_index(nb), end_index(nb) 3278 ! 3279 !-- Interpolation of the u velocity component onto particle position. 3280 !-- Particles are interpolation bi-linearly in the horizontal and a 3281 !-- linearly in the vertical. An exception is made for particles below 3282 !-- the first vertical grid level in case of a prandtl layer. In this 3283 !-- case the horizontal particle velocity components are determined using 3284 !-- Monin-Obukhov relations (if branch). 3285 !-- First, check if particle is located below first vertical grid level 3286 !-- above topography (Prandtl-layer height) 3287 !-- Determine vertical index of topography top 3288 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 3289 3290 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3291 ! 3292 !-- Resolved-scale horizontal particle velocity is zero below z0. 3293 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3294 u_int(n) = 0.0_wp 3295 ELSE 3296 ! 3297 !-- Determine the sublayer. Further used as index. 3298 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3299 * REAL( number_of_sublayers, KIND=wp ) & 3300 * d_z_p_z0 3301 ! 3302 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3303 !-- interpolate linearly between precalculated logarithm. 3304 log_z_z0_int = log_z_z0(INT(height_p)) & 3305 + ( height_p - INT(height_p) ) & 3306 * ( log_z_z0(INT(height_p)+1) & 3307 - log_z_z0(INT(height_p)) & 3308 ) 3309 ! 3310 !-- Get friction velocity and momentum flux from new surface data 3311 !-- types. 3312 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3313 surf_def_h(0)%end_index(jp,ip) ) THEN 3314 surf_start = surf_def_h(0)%start_index(jp,ip) 3315 !-- Limit friction velocity. In narrow canyons or holes the 3316 !-- friction velocity can become very small, resulting in a too 3317 !-- large particle speed. 3318 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3319 usws_int = surf_def_h(0)%usws(surf_start) 3320 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3321 surf_lsm_h%end_index(jp,ip) ) THEN 3322 surf_start = surf_lsm_h%start_index(jp,ip) 3323 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3324 usws_int = surf_lsm_h%usws(surf_start) 3325 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3326 surf_usm_h%end_index(jp,ip) ) THEN 3327 surf_start = surf_usm_h%start_index(jp,ip) 3328 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3329 usws_int = surf_usm_h%usws(surf_start) 3304 dt_particle = dt_3d 3305 3306 3307 SELECT CASE ( interpolation_method ) 3308 3309 ! 3310 !-- This case uses a simple interpolation method for the particle velocites, 3311 !-- and applying a predictor-corrector method. @attention: for the corrector 3312 !-- step the velocities of t(n+1) are required. However, at this moment of 3313 !-- the time integration they are not free of divergence. This interpolation 3314 !-- method is described in more detail in Grabowski et al., 2018 (GMD). 3315 CASE ( 'simple_corrector' ) 3316 ! 3317 !-- Predictor step 3318 kkw = kp - 1 3319 DO n = 1, number_of_particles 3320 3321 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3322 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3323 3324 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3325 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3326 3327 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3328 ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp ) 3329 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3330 3331 ENDDO 3332 3333 ! 3334 !-- Corrector step 3335 DO n = 1, number_of_particles 3336 3337 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3338 3339 DO nn = 1, iteration_steps 3340 3341 ! 3342 !-- Guess new position 3343 xp = particles(n)%x + u_int(n) * dt_particle(n) 3344 yp = particles(n)%y + v_int(n) * dt_particle(n) 3345 zp = particles(n)%z + w_int(n) * dt_particle(n) 3346 ! 3347 !-- x direction 3348 i_next = FLOOR( xp * ddx , KIND=iwp) 3349 alpha = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3350 ! 3351 !-- y direction 3352 j_next = FLOOR( yp * ddy ) 3353 beta = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3354 ! 3355 !-- z_direction 3356 k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0) 3357 gamma = MAX( MIN( ( zp - zw(k_next) ) / & 3358 ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp ) 3359 ! 3360 !-- Calculate part of the corrector step 3361 unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) + & 3362 u_p(k_next+1, j_next, i_next+1) * alpha 3363 3364 vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta ) + & 3365 v_p(k_next+1, j_next+1, i_next ) * beta 3366 3367 wnext = w_p(k_next, j_next, i_next) * ( 1.0_wp - gamma ) + & 3368 w_p(k_next+1, j_next, i_next ) * gamma 3369 3370 ! 3371 !-- Calculate interpolated particle velocity with predictor 3372 !-- corrector step. u_int, v_int and w_int describes the part of 3373 !-- the predictor step. unext, vnext and wnext is the part of the 3374 !-- corrector step. The resulting new position is set below. The 3375 !-- implementation is based on Grabowski et al., 2018 (GMD). 3376 u_int(n) = 0.5_wp * ( u_int(n) + unext ) 3377 v_int(n) = 0.5_wp * ( v_int(n) + vnext ) 3378 w_int(n) = 0.5_wp * ( w_int(n) + wnext ) 3379 3380 ENDDO 3381 ENDDO 3382 3383 3384 ! 3385 !-- This case uses a simple interpolation method for the particle velocites, 3386 !-- and applying a predictor. 3387 CASE ( 'simple_predictor' ) 3388 ! 3389 !-- The particle position for the w velociy is based on the value of kp and kp-1 3390 kkw = kp - 1 3391 DO n = 1, number_of_particles 3392 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3393 3394 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3395 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3396 3397 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3398 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3399 3400 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3401 ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp ) 3402 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3403 ENDDO 3404 ! 3405 !-- The trilinear interpolation. 3406 CASE ( 'trilinear' ) 3407 3408 start_index = grid_particles(kp,jp,ip)%start_index 3409 end_index = grid_particles(kp,jp,ip)%end_index 3410 3411 DO nb = 0, 7 3412 ! 3413 !-- Interpolate u velocity-component 3414 i = ip 3415 j = jp + block_offset(nb)%j_off 3416 k = kp + block_offset(nb)%k_off 3417 3418 DO n = start_index(nb), end_index(nb) 3419 ! 3420 !-- Interpolation of the u velocity component onto particle position. 3421 !-- Particles are interpolation bi-linearly in the horizontal and a 3422 !-- linearly in the vertical. An exception is made for particles below 3423 !-- the first vertical grid level in case of a prandtl layer. In this 3424 !-- case the horizontal particle velocity components are determined using 3425 !-- Monin-Obukhov relations (if branch). 3426 !-- First, check if particle is located below first vertical grid level 3427 !-- above topography (Prandtl-layer height) 3428 !-- Determine vertical index of topography top 3429 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 3430 3431 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3432 ! 3433 !-- Resolved-scale horizontal particle velocity is zero below z0. 3434 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3435 u_int(n) = 0.0_wp 3436 ELSE 3437 ! 3438 !-- Determine the sublayer. Further used as index. 3439 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3440 * REAL( number_of_sublayers, KIND=wp ) & 3441 * d_z_p_z0 3442 ! 3443 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3444 !-- interpolate linearly between precalculated logarithm. 3445 log_z_z0_int = log_z_z0(INT(height_p)) & 3446 + ( height_p - INT(height_p) ) & 3447 * ( log_z_z0(INT(height_p)+1) & 3448 - log_z_z0(INT(height_p)) & 3449 ) 3450 ! 3451 !-- Get friction velocity and momentum flux from new surface data 3452 !-- types. 3453 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3454 surf_def_h(0)%end_index(jp,ip) ) THEN 3455 surf_start = surf_def_h(0)%start_index(jp,ip) 3456 !-- Limit friction velocity. In narrow canyons or holes the 3457 !-- friction velocity can become very small, resulting in a too 3458 !-- large particle speed. 3459 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3460 usws_int = surf_def_h(0)%usws(surf_start) 3461 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3462 surf_lsm_h%end_index(jp,ip) ) THEN 3463 surf_start = surf_lsm_h%start_index(jp,ip) 3464 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3465 usws_int = surf_lsm_h%usws(surf_start) 3466 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3467 surf_usm_h%end_index(jp,ip) ) THEN 3468 surf_start = surf_usm_h%start_index(jp,ip) 3469 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3470 usws_int = surf_usm_h%usws(surf_start) 3471 ENDIF 3472 ! 3473 !-- Neutral solution is applied for all situations, e.g. also for 3474 !-- unstable and stable situations. Even though this is not exact 3475 !-- this saves a lot of CPU time since several calls of intrinsic 3476 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3477 !-- as sensitivity studies revealed no significant effect of 3478 !-- using the neutral solution also for un/stable situations. 3479 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 3480 * log_z_z0_int - u_gtrans 3481 ENDIF 3482 ! 3483 !-- Particle above the first grid level. Bi-linear interpolation in the 3484 !-- horizontal and linear interpolation in the vertical direction. 3485 ELSE 3486 x = xv(n) - i * dx 3487 y = yv(n) + ( 0.5_wp - j ) * dy 3488 aa = x**2 + y**2 3489 bb = ( dx - x )**2 + y**2 3490 cc = x**2 + ( dy - y )**2 3491 dd = ( dx - x )**2 + ( dy - y )**2 3492 gg = aa + bb + cc + dd 3493 3494 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 3495 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * & 3496 u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3497 3498 IF ( k == nzt ) THEN 3499 u_int(n) = u_int_l 3500 ELSE 3501 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 3502 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * & 3503 u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3504 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3505 ( u_int_u - u_int_l ) 3506 ENDIF 3330 3507 ENDIF 3331 3332 ! 3333 !-- Neutral solution is applied for all situations, e.g. also for 3334 !-- unstable and stable situations. Even though this is not exact 3335 !-- this saves a lot of CPU time since several calls of intrinsic 3336 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3337 !-- as sensitivity studies revealed no significant effect of 3338 !-- using the neutral solution also for un/stable situations. 3339 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 3340 * log_z_z0_int - u_gtrans 3341 3342 ENDIF 3343 ! 3344 !-- Particle above the first grid level. Bi-linear interpolation in the 3345 !-- horizontal and linear interpolation in the vertical direction. 3346 ELSE 3347 3348 x = xv(n) - i * dx 3349 y = yv(n) + ( 0.5_wp - j ) * dy 3350 aa = x**2 + y**2 3351 bb = ( dx - x )**2 + y**2 3352 cc = x**2 + ( dy - y )**2 3353 dd = ( dx - x )**2 + ( dy - y )**2 3354 gg = aa + bb + cc + dd 3355 3356 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 3357 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * & 3358 u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3359 3360 IF ( k == nzt ) THEN 3361 u_int(n) = u_int_l 3362 ELSE 3363 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 3364 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * & 3365 u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3366 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3367 ( u_int_u - u_int_l ) 3368 ENDIF 3369 3370 ENDIF 3371 3372 ENDDO 3373 ! 3374 !-- Same procedure for interpolation of the v velocity-component 3375 i = ip + block_offset(nb)%i_off 3376 j = jp 3377 k = kp + block_offset(nb)%k_off 3378 3379 DO n = start_index(nb), end_index(nb) 3380 3381 ! 3382 !-- Determine vertical index of topography top 3383 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 3384 3385 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3386 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3387 ! 3388 !-- Resolved-scale horizontal particle velocity is zero below z0. 3389 v_int(n) = 0.0_wp 3390 ELSE 3391 ! 3392 !-- Determine the sublayer. Further used as index. Please note, 3393 !-- logarithmus can not be reused from above, as in in case of 3394 !-- topography particle on u-grid can be above surface-layer height, 3395 !-- whereas it can be below on v-grid. 3396 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3397 * REAL( number_of_sublayers, KIND=wp ) & 3398 * d_z_p_z0 3399 ! 3400 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3401 !-- interpolate linearly between precalculated logarithm. 3402 log_z_z0_int = log_z_z0(INT(height_p)) & 3403 + ( height_p - INT(height_p) ) & 3404 * ( log_z_z0(INT(height_p)+1) & 3405 - log_z_z0(INT(height_p)) & 3406 ) 3407 ! 3408 !-- Get friction velocity and momentum flux from new surface data 3409 !-- types. 3410 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3411 surf_def_h(0)%end_index(jp,ip) ) THEN 3412 surf_start = surf_def_h(0)%start_index(jp,ip) 3413 !-- Limit friction velocity. In narrow canyons or holes the 3414 !-- friction velocity can become very small, resulting in a too 3415 !-- large particle speed. 3416 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3417 vsws_int = surf_def_h(0)%vsws(surf_start) 3418 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3419 surf_lsm_h%end_index(jp,ip) ) THEN 3420 surf_start = surf_lsm_h%start_index(jp,ip) 3421 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3422 vsws_int = surf_lsm_h%vsws(surf_start) 3423 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3424 surf_usm_h%end_index(jp,ip) ) THEN 3425 surf_start = surf_usm_h%start_index(jp,ip) 3426 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3427 vsws_int = surf_usm_h%vsws(surf_start) 3428 ENDIF 3429 ! 3430 !-- Neutral solution is applied for all situations, e.g. also for 3431 !-- unstable and stable situations. Even though this is not exact 3432 !-- this saves a lot of CPU time since several calls of intrinsic 3433 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3434 !-- as sensitivity studies revealed no significant effect of 3435 !-- using the neutral solution also for un/stable situations. 3436 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 3437 * log_z_z0_int - v_gtrans 3438 3439 ENDIF 3440 3441 ELSE 3442 x = xv(n) + ( 0.5_wp - i ) * dx 3443 y = yv(n) - j * dy 3444 aa = x**2 + y**2 3445 bb = ( dx - x )**2 + y**2 3446 cc = x**2 + ( dy - y )**2 3447 dd = ( dx - x )**2 + ( dy - y )**2 3448 gg = aa + bb + cc + dd 3449 3450 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 3451 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 3452 ) / ( 3.0_wp * gg ) - v_gtrans 3453 3454 IF ( k == nzt ) THEN 3455 v_int(n) = v_int_l 3456 ELSE 3457 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 3458 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 3459 ) / ( 3.0_wp * gg ) - v_gtrans 3460 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3461 ( v_int_u - v_int_l ) 3462 ENDIF 3463 3464 ENDIF 3465 3466 ENDDO 3467 ! 3468 !-- Same procedure for interpolation of the w velocity-component 3469 i = ip + block_offset(nb)%i_off 3470 j = jp + block_offset(nb)%j_off 3471 k = kp - 1 3472 3473 DO n = start_index(nb), end_index(nb) 3474 3475 IF ( vertical_particle_advection(particles(n)%group) ) THEN 3476 3477 x = xv(n) + ( 0.5_wp - i ) * dx 3478 y = yv(n) + ( 0.5_wp - j ) * dy 3479 aa = x**2 + y**2 3480 bb = ( dx - x )**2 + y**2 3481 cc = x**2 + ( dy - y )**2 3482 dd = ( dx - x )**2 + ( dy - y )**2 3483 gg = aa + bb + cc + dd 3484 3485 w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & 3486 + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & 3487 ) / ( 3.0_wp * gg ) 3488 3489 IF ( k == nzt ) THEN 3490 w_int(n) = w_int_l 3491 ELSE 3492 w_int_u = ( ( gg-aa ) * w(k+1,j,i) + & 3493 ( gg-bb ) * w(k+1,j,i+1) + & 3494 ( gg-cc ) * w(k+1,j+1,i) + & 3495 ( gg-dd ) * w(k+1,j+1,i+1) & 3496 ) / ( 3.0_wp * gg ) 3497 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) * & 3498 ( w_int_u - w_int_l ) 3499 ENDIF 3500 3501 ELSE 3502 3503 w_int(n) = 0.0_wp 3504 3505 ENDIF 3506 3507 ENDDO 3508 3509 ENDDO 3508 ENDDO 3509 ! 3510 !-- Same procedure for interpolation of the v velocity-component 3511 i = ip + block_offset(nb)%i_off 3512 j = jp 3513 k = kp + block_offset(nb)%k_off 3514 3515 DO n = start_index(nb), end_index(nb) 3516 ! 3517 !-- Determine vertical index of topography top 3518 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 3519 3520 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3521 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3522 ! 3523 !-- Resolved-scale horizontal particle velocity is zero below z0. 3524 v_int(n) = 0.0_wp 3525 ELSE 3526 ! 3527 !-- Determine the sublayer. Further used as index. Please note, 3528 !-- logarithmus can not be reused from above, as in in case of 3529 !-- topography particle on u-grid can be above surface-layer height, 3530 !-- whereas it can be below on v-grid. 3531 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3532 * REAL( number_of_sublayers, KIND=wp ) & 3533 * d_z_p_z0 3534 ! 3535 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3536 !-- interpolate linearly between precalculated logarithm. 3537 log_z_z0_int = log_z_z0(INT(height_p)) & 3538 + ( height_p - INT(height_p) ) & 3539 * ( log_z_z0(INT(height_p)+1) & 3540 - log_z_z0(INT(height_p)) & 3541 ) 3542 ! 3543 !-- Get friction velocity and momentum flux from new surface data 3544 !-- types. 3545 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3546 surf_def_h(0)%end_index(jp,ip) ) THEN 3547 surf_start = surf_def_h(0)%start_index(jp,ip) 3548 !-- Limit friction velocity. In narrow canyons or holes the 3549 !-- friction velocity can become very small, resulting in a too 3550 !-- large particle speed. 3551 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3552 vsws_int = surf_def_h(0)%vsws(surf_start) 3553 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3554 surf_lsm_h%end_index(jp,ip) ) THEN 3555 surf_start = surf_lsm_h%start_index(jp,ip) 3556 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3557 vsws_int = surf_lsm_h%vsws(surf_start) 3558 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3559 surf_usm_h%end_index(jp,ip) ) THEN 3560 surf_start = surf_usm_h%start_index(jp,ip) 3561 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3562 vsws_int = surf_usm_h%vsws(surf_start) 3563 ENDIF 3564 ! 3565 !-- Neutral solution is applied for all situations, e.g. also for 3566 !-- unstable and stable situations. Even though this is not exact 3567 !-- this saves a lot of CPU time since several calls of intrinsic 3568 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3569 !-- as sensitivity studies revealed no significant effect of 3570 !-- using the neutral solution also for un/stable situations. 3571 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 3572 * log_z_z0_int - v_gtrans 3573 3574 ENDIF 3575 ELSE 3576 x = xv(n) + ( 0.5_wp - i ) * dx 3577 y = yv(n) - j * dy 3578 aa = x**2 + y**2 3579 bb = ( dx - x )**2 + y**2 3580 cc = x**2 + ( dy - y )**2 3581 dd = ( dx - x )**2 + ( dy - y )**2 3582 gg = aa + bb + cc + dd 3583 3584 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 3585 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 3586 ) / ( 3.0_wp * gg ) - v_gtrans 3587 3588 IF ( k == nzt ) THEN 3589 v_int(n) = v_int_l 3590 ELSE 3591 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 3592 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 3593 ) / ( 3.0_wp * gg ) - v_gtrans 3594 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3595 ( v_int_u - v_int_l ) 3596 ENDIF 3597 ENDIF 3598 ENDDO 3599 ! 3600 !-- Same procedure for interpolation of the w velocity-component 3601 i = ip + block_offset(nb)%i_off 3602 j = jp + block_offset(nb)%j_off 3603 k = kp - 1 3604 3605 DO n = start_index(nb), end_index(nb) 3606 IF ( vertical_particle_advection(particles(n)%group) ) THEN 3607 x = xv(n) + ( 0.5_wp - i ) * dx 3608 y = yv(n) + ( 0.5_wp - j ) * dy 3609 aa = x**2 + y**2 3610 bb = ( dx - x )**2 + y**2 3611 cc = x**2 + ( dy - y )**2 3612 dd = ( dx - x )**2 + ( dy - y )**2 3613 gg = aa + bb + cc + dd 3614 3615 w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & 3616 + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & 3617 ) / ( 3.0_wp * gg ) 3618 3619 IF ( k == nzt ) THEN 3620 w_int(n) = w_int_l 3621 ELSE 3622 w_int_u = ( ( gg-aa ) * w(k+1,j,i) + & 3623 ( gg-bb ) * w(k+1,j,i+1) + & 3624 ( gg-cc ) * w(k+1,j+1,i) + & 3625 ( gg-dd ) * w(k+1,j+1,i+1) & 3626 ) / ( 3.0_wp * gg ) 3627 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) * & 3628 ( w_int_u - w_int_l ) 3629 ENDIF 3630 ELSE 3631 w_int(n) = 0.0_wp 3632 ENDIF 3633 ENDDO 3634 ENDDO 3635 3636 CASE DEFAULT 3637 WRITE( message_string, * ) 'unknown particle velocity interpolation method = "', & 3638 TRIM( interpolation_method ), '"' 3639 CALL message( 'lpm_advec', 'PA0660', 1, 2, 0, 6, 0 ) 3640 3641 END SELECT 3510 3642 3511 3643 !-- Interpolate and calculate quantities needed for calculating the SGS … … 3819 3951 !-- and calculate SGS velocities again 3820 3952 dz_temp = zw(kp)-zw(kp-1) 3821 3953 3822 3954 DO nb = 0, 7 3823 3955 DO n = start_index(nb), end_index(nb) … … 3826 3958 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n)) .OR. & 3827 3959 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN 3828 3960 3829 3961 dt_particle(n) = 0.9_wp * MIN( & 3830 3962 ( dx / ABS( u_int(n) + rvar1_temp(n) ) ), & … … 3837 3969 rvar2_temp(n) = particles(n)%rvar2 3838 3970 rvar3_temp(n) = particles(n)%rvar3 3839 3971 3840 3972 de_dt_min = - e_int(n) / dt_particle(n) 3841 3973 … … 3857 3989 de_dz_int(n), de_dt, diss_int(n), & 3858 3990 dt_particle(n), rg(n,3), term_1_2(n) ) 3859 ENDIF 3860 3991 ENDIF 3992 3861 3993 ! 3862 3994 !-- Update particle velocites … … 3873 4005 ENDDO 3874 4006 ENDDO 3875 4007 3876 4008 ELSE 3877 4009 ! … … 3883 4015 3884 4016 dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio 3885 3886 4017 IF ( ANY( dens_ratio == 0.0_wp ) ) THEN 3887 DO nb = 0, 7 3888 DO n = start_index(nb), end_index(nb) 4018 ! 4019 !-- Decide whether the particle loop runs over the subboxes or only over 1, 4020 !-- number_of_particles. This depends on the selected interpolation method. 4021 !-- If particle interpolation method is not trilinear, then the sorting within 4022 !-- subboxes is not required. However, therefore the index start_index(nb) and 4023 !-- end_index(nb) are not defined and the loops are still over 4024 !-- number_of_particles. @todo find a more generic way to write this loop or 4025 !-- delete trilinear interpolation 4026 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4027 subbox_start = 0 4028 subbox_end = 7 4029 ELSE 4030 subbox_start = 1 4031 subbox_end = 1 4032 ENDIF 4033 ! 4034 !-- loop over subboxes. In case of simple interpolation scheme no subboxes 4035 !-- are introduced, as they are not required. Accordingly, this loops goes 4036 !-- from 1 to 1. 4037 DO nb = subbox_start, subbox_end 4038 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4039 particle_start = start_index(nb) 4040 particle_end = end_index(nb) 4041 ELSE 4042 particle_start = 1 4043 particle_end = number_of_particles 4044 ENDIF 4045 ! 4046 !-- Loop from particle start to particle end 4047 DO n = particle_start, particle_end 3889 4048 3890 4049 ! … … 3915 4074 IF ( cloud_droplets ) THEN 3916 4075 ! 3917 !-- Terminal velocity is computed for vertical direction (Rogers et 4076 !-- Terminal velocity is computed for vertical direction (Rogers et 3918 4077 !-- al., 1993, J. Appl. Meteorol.) 3919 4078 diameter = particles(n)%radius * 2000.0_wp !diameter in mm … … 3974 4133 ENDDO 3975 4134 ENDDO 3976 4135 3977 4136 ELSE 3978 3979 DO nb = 0, 7 3980 DO n = start_index(nb), end_index(nb) 4137 ! 4138 !-- Decide whether the particle loop runs over the subboxes or only over 1, 4139 !-- number_of_particles. This depends on the selected interpolation method. 4140 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4141 subbox_start = 0 4142 subbox_end = 7 4143 ELSE 4144 subbox_start = 1 4145 subbox_end = 1 4146 ENDIF 4147 !-- loop over subboxes. In case of simple interpolation scheme no subboxes 4148 !-- are introduced, as they are not required. Accordingly, this loops goes 4149 !-- from 1 to 1. 4150 DO nb = subbox_start, subbox_end 4151 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4152 particle_start = start_index(nb) 4153 particle_end = end_index(nb) 4154 ELSE 4155 particle_start = 1 4156 particle_end = number_of_particles 4157 ENDIF 4158 ! 4159 !-- Loop from particle start to particle end 4160 DO n = particle_start, particle_end 4161 3981 4162 ! 3982 4163 !-- Transport of particles with inertia … … 3988 4169 IF ( cloud_droplets ) THEN 3989 4170 ! 3990 !-- Terminal velocity is computed for vertical direction (Rogers et al., 4171 !-- Terminal velocity is computed for vertical direction (Rogers et al., 3991 4172 !-- 1993, J. Appl. Meteorol.) 3992 4173 diameter = particles(n)%radius * 2000.0_wp !diameter in mm … … 4054 4235 particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age 4055 4236 4056 DO nb = 0, 7 4057 DO n = start_index(nb), end_index(nb) 4237 ! 4238 !-- loop over subboxes. In case of simple interpolation scheme no subboxes 4239 !-- are introduced, as they are not required. Accordingly, this loops goes 4240 !-- from 1 to 1. 4241 ! 4242 !-- Decide whether the particle loop runs over the subboxes or only over 1, 4243 !-- number_of_particles. This depends on the selected interpolation method. 4244 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4245 subbox_start = 0 4246 subbox_end = 7 4247 ELSE 4248 subbox_start = 1 4249 subbox_end = 1 4250 ENDIF 4251 DO nb = subbox_start, subbox_end 4252 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 4253 particle_start = start_index(nb) 4254 particle_end = end_index(nb) 4255 ELSE 4256 particle_start = 1 4257 particle_end = number_of_particles 4258 ENDIF 4259 ! 4260 !-- Loop from particle start to particle end 4261 DO n = particle_start, particle_end 4058 4262 ! 4059 4263 !-- Increment the particle age and the total time that the particle … … 4063 4267 4064 4268 ! 4065 !-- Check whether there is still a particle that has not yet completed 4269 !-- Check whether there is still a particle that has not yet completed 4066 4270 !-- the total LES timestep 4067 4271 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN … … 4086 4290 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n, & 4087 4291 dt_n, rg_n, fac ) 4088 4292 4089 4293 REAL(wp) :: a1 !< dummy argument 4090 4294 REAL(wp) :: dedt_n !< time derivative of TKE at particle position … … 4161 4365 INTEGER(iwp), INTENT(IN) :: j !< grid index of particle box along y 4162 4366 INTEGER(iwp), INTENT(IN) :: k !< grid index of particle box along z 4163 4367 4164 4368 INTEGER(iwp) :: inc !< dummy for sorting algorithmus 4165 4369 INTEGER(iwp) :: ir !< dummy for sorting algorithmus … … 4184 4388 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions 4185 4389 INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions 4186 4390 4187 4391 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary 4188 4392 LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary … … 4225 4429 CASE ( 'bottom/top' ) 4226 4430 4227 ! IF ( location_bc == 'bottom/top' ) THEN4228 4229 4431 ! 4230 4432 !-- Apply boundary conditions to those particles that have crossed the top or … … 4270 4472 ENDIF 4271 4473 ENDIF 4272 4474 4273 4475 IF ( particles(n)%z < zw(0) .AND. particles(n)%particle_mask ) THEN 4274 4476 IF ( ibc_par_b == 1 ) THEN … … 4290 4492 ENDDO 4291 4493 4292 ! ELSEIF ( location_bc == 'walls' ) THEN4293 4494 CASE ( 'walls' ) 4294 4495 … … 4347 4548 ELSE 4348 4549 zwall = zw(k-1) 4349 ENDIF 4550 ENDIF 4350 4551 ! 4351 4552 !-- Initialize flags to check if particle reflection is necessary … … 4418 4619 MIN( prt_z - pos_z_old, -1E-30_wp ), & 4419 4620 prt_z > pos_z_old ) 4420 4621 4421 4622 x_ind(t_index) = i1 4422 4623 y_ind(t_index) = j1 … … 4429 4630 cross_wall_z = .TRUE. 4430 4631 ENDIF 4431 4632 4432 4633 t_index_number = t_index - 1 4433 4634 ! … … 4484 4685 t_old = 0.0_wp 4485 4686 DO t_index = 1, t_index_number 4486 ! 4687 ! 4487 4688 !-- Calculate intermediate particle position according to the 4488 4689 !-- timesteps a particle reaches any wall. … … 4583 4784 y_ind(t_index:t_index_number) = j2 4584 4785 ENDIF !particle reflection in y direction done 4585 4786 4586 4787 ! 4587 4788 !-- Check if a particle needs to be reflected at any xy-wall. If … … 4623 4824 z_ind(t_index:t_index_number) = k2 4624 4825 ENDIF !particle reflection in z direction done 4625 4826 4626 4827 ! 4627 4828 !-- Swap time … … 4655 4856 4656 4857 END SUBROUTINE lpm_boundary_conds 4657 4658 4858 4859 4860 !------------------------------------------------------------------------------! 4861 ! Description: 4862 ! ------------ 4863 !> Calculates change in droplet radius by condensation/evaporation, using 4864 !> either an analytic formula or by numerically integrating the radius growth 4865 !> equation including curvature and solution effects using Rosenbrocks method 4866 !> (see Numerical recipes in FORTRAN, 2nd edition, p. 731). 4867 !> The analytical formula and growth equation follow those given in 4868 !> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103). 4869 !------------------------------------------------------------------------------! 4659 4870 SUBROUTINE lpm_droplet_condensation (i,j,k) 4660 4871 … … 4812 5023 DO 4813 5024 4814 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -&5025 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 4815 5026 afactor / r_ros + & 4816 5027 bfactor / r_ros**3 & … … 4818 5029 4819 5030 d2rdtdr = -ddenom * ventilation_effect(n) * ( & 4820 (e_a / e_s - 1.0) * r_ros**4 - &4821 afactor * r0 * r_ros**2 -&4822 2.0 * afactor * r_ros**3 +&4823 3.0 * bfactor * r0 +&4824 4.0 * bfactor * r_ros&5031 (e_a / e_s - 1.0_wp ) * r_ros**4 - & 5032 afactor * r0 * r_ros**2 - & 5033 2.0_wp * afactor * r_ros**3 + & 5034 3.0_wp * bfactor * r0 + & 5035 4.0_wp * bfactor * r_ros & 4825 5036 ) & 4826 5037 / ( r_ros**4 * ( r_ros + r0 )**2 ) 4827 5038 4828 k1 = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )5039 k1 = drdt / ( 1.0_wp - gamma * dt_ros * d2rdtdr ) 4829 5040 4830 5041 r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1) 4831 5042 r_err = r_ros 4832 5043 4833 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -&4834 afactor / r_ros + &4835 bfactor / r_ros**3 &5044 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 5045 afactor / r_ros + & 5046 bfactor / r_ros**3 & 4836 5047 ) / ( r_ros + r0 ) 4837 5048 4838 5049 k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / & 4839 ( 1.0 - dt_ros * gamma * d2rdtdr )4840 4841 r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5* k2), particles(n)%aux1)5050 ( 1.0_wp - dt_ros * gamma * d2rdtdr ) 5051 5052 r_ros = MAX(r_ros_ini + dt_ros * ( 1.5_wp * k1 + 0.5_wp * k2), particles(n)%aux1) 4842 5053 ! 4843 5054 !-- Check error of the solution, and reduce dt_ros if necessary. … … 4909 5120 particles(n)%class = MAX( particles(n)%class, 1 ) 4910 5121 ENDIF 4911 4912 !--Store new radius to particle features5122 ! 5123 !-- Store new radius to particle features 4913 5124 particles(n)%radius = new_r(n) 4914 5125 … … 5134 5345 aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * & 5135 5346 particles(1:number_of_particles)%aux1**3 * & 5136 4.0 / 3.0* pi * rho_s5347 4.0_wp / 3.0_wp * pi * rho_s 5137 5348 ENDIF 5138 5349 ! … … 5144 5355 !-- For collisions, the weighting factor of at least one super-droplet 5145 5356 !-- needs to be larger or equal to one. 5146 IF ( MIN( weight(n), weight(m) ) .LT. 1.0 ) CYCLE5357 IF ( MIN( weight(n), weight(m) ) .LT. 1.0_wp ) CYCLE 5147 5358 ! 5148 5359 !-- Get mass of individual droplets (aerosols) … … 5175 5386 ENDIF 5176 5387 5177 IF ( collection_probability .GT. 0.0 ) THEN5388 IF ( collection_probability .GT. 0.0_wp ) THEN 5178 5389 ! 5179 5390 !-- Super-droplet n collects droplets of super-droplet m … … 5363 5574 ( hkernel(i,j), i = 1,radius_classes ) 5364 5575 ENDDO 5365 5576 5366 5577 DO k = 1, dissipation_classes 5367 5578 DO i = 1, radius_classes … … 6154 6365 ql(k,j,i) < ql_crit ) CYCLE 6155 6366 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 6156 ! 6367 ! 6157 6368 !-- Start splitting operations. Each particle is checked if it 6158 6369 !-- fulfilled the splitting criterion's. In splitting mode 'const' … … 6168 6379 particles(n)%radius >= radius_split .AND. & 6169 6380 particles(n)%weight_factor >= weight_factor_split ) & 6170 THEN 6381 THEN 6171 6382 ! 6172 6383 !-- Calculate the new number of particles. 6173 new_size = prt_count(k,j,i) + splitting_factor - 1 6384 new_size = prt_count(k,j,i) + splitting_factor - 1 6174 6385 ! 6175 6386 !-- Cycle if maximum number of particles per grid box 6176 6387 !-- is greater than the allowed maximum number. 6177 IF ( new_size >= max_number_particles_per_gridbox ) CYCLE 6388 IF ( new_size >= max_number_particles_per_gridbox ) CYCLE 6178 6389 ! 6179 6390 !-- Reallocate particle array if necessary. … … 6191 6402 DO jpp = 1, splitting_factor-1 6192 6403 grid_particles(k,j,i)%particles(jpp+old_size) = & 6193 tmp_particle 6404 tmp_particle 6194 6405 ENDDO 6195 6406 new_particles_gb = new_particles_gb + splitting_factor - 1 … … 6197 6408 !-- Save the new number of super droplets for every grid box. 6198 6409 prt_count(k,j,i) = prt_count(k,j,i) + & 6199 splitting_factor - 1 6410 splitting_factor - 1 6200 6411 ENDIF 6201 6412 ENDDO … … 6216 6427 m3 = 0.0_wp 6217 6428 m3_total = 0.0_wp 6218 nr = 0.0_wp 6429 nr = 0.0_wp 6219 6430 nrclgb = 0.0_wp 6220 nrclgb_total = 0.0_wp 6431 nrclgb_total = 0.0_wp 6221 6432 nr_total = 0.0_wp 6222 6433 rm = 0.0_wp 6223 6434 rm_total = 0.0_wp 6224 6435 6225 6436 DO i = nxl, nxr 6226 6437 DO j = nys, nyn … … 6231 6442 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 6232 6443 nrclgb = nrclgb + 1.0_wp 6233 ! 6234 !-- Calculate moments of DSD. 6444 ! 6445 !-- Calculate moments of DSD. 6235 6446 DO n = 1, number_of_particles 6236 6447 IF ( particles(n)%particle_mask .AND. & … … 6242 6453 particles(n)%weight_factor 6243 6454 IF ( isf == 1 ) THEN 6244 diameter = particles(n)%radius * 2.0_wp 6455 diameter = particles(n)%radius * 2.0_wp 6245 6456 lwc = lwc + factor_volume_to_mass * & 6246 6457 particles(n)%radius**3 * & 6247 6458 particles(n)%weight_factor 6248 m1 = m1 + particles(n)%weight_factor * diameter 6249 m2 = m2 + particles(n)%weight_factor * diameter**2 6459 m1 = m1 + particles(n)%weight_factor * diameter 6460 m2 = m2 + particles(n)%weight_factor * diameter**2 6250 6461 m3 = m3 + particles(n)%weight_factor * diameter**3 6251 ENDIF 6252 ENDIF 6462 ENDIF 6463 ENDIF 6253 6464 ENDDO 6254 6465 ENDDO … … 6287 6498 ( nr_total * factor_volume_to_mass ) & 6288 6499 )**0.3333333_wp, 0.0_wp, nrclgb_total > 0.0_wp & 6289 ) 6500 ) 6290 6501 ! 6291 6502 !-- Check which function should be used to approximate the DSD. … … 6303 6514 0.0_wp, nrclgb_total > 0.0_wp & 6304 6515 ) 6305 zeta = m1_total * m3_total / m2_total**2 6516 zeta = m1_total * m3_total / m2_total**2 6306 6517 mu = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) / & 6307 6518 ( zeta - 1.0_wp ), 0.0_wp & … … 6309 6520 6310 6521 lambda = ( pirho_l * nr_total / lwc_total * & 6311 ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp ) & 6522 ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp ) & 6312 6523 )**0.3333333_wp 6313 6524 nr0 = nr_total / gamma( mu + 1.0_wp ) * lambda**( mu + 1.0_wp ) 6314 6525 6315 6526 DO n = 0, n_max-1 6316 diameter = r_bin_mid(n) * 2.0_wp 6527 diameter = r_bin_mid(n) * 2.0_wp 6317 6528 an_spl(n) = nr0 * diameter**mu * EXP( -lambda * diameter ) * & 6318 6529 ( r_bin(n+1) - r_bin(n) ) * 2.0_wp … … 6338 6549 !-- Criterion to avoid super droplets with a weighting factor < 1.0. 6339 6550 an_spl = MAX(an_spl, 1.0_wp) 6340 6551 6341 6552 DO i = nxl, nxr 6342 6553 DO j = nys, nyn 6343 6554 DO k = nzb+1, nzt 6344 6555 number_of_particles = prt_count(k,j,i) 6345 IF ( number_of_particles <= 0 .OR. & 6556 IF ( number_of_particles <= 0 .OR. & 6346 6557 ql(k,j,i) < ql_crit ) CYCLE 6347 6558 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 6348 new_particles_gb = 0 6349 ! 6559 new_particles_gb = 0 6560 ! 6350 6561 !-- Start splitting operations. Each particle is checked if it 6351 !-- fulfilled the splitting criterion's. In splitting mode 'cl_av' 6352 !-- a critical radius (radius_split) and a splitting function must 6562 !-- fulfilled the splitting criterion's. In splitting mode 'cl_av' 6563 !-- a critical radius (radius_split) and a splitting function must 6353 6564 !-- be prescribed (see particles_par). The critical weighting factor 6354 6565 !-- is calculated while approximating a 'gamma', 'log' or 'exp'- 6355 !-- drop size distribution. In this mode the DSD is calculated as 6356 !-- an average over all cloudy grid boxes. Super droplets which 6566 !-- drop size distribution. In this mode the DSD is calculated as 6567 !-- an average over all cloudy grid boxes. Super droplets which 6357 6568 !-- have a larger radius and larger weighting factor are split into 6358 !-- 'splitting_factor' super droplets. In this case the splitting 6359 !-- factor is calculated of weighting factor of the super droplet 6360 !-- and the approximated number concentration for droplet of such 6361 !-- a size. Due to the splitting, the weighting factor of the 6362 !-- super droplet and all created clones is reduced by the factor 6569 !-- 'splitting_factor' super droplets. In this case the splitting 6570 !-- factor is calculated of weighting factor of the super droplet 6571 !-- and the approximated number concentration for droplet of such 6572 !-- a size. Due to the splitting, the weighting factor of the 6573 !-- super droplet and all created clones is reduced by the factor 6363 6574 !-- of 'splitting_facor'. 6364 6575 DO n = 1, number_of_particles … … 6379 6590 IF ( splitting_factor < 2 ) CYCLE 6380 6591 ! 6381 !-- Calculate the new number of particles. 6592 !-- Calculate the new number of particles. 6382 6593 new_size = prt_count(k,j,i) + splitting_factor - 1 6383 6594 ! 6384 6595 !-- Cycle if maximum number of particles per grid box 6385 !-- is greater than the allowed maximum number. 6596 !-- is greater than the allowed maximum number. 6386 6597 IF ( new_size >= max_number_particles_per_gridbox ) & 6387 6598 CYCLE … … 6391 6602 CALL realloc_particles_array(i,j,k,new_size) 6392 6603 ENDIF 6393 old_size = prt_count(k,j,i) 6604 old_size = prt_count(k,j,i) 6394 6605 new_particles_gb = new_particles_gb + & 6395 6606 splitting_factor - 1 … … 6400 6611 tmp_particle = particles(n) 6401 6612 ! 6402 !-- Create splitting_factor-1 new particles. 6613 !-- Create splitting_factor-1 new particles. 6403 6614 DO jpp = 1, splitting_factor-1 6404 6615 grid_particles(k,j,i)%particles(jpp+old_size) = & 6405 tmp_particle 6616 tmp_particle 6406 6617 ENDDO 6407 ! 6618 ! 6408 6619 !-- Save the new number of super droplets. 6409 6620 prt_count(k,j,i) = prt_count(k,j,i) + & … … 6412 6623 ENDDO 6413 6624 ENDDO 6414 6625 6415 6626 ENDDO 6416 6627 ENDDO … … 6422 6633 DO j = nys, nyn 6423 6634 DO k = nzb+1, nzt 6424 6425 ! 6426 !-- Initialize summing variables. 6635 6636 ! 6637 !-- Initialize summing variables. 6427 6638 lwc = 0.0_wp 6428 6639 m1 = 0.0_wp … … 6431 6642 nr = 0.0_wp 6432 6643 rm = 0.0_wp 6433 6644 6434 6645 new_particles_gb = 0 6435 6646 number_of_particles = prt_count(k,j,i) … … 6437 6648 ql(k,j,i) < ql_crit ) CYCLE 6438 6649 particles => grid_particles(k,j,i)%particles 6439 ! 6440 !-- Calculate moments of DSD. 6650 ! 6651 !-- Calculate moments of DSD. 6441 6652 DO n = 1, number_of_particles 6442 6653 IF ( particles(n)%particle_mask .AND. & … … 6447 6658 particles(n)%radius**3 * & 6448 6659 particles(n)%weight_factor 6449 IF ( isf == 1 ) THEN 6450 diameter = particles(n)%radius * 2.0_wp 6660 IF ( isf == 1 ) THEN 6661 diameter = particles(n)%radius * 2.0_wp 6451 6662 lwc = lwc + factor_volume_to_mass * & 6452 particles(n)%radius**3 * & 6663 particles(n)%radius**3 * & 6453 6664 particles(n)%weight_factor 6454 6665 m1 = m1 + particles(n)%weight_factor * diameter 6455 6666 m2 = m2 + particles(n)%weight_factor * diameter**2 6456 6667 m3 = m3 + particles(n)%weight_factor * diameter**3 6457 ENDIF 6458 ENDIF 6459 ENDDO 6460 6461 IF ( nr <= 0.0 .OR. rm <= 0.0_wp ) CYCLE6462 ! 6463 !-- Calculate mean volume averaged radius. 6668 ENDIF 6669 ENDIF 6670 ENDDO 6671 6672 IF ( nr <= 0.0_wp .OR. rm <= 0.0_wp ) CYCLE 6673 ! 6674 !-- Calculate mean volume averaged radius. 6464 6675 rm = ( rm / ( nr * factor_volume_to_mass ) )**0.3333333_wp 6465 6676 ! 6466 !-- Check which function should be used to approximate the DSD. 6677 !-- Check which function should be used to approximate the DSD. 6467 6678 IF ( isf == 1 ) THEN 6468 6679 ! … … 6475 6686 lambda = ( pirho_l * nr / lwc * & 6476 6687 ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * & 6477 ( mu + 1.0_wp ) & 6688 ( mu + 1.0_wp ) & 6478 6689 )**0.3333333_wp 6479 6690 nr0 = ( nr / (gamma( mu + 1.0_wp ) ) ) * & … … 6481 6692 6482 6693 DO n = 0, n_max-1 6483 diameter = r_bin_mid(n) * 2.0_wp 6694 diameter = r_bin_mid(n) * 2.0_wp 6484 6695 an_spl(n) = nr0 * diameter**mu * & 6485 6696 EXP( -lambda * diameter ) * & … … 6489 6700 ! 6490 6701 !-- Lognormal size distribution to calculate critical 6491 !-- weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974). 6702 !-- weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974). 6492 6703 DO n = 0, n_max-1 6493 6704 an_spl(n) = nr / ( SQRT( 2.0_wp * pi ) * & … … 6495 6706 ) * & 6496 6707 EXP( -( LOG( r_bin_mid(n) / rm )**2 ) / & 6497 ( 2.0_wp * LOG(sigma_log)**2 ) & 6498 ) * & 6708 ( 2.0_wp * LOG(sigma_log)**2 ) & 6709 ) * & 6499 6710 ( r_bin(n+1) - r_bin(n) ) 6500 6711 ENDDO … … 6509 6720 ENDDO 6510 6721 ENDIF 6511 6512 ! 6513 !-- Criterion to avoid super droplets with a weighting factor < 1.0. 6722 6723 ! 6724 !-- Criterion to avoid super droplets with a weighting factor < 1.0. 6514 6725 an_spl = MAX(an_spl, 1.0_wp) 6515 ! 6726 ! 6516 6727 !-- Start splitting operations. Each particle is checked if it 6517 !-- fulfilled the splitting criterion's. In splitting mode 'gb_av' 6728 !-- fulfilled the splitting criterion's. In splitting mode 'gb_av' 6518 6729 !-- a critical radius (radius_split) and a splitting function must 6519 6730 !-- be prescribed (see particles_par). The critical weighting factor … … 6546 6757 6547 6758 ! 6548 !-- Calculate the new number of particles. 6759 !-- Calculate the new number of particles. 6549 6760 new_size = prt_count(k,j,i) + splitting_factor - 1 6550 6761 ! 6551 6762 !-- Cycle if maximum number of particles per grid box 6552 !-- is greater than the allowed maximum number. 6763 !-- is greater than the allowed maximum number. 6553 6764 IF ( new_size >= max_number_particles_per_gridbox ) & 6554 6765 CYCLE 6555 6766 ! 6556 !-- Reallocate particle array if necessary. 6767 !-- Reallocate particle array if necessary. 6557 6768 IF ( new_size > SIZE(particles) ) THEN 6558 6769 CALL realloc_particles_array(i,j,k,new_size) … … 6575 6786 splitting_factor - 1 6576 6787 new_particles_gb = new_particles_gb + & 6577 splitting_factor - 1 6788 splitting_factor - 1 6578 6789 ENDIF 6579 6790 ENDDO … … 6602 6813 SUBROUTINE lpm_merging 6603 6814 6604 INTEGER(iwp) :: i !< 6605 INTEGER(iwp) :: j !< 6606 INTEGER(iwp) :: k !< 6815 INTEGER(iwp) :: i !< 6816 INTEGER(iwp) :: j !< 6817 INTEGER(iwp) :: k !< 6607 6818 INTEGER(iwp) :: n !< 6608 6819 INTEGER(iwp) :: merge_drp = 0 !< number of merged droplets 6609 6610 6820 6821 6611 6822 REAL(wp) :: ql_crit = 1.0E-5_wp !< threshold lwc for cloudy grid cells 6612 6823 !< (e.g. Siebesma et al 2003, JAS, 60) 6613 6824 6614 6825 CALL cpu_log( log_point_s(81), 'lpm_merging', 'start' ) 6615 6826 6616 6827 merge_drp = 0 6617 6828 6618 6829 IF ( weight_factor_merge == -1.0_wp ) THEN 6619 6830 weight_factor_merge = 0.5_wp * initial_weighting_factor … … 6623 6834 DO j = nys, nyn 6624 6835 DO k = nzb+1, nzt 6625 6836 6626 6837 number_of_particles = prt_count(k,j,i) 6627 6838 IF ( number_of_particles <= 0 .OR. & … … 6635 6846 !-- particle array. Tests showed that this simplified method can be 6636 6847 !-- used because it will only take place outside of cloudy grid 6637 !-- boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned 6848 !-- boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned 6638 6849 !-- and subsequent evaporated super droplets will be merged. 6639 6850 DO n = 1, number_of_particles-1 … … 6641 6852 particles(n+1)%particle_mask .AND. & 6642 6853 particles(n)%radius <= radius_merge .AND. & 6643 particles(n)%weight_factor <= weight_factor_merge ) & 6854 particles(n)%weight_factor <= weight_factor_merge ) & 6644 6855 THEN 6645 6856 particles(n+1)%weight_factor = & … … 6652 6863 deleted_particles = deleted_particles + 1 6653 6864 merge_drp = merge_drp + 1 6654 6865 6655 6866 ENDIF 6656 6867 ENDDO … … 6658 6869 ENDDO 6659 6870 ENDDO 6660 6661 6871 6872 6662 6873 CALL cpu_log( log_point_s(81), 'lpm_merging', 'stop' ) 6663 6874 … … 7145 7356 7146 7357 ALLOCATE(rvnp(MAX(1,trnp_count_recv))) 7147 ! 7358 ! 7148 7359 !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit 7149 7360 !-- variables in structure particle_type (due to the calculation of par_size) … … 7164 7375 7165 7376 ALLOCATE(rvsp(MAX(1,trsp_count_recv))) 7166 ! 7377 ! 7167 7378 !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit 7168 7379 !-- variables in structure particle_type (due to the calculation of par_size) … … 7516 7727 IF ( np_before_move <= 0 ) CYCLE 7517 7728 particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move) 7518 7729 7519 7730 DO n = 1, np_before_move 7520 7731 i = particles_before_move(n)%x * ddx … … 7578 7789 !------------------------------------------------------------------------------! 7579 7790 SUBROUTINE lpm_check_cfl 7580 7791 7581 7792 IMPLICIT NONE 7582 7793 7583 7794 INTEGER(iwp) :: i !< running index, x-direction 7584 7795 INTEGER(iwp) :: j !< running index, y-direction … … 7759 7970 ! ----------- 7760 7971 !> Routine for the whole processor 7761 !> Sort all particles into the 8 respective subgrid boxes 7972 !> Sort all particles into the 8 respective subgrid boxes and free space of 7973 !> particles which has been marked for deletion 7762 7974 !------------------------------------------------------------------------------! 7763 SUBROUTINE lpm_sort_ in_subboxes7975 SUBROUTINE lpm_sort_and_delete 7764 7976 7765 7977 INTEGER(iwp) :: i !< … … 7776 7988 INTEGER(iwp), DIMENSION(0:7) :: sort_count !< 7777 7989 7778 TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE :: sort_particles !< 7779 7780 CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'start' ) 7781 DO ip = nxl, nxr 7782 DO jp = nys, nyn 7783 DO kp = nzb+1, nzt 7784 number_of_particles = prt_count(kp,jp,ip) 7785 IF ( number_of_particles <= 0 ) CYCLE 7786 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 7787 7788 nn = 0 7789 sort_count = 0 7790 ALLOCATE( sort_particles(number_of_particles, 0:7) ) 7791 7792 DO n = 1, number_of_particles 7793 sort_index = 0 7794 7795 IF ( particles(n)%particle_mask ) THEN 7796 nn = nn + 1 7797 ! 7798 !-- Sorting particles with a binary scheme 7799 !-- sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox 7800 !-- sort_index=000_2=0_10 -> particle at the right,north,top subgridbox 7801 !-- For this the center of the gridbox is calculated 7802 i = (particles(n)%x + 0.5_wp * dx) * ddx 7803 j = (particles(n)%y + 0.5_wp * dy) * ddy 7804 7805 IF ( i == ip ) sort_index = sort_index + 4 7806 IF ( j == jp ) sort_index = sort_index + 2 7807 IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1 7808 7809 sort_count(sort_index) = sort_count(sort_index) + 1 7810 m = sort_count(sort_index) 7811 sort_particles(m,sort_index) = particles(n) 7812 sort_particles(m,sort_index)%block_nr = sort_index 7813 ENDIF 7990 TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE :: sort_particles !< 7991 7992 CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' ) 7993 IF ( TRIM(particle_interpolation) == 'trilinear' ) THEN 7994 DO ip = nxl, nxr 7995 DO jp = nys, nyn 7996 DO kp = nzb+1, nzt 7997 number_of_particles = prt_count(kp,jp,ip) 7998 IF ( number_of_particles <= 0 ) CYCLE 7999 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 8000 nn = 0 8001 sort_count = 0 8002 ALLOCATE( sort_particles(number_of_particles, 0:7) ) 8003 8004 DO n = 1, number_of_particles 8005 sort_index = 0 8006 8007 IF ( particles(n)%particle_mask ) THEN 8008 nn = nn + 1 8009 ! 8010 !-- Sorting particles with a binary scheme 8011 !-- sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox 8012 !-- sort_index=000_2=0_10 -> particle at the right,north,top subgridbox 8013 !-- For this the center of the gridbox is calculated 8014 i = (particles(n)%x + 0.5_wp * dx) * ddx 8015 j = (particles(n)%y + 0.5_wp * dy) * ddy 8016 8017 IF ( i == ip ) sort_index = sort_index + 4 8018 IF ( j == jp ) sort_index = sort_index + 2 8019 IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1 8020 8021 sort_count(sort_index) = sort_count(sort_index) + 1 8022 m = sort_count(sort_index) 8023 sort_particles(m,sort_index) = particles(n) 8024 sort_particles(m,sort_index)%block_nr = sort_index 8025 ENDIF 8026 ENDDO 8027 8028 nn = 0 8029 DO is = 0,7 8030 grid_particles(kp,jp,ip)%start_index(is) = nn + 1 8031 DO n = 1,sort_count(is) 8032 nn = nn + 1 8033 particles(nn) = sort_particles(n,is) 8034 ENDDO 8035 grid_particles(kp,jp,ip)%end_index(is) = nn 8036 ENDDO 8037 8038 number_of_particles = nn 8039 prt_count(kp,jp,ip) = number_of_particles 8040 DEALLOCATE(sort_particles) 7814 8041 ENDDO 7815 7816 nn = 07817 DO is = 0,77818 grid_particles(kp,jp,ip)%start_index(is) = nn + 17819 DO n = 1,sort_count(is)7820 nn = nn + 17821 particles(nn) = sort_particles(n,is)7822 ENDDO7823 grid_particles(kp,jp,ip)%end_index(is) = nn7824 ENDDO7825 7826 number_of_particles = nn7827 prt_count(kp,jp,ip) = number_of_particles7828 DEALLOCATE(sort_particles)7829 8042 ENDDO 7830 8043 ENDDO 7831 ENDDO 7832 CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'stop' ) 7833 7834 END SUBROUTINE lpm_sort_in_subboxes 8044 8045 !-- In case of the simple interpolation method the particles must not 8046 !-- be sorted in subboxes but particles marked for deletion must be 8047 !-- deleted and number of particles must be recalculated 8048 ELSE 8049 8050 DO ip = nxl, nxr 8051 DO jp = nys, nyn 8052 DO kp = nzb+1, nzt 8053 8054 number_of_particles = prt_count(kp,jp,ip) 8055 IF ( number_of_particles <= 0 ) CYCLE 8056 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 8057 ! 8058 !-- repack particles array, i.e. delete particles and recalculate 8059 !-- number of particles 8060 CALL lpm_pack 8061 prt_count(kp,jp,ip) = number_of_particles 8062 ENDDO 8063 ENDDO 8064 ENDDO 8065 ENDIF 8066 CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'stop' ) 8067 8068 END SUBROUTINE lpm_sort_and_delete 7835 8069 7836 8070 … … 7841 8075 !------------------------------------------------------------------------------! 7842 8076 SUBROUTINE lpm_pack 7843 8077 7844 8078 INTEGER(iwp) :: n !< 7845 8079 INTEGER(iwp) :: nn !< … … 7874 8108 7875 8109 END SUBROUTINE lpm_pack 7876 7877 8110 8111 7878 8112 !------------------------------------------------------------------------------! 7879 8113 ! Description: … … 7904 8138 number_of_particles = prt_count(k,j,i) 7905 8139 IF ( number_of_particles <= 0 ) CYCLE 7906 7907 8140 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 7908 8141
Note: See TracChangeset
for help on using the changeset viewer.
