[1873] | 1 | !> @file lpm_init.f90 |
---|
[2000] | 2 | !------------------------------------------------------------------------------! |
---|
[1036] | 3 | ! This file is part of PALM. |
---|
| 4 | ! |
---|
[2000] | 5 | ! PALM is free software: you can redistribute it and/or modify it under the |
---|
| 6 | ! terms of the GNU General Public License as published by the Free Software |
---|
| 7 | ! Foundation, either version 3 of the License, or (at your option) any later |
---|
| 8 | ! version. |
---|
[1036] | 9 | ! |
---|
| 10 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 11 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 12 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 13 | ! |
---|
| 14 | ! You should have received a copy of the GNU General Public License along with |
---|
| 15 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 16 | ! |
---|
[2101] | 17 | ! Copyright 1997-2017 Leibniz Universitaet Hannover |
---|
[2000] | 18 | !------------------------------------------------------------------------------! |
---|
[1036] | 19 | ! |
---|
[254] | 20 | ! Current revisions: |
---|
[1] | 21 | ! ----------------- |
---|
[2318] | 22 | ! |
---|
| 23 | ! |
---|
[1930] | 24 | ! Former revisions: |
---|
| 25 | ! ----------------- |
---|
| 26 | ! $Id: lpm_init.f90 2318 2017-07-20 17:27:44Z suehring $ |
---|
[2318] | 27 | ! Get topography top index via Function call |
---|
| 28 | ! |
---|
| 29 | ! 2317 2017-07-20 17:27:19Z suehring |
---|
[2312] | 30 | ! Extended particle data type. Aerosol initialization improved. |
---|
| 31 | ! |
---|
| 32 | ! 2305 2017-07-06 11:18:47Z hoffmann |
---|
[2305] | 33 | ! Improved calculation of particle IDs. |
---|
[2312] | 34 | ! |
---|
[2305] | 35 | ! 2274 2017-06-09 13:27:48Z Giersch |
---|
[2274] | 36 | ! Changed error messages |
---|
[2312] | 37 | ! |
---|
[2274] | 38 | ! 2265 2017-06-08 16:58:28Z schwenkel |
---|
[2265] | 39 | ! Unused variables removed. |
---|
[2312] | 40 | ! |
---|
[2265] | 41 | ! 2263 2017-06-08 14:59:01Z schwenkel |
---|
[2263] | 42 | ! Implemented splitting and merging algorithm |
---|
[2312] | 43 | ! |
---|
[2263] | 44 | ! 2233 2017-05-30 18:08:54Z suehring |
---|
[1930] | 45 | ! |
---|
[2233] | 46 | ! 2232 2017-05-30 17:47:52Z suehring |
---|
| 47 | ! Adjustments according to new topography realization |
---|
[2312] | 48 | ! |
---|
| 49 | ! |
---|
[2224] | 50 | ! 2223 2017-05-15 16:38:09Z suehring |
---|
| 51 | ! Add check for particle release at model top |
---|
[2312] | 52 | ! |
---|
[2183] | 53 | ! 2182 2017-03-17 14:27:40Z schwenkel |
---|
| 54 | ! Added parameters for simplified particle initialization. |
---|
[2305] | 55 | ! |
---|
[2123] | 56 | ! 2122 2017-01-18 12:22:54Z hoffmann |
---|
[2312] | 57 | ! Improved initialization of equilibrium aerosol radii |
---|
[2123] | 58 | ! Calculation of particle ID |
---|
| 59 | ! |
---|
[2001] | 60 | ! 2000 2016-08-20 18:09:15Z knoop |
---|
| 61 | ! Forced header and separation lines into 80 columns |
---|
[2312] | 62 | ! |
---|
[1930] | 63 | ! 2016-06-09 16:25:25Z suehring |
---|
[2312] | 64 | ! Bugfix in determining initial particle height and grid index in case of |
---|
[1929] | 65 | ! seed_follows_topography. |
---|
[2312] | 66 | ! Bugfix concerning random positions, ensure that particles do not move more |
---|
[1929] | 67 | ! than one grid length. |
---|
| 68 | ! Bugfix logarithmic interpolation. |
---|
| 69 | ! Initial setting of sgs_wf_part. |
---|
[1321] | 70 | ! |
---|
[1891] | 71 | ! 1890 2016-04-22 08:52:11Z hoffmann |
---|
[2312] | 72 | ! Initialization of aerosol equilibrium radius not possible in supersaturated |
---|
[1891] | 73 | ! environments. Therefore, a maximum supersaturation of -1 % is assumed during |
---|
| 74 | ! initialization. |
---|
| 75 | ! |
---|
[1874] | 76 | ! 1873 2016-04-18 14:50:06Z maronga |
---|
[1929] | 77 | ! Module renamed (removed _mod |
---|
[2312] | 78 | ! |
---|
[1872] | 79 | ! 1871 2016-04-15 11:46:09Z hoffmann |
---|
| 80 | ! Initialization of aerosols added. |
---|
| 81 | ! |
---|
[1851] | 82 | ! 1850 2016-04-08 13:29:27Z maronga |
---|
| 83 | ! Module renamed |
---|
| 84 | ! |
---|
[1832] | 85 | ! 1831 2016-04-07 13:15:51Z hoffmann |
---|
| 86 | ! curvature_solution_effects moved to particle_attributes |
---|
| 87 | ! |
---|
[1823] | 88 | ! 1822 2016-04-07 07:49:42Z hoffmann |
---|
| 89 | ! Unused variables removed. |
---|
| 90 | ! |
---|
[1784] | 91 | ! 1783 2016-03-06 18:36:17Z raasch |
---|
| 92 | ! netcdf module added |
---|
| 93 | ! |
---|
[2312] | 94 | ! 1725 2015-11-17 13:01:51Z hoffmann |
---|
| 95 | ! Bugfix: Processor-dependent seed for random function is generated before it is |
---|
[1726] | 96 | ! used. |
---|
| 97 | ! |
---|
[1692] | 98 | ! 1691 2015-10-26 16:17:44Z maronga |
---|
| 99 | ! Renamed prandtl_layer to constant_flux_layer. |
---|
| 100 | ! |
---|
[1686] | 101 | ! 1685 2015-10-08 07:32:13Z raasch |
---|
| 102 | ! bugfix concerning vertical index offset in case of ocean |
---|
| 103 | ! |
---|
[1683] | 104 | ! 1682 2015-10-07 23:56:08Z knoop |
---|
[2312] | 105 | ! Code annotations made doxygen readable |
---|
[1683] | 106 | ! |
---|
[1576] | 107 | ! 1575 2015-03-27 09:56:27Z raasch |
---|
| 108 | ! initial vertical particle position is allowed to follow the topography |
---|
| 109 | ! |
---|
[1360] | 110 | ! 1359 2014-04-11 17:15:14Z hoffmann |
---|
[2312] | 111 | ! New particle structure integrated. |
---|
[1360] | 112 | ! Kind definition added to all floating point numbers. |
---|
| 113 | ! lpm_init changed form a subroutine to a module. |
---|
[2312] | 114 | ! |
---|
[1329] | 115 | ! 1327 2014-03-21 11:00:16Z raasch |
---|
| 116 | ! -netcdf_output |
---|
| 117 | ! |
---|
[1323] | 118 | ! 1322 2014-03-20 16:38:49Z raasch |
---|
| 119 | ! REAL functions provided with KIND-attribute |
---|
| 120 | ! |
---|
[1321] | 121 | ! 1320 2014-03-20 08:40:49Z raasch |
---|
[1320] | 122 | ! ONLY-attribute added to USE-statements, |
---|
| 123 | ! kind-parameters added to all INTEGER and REAL declaration statements, |
---|
| 124 | ! kinds are defined in new module kinds, |
---|
| 125 | ! revision history before 2012 removed, |
---|
| 126 | ! comment fields (!:) to be used for variable explanations added to |
---|
| 127 | ! all variable declaration statements |
---|
| 128 | ! bugfix: #if defined( __parallel ) added |
---|
[850] | 129 | ! |
---|
[1315] | 130 | ! 1314 2014-03-14 18:25:17Z suehring |
---|
[2312] | 131 | ! Vertical logarithmic interpolation of horizontal particle speed for particles |
---|
[1315] | 132 | ! between roughness height and first vertical grid level. |
---|
| 133 | ! |
---|
[1093] | 134 | ! 1092 2013-02-02 11:24:22Z raasch |
---|
| 135 | ! unused variables removed |
---|
| 136 | ! |
---|
[1037] | 137 | ! 1036 2012-10-22 13:43:42Z raasch |
---|
| 138 | ! code put under GPL (PALM 3.9) |
---|
| 139 | ! |
---|
[850] | 140 | ! 849 2012-03-15 10:35:09Z raasch |
---|
[849] | 141 | ! routine renamed: init_particles -> lpm_init |
---|
| 142 | ! de_dx, de_dy, de_dz are allocated here (instead of automatic arrays in |
---|
| 143 | ! advec_particles), |
---|
| 144 | ! sort_particles renamed lpm_sort_arrays, user_init_particles renamed lpm_init |
---|
[392] | 145 | ! |
---|
[829] | 146 | ! 828 2012-02-21 12:00:36Z raasch |
---|
| 147 | ! call of init_kernels, particle feature color renamed class |
---|
| 148 | ! |
---|
[826] | 149 | ! 824 2012-02-17 09:09:57Z raasch |
---|
| 150 | ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3, |
---|
| 151 | ! array particles implemented as pointer |
---|
| 152 | ! |
---|
[668] | 153 | ! 667 2010-12-23 12:06:00Z suehring/gryschka |
---|
[2312] | 154 | ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation |
---|
[668] | 155 | ! of arrays. |
---|
| 156 | ! |
---|
[1] | 157 | ! Revision 1.1 1999/11/25 16:22:38 raasch |
---|
| 158 | ! Initial revision |
---|
| 159 | ! |
---|
| 160 | ! |
---|
| 161 | ! Description: |
---|
| 162 | ! ------------ |
---|
[1682] | 163 | !> This routine initializes a set of particles and their attributes (position, |
---|
| 164 | !> radius, ..) which are used by the Lagrangian particle model (see lpm). |
---|
[1] | 165 | !------------------------------------------------------------------------------! |
---|
[1682] | 166 | MODULE lpm_init_mod |
---|
[2312] | 167 | |
---|
[2305] | 168 | USE, INTRINSIC :: ISO_C_BINDING |
---|
[1] | 169 | |
---|
[1320] | 170 | USE arrays_3d, & |
---|
[2232] | 171 | ONLY: de_dx, de_dy, de_dz, zu, zw |
---|
[1320] | 172 | |
---|
| 173 | USE control_parameters, & |
---|
[1691] | 174 | ONLY: cloud_droplets, constant_flux_layer, current_timestep_number, & |
---|
[2312] | 175 | dt_3d, dz, initializing_actions, message_string, ocean, & |
---|
| 176 | simulated_time |
---|
[1320] | 177 | |
---|
| 178 | USE grid_variables, & |
---|
[1359] | 179 | ONLY: ddx, dx, ddy, dy |
---|
[1320] | 180 | |
---|
| 181 | USE indices, & |
---|
[1575] | 182 | ONLY: nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, & |
---|
[2317] | 183 | nzt, wall_flags_0 |
---|
[1320] | 184 | |
---|
| 185 | USE kinds |
---|
| 186 | |
---|
| 187 | USE lpm_collision_kernels_mod, & |
---|
| 188 | ONLY: init_kernels |
---|
| 189 | |
---|
[1783] | 190 | USE netcdf_interface, & |
---|
| 191 | ONLY: netcdf_data_format |
---|
| 192 | |
---|
[1320] | 193 | USE particle_attributes, & |
---|
[1359] | 194 | ONLY: alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & |
---|
| 195 | block_offset, block_offset_def, collision_kernel, & |
---|
[2265] | 196 | curvature_solution_effects, density_ratio, grid_particles, & |
---|
| 197 | isf,i_splitting_mode, initial_weighting_factor, ibc_par_b, & |
---|
| 198 | ibc_par_lr, ibc_par_ns, ibc_par_t, iran_part, log_z_z0, & |
---|
| 199 | max_number_of_particle_groups, min_nr_particle, & |
---|
[2305] | 200 | number_concentration, & |
---|
[2265] | 201 | number_particles_per_gridbox, number_of_particles, & |
---|
[1320] | 202 | number_of_particle_groups, number_of_sublayers, & |
---|
[1822] | 203 | offset_ocean_nzt, offset_ocean_nzt_m1, & |
---|
[1359] | 204 | particles, particle_advection_start, particle_groups, & |
---|
| 205 | particle_groups_type, particles_per_point, & |
---|
[2312] | 206 | particle_type, pdx, pdy, pdz, prt_count, psb, psl, psn, psr, & |
---|
| 207 | pss, pst, radius, random_start_position, & |
---|
[2265] | 208 | read_particles_from_restartfile, seed_follows_topography, & |
---|
| 209 | sgs_wf_part, sort_count, splitting_function, splitting_mode, & |
---|
| 210 | total_number_of_particles, use_sgs_for_particles, & |
---|
| 211 | write_particle_statistics, zero_particle, z0_av_global |
---|
[1320] | 212 | |
---|
[1] | 213 | USE pegrid |
---|
| 214 | |
---|
[1320] | 215 | USE random_function_mod, & |
---|
| 216 | ONLY: random_function |
---|
[1] | 217 | |
---|
[2232] | 218 | USE surface_mod, & |
---|
[2317] | 219 | ONLY: get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h |
---|
[2232] | 220 | |
---|
[1359] | 221 | IMPLICIT NONE |
---|
[1320] | 222 | |
---|
[1359] | 223 | PRIVATE |
---|
| 224 | |
---|
[1682] | 225 | INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !< |
---|
| 226 | INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !< |
---|
[1359] | 227 | |
---|
| 228 | INTERFACE lpm_init |
---|
| 229 | MODULE PROCEDURE lpm_init |
---|
| 230 | END INTERFACE lpm_init |
---|
| 231 | |
---|
| 232 | INTERFACE lpm_create_particle |
---|
| 233 | MODULE PROCEDURE lpm_create_particle |
---|
| 234 | END INTERFACE lpm_create_particle |
---|
| 235 | |
---|
| 236 | PUBLIC lpm_init, lpm_create_particle |
---|
| 237 | |
---|
[1929] | 238 | CONTAINS |
---|
[1359] | 239 | |
---|
[1682] | 240 | !------------------------------------------------------------------------------! |
---|
| 241 | ! Description: |
---|
| 242 | ! ------------ |
---|
| 243 | !> @todo Missing subroutine description. |
---|
| 244 | !------------------------------------------------------------------------------! |
---|
[1359] | 245 | SUBROUTINE lpm_init |
---|
| 246 | |
---|
| 247 | USE lpm_collision_kernels_mod, & |
---|
| 248 | ONLY: init_kernels |
---|
| 249 | |
---|
[1] | 250 | IMPLICIT NONE |
---|
| 251 | |
---|
[1682] | 252 | INTEGER(iwp) :: i !< |
---|
| 253 | INTEGER(iwp) :: j !< |
---|
| 254 | INTEGER(iwp) :: k !< |
---|
[1320] | 255 | |
---|
[2312] | 256 | REAL(wp) :: div !< |
---|
[1682] | 257 | REAL(wp) :: height_int !< |
---|
| 258 | REAL(wp) :: height_p !< |
---|
| 259 | REAL(wp) :: z_p !< |
---|
| 260 | REAL(wp) :: z0_av_local !< |
---|
[1] | 261 | |
---|
[1359] | 262 | |
---|
[1] | 263 | ! |
---|
[150] | 264 | !-- In case of oceans runs, the vertical index calculations need an offset, |
---|
| 265 | !-- because otherwise the k indices will become negative |
---|
| 266 | IF ( ocean ) THEN |
---|
| 267 | offset_ocean_nzt = nzt |
---|
| 268 | offset_ocean_nzt_m1 = nzt - 1 |
---|
| 269 | ENDIF |
---|
| 270 | |
---|
[1359] | 271 | ! |
---|
| 272 | !-- Define block offsets for dividing a gridcell in 8 sub cells |
---|
[150] | 273 | |
---|
[1359] | 274 | block_offset(0) = block_offset_def (-1,-1,-1) |
---|
| 275 | block_offset(1) = block_offset_def (-1,-1, 0) |
---|
| 276 | block_offset(2) = block_offset_def (-1, 0,-1) |
---|
| 277 | block_offset(3) = block_offset_def (-1, 0, 0) |
---|
| 278 | block_offset(4) = block_offset_def ( 0,-1,-1) |
---|
| 279 | block_offset(5) = block_offset_def ( 0,-1, 0) |
---|
| 280 | block_offset(6) = block_offset_def ( 0, 0,-1) |
---|
| 281 | block_offset(7) = block_offset_def ( 0, 0, 0) |
---|
[150] | 282 | ! |
---|
[1] | 283 | !-- Check the number of particle groups. |
---|
| 284 | IF ( number_of_particle_groups > max_number_of_particle_groups ) THEN |
---|
[274] | 285 | WRITE( message_string, * ) 'max_number_of_particle_groups =', & |
---|
| 286 | max_number_of_particle_groups , & |
---|
[254] | 287 | '&number_of_particle_groups reset to ', & |
---|
| 288 | max_number_of_particle_groups |
---|
[849] | 289 | CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 ) |
---|
[1] | 290 | number_of_particle_groups = max_number_of_particle_groups |
---|
| 291 | ENDIF |
---|
[2232] | 292 | ! |
---|
| 293 | !-- Check if downward-facing walls exist. This case, reflection boundary |
---|
| 294 | !-- conditions (as well as subgrid-scale velocities) may do not work |
---|
[2312] | 295 | !-- propably (not realized so far). |
---|
[2232] | 296 | IF ( surf_def_h(1)%ns >= 1 ) THEN |
---|
[2312] | 297 | WRITE( message_string, * ) 'Overhanging topography do not work '// & |
---|
| 298 | 'with particles' |
---|
[2232] | 299 | CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 ) |
---|
[1] | 300 | |
---|
[2312] | 301 | ENDIF |
---|
[2232] | 302 | |
---|
[1] | 303 | ! |
---|
| 304 | !-- Set default start positions, if necessary |
---|
[1359] | 305 | IF ( psl(1) == 9999999.9_wp ) psl(1) = -0.5_wp * dx |
---|
| 306 | IF ( psr(1) == 9999999.9_wp ) psr(1) = ( nx + 0.5_wp ) * dx |
---|
| 307 | IF ( pss(1) == 9999999.9_wp ) pss(1) = -0.5_wp * dy |
---|
| 308 | IF ( psn(1) == 9999999.9_wp ) psn(1) = ( ny + 0.5_wp ) * dy |
---|
| 309 | IF ( psb(1) == 9999999.9_wp ) psb(1) = zu(nz/2) |
---|
| 310 | IF ( pst(1) == 9999999.9_wp ) pst(1) = psb(1) |
---|
[1] | 311 | |
---|
[1359] | 312 | IF ( pdx(1) == 9999999.9_wp .OR. pdx(1) == 0.0_wp ) pdx(1) = dx |
---|
| 313 | IF ( pdy(1) == 9999999.9_wp .OR. pdy(1) == 0.0_wp ) pdy(1) = dy |
---|
| 314 | IF ( pdz(1) == 9999999.9_wp .OR. pdz(1) == 0.0_wp ) pdz(1) = zu(2) - zu(1) |
---|
[1] | 315 | |
---|
[2182] | 316 | ! |
---|
| 317 | !-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are |
---|
| 318 | !-- calculated diagnostically. Therfore an isotropic distribution is prescribed. |
---|
[2312] | 319 | IF ( number_particles_per_gridbox /= -1 .AND. & |
---|
[2182] | 320 | number_particles_per_gridbox >= 1 ) THEN |
---|
[2312] | 321 | pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) / & |
---|
[2182] | 322 | REAL(number_particles_per_gridbox))**0.3333333_wp |
---|
| 323 | ! |
---|
[2312] | 324 | !-- Ensure a smooth value (two significant digits) of distance between |
---|
[2182] | 325 | !-- particles (pdx, pdy, pdz). |
---|
| 326 | div = 1000.0_wp |
---|
| 327 | DO WHILE ( pdx(1) < div ) |
---|
| 328 | div = div / 10.0_wp |
---|
| 329 | ENDDO |
---|
| 330 | pdx(1) = NINT( pdx(1) * 100.0_wp / div ) * div / 100.0_wp |
---|
| 331 | pdy(1) = pdx(1) |
---|
| 332 | pdz(1) = pdx(1) |
---|
| 333 | |
---|
| 334 | ENDIF |
---|
| 335 | |
---|
[1] | 336 | DO j = 2, number_of_particle_groups |
---|
[1359] | 337 | IF ( psl(j) == 9999999.9_wp ) psl(j) = psl(j-1) |
---|
| 338 | IF ( psr(j) == 9999999.9_wp ) psr(j) = psr(j-1) |
---|
| 339 | IF ( pss(j) == 9999999.9_wp ) pss(j) = pss(j-1) |
---|
| 340 | IF ( psn(j) == 9999999.9_wp ) psn(j) = psn(j-1) |
---|
| 341 | IF ( psb(j) == 9999999.9_wp ) psb(j) = psb(j-1) |
---|
| 342 | IF ( pst(j) == 9999999.9_wp ) pst(j) = pst(j-1) |
---|
| 343 | IF ( pdx(j) == 9999999.9_wp .OR. pdx(j) == 0.0_wp ) pdx(j) = pdx(j-1) |
---|
| 344 | IF ( pdy(j) == 9999999.9_wp .OR. pdy(j) == 0.0_wp ) pdy(j) = pdy(j-1) |
---|
| 345 | IF ( pdz(j) == 9999999.9_wp .OR. pdz(j) == 0.0_wp ) pdz(j) = pdz(j-1) |
---|
[1] | 346 | ENDDO |
---|
| 347 | |
---|
| 348 | ! |
---|
[2312] | 349 | !-- Allocate arrays required for calculating particle SGS velocities. |
---|
[1929] | 350 | !-- Initialize prefactor required for stoachastic Weil equation. |
---|
[1822] | 351 | IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN |
---|
[849] | 352 | ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
| 353 | de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
| 354 | de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1929] | 355 | |
---|
[2312] | 356 | sgs_wf_part = 1.0_wp / 3.0_wp |
---|
[849] | 357 | ENDIF |
---|
| 358 | |
---|
| 359 | ! |
---|
[2312] | 360 | !-- Allocate array required for logarithmic vertical interpolation of |
---|
[1314] | 361 | !-- horizontal particle velocities between the surface and the first vertical |
---|
[2312] | 362 | !-- grid level. In order to avoid repeated CPU cost-intensive CALLS of |
---|
| 363 | !-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for |
---|
| 364 | !-- several heights. Splitting into 20 sublayers turned out to be sufficient. |
---|
[1314] | 365 | !-- To obtain exact height levels of particles, linear interpolation is applied |
---|
| 366 | !-- (see lpm_advec.f90). |
---|
[1691] | 367 | IF ( constant_flux_layer ) THEN |
---|
[2312] | 368 | |
---|
| 369 | ALLOCATE ( log_z_z0(0:number_of_sublayers) ) |
---|
[2232] | 370 | z_p = zu(nzb+1) - zw(nzb) |
---|
[1314] | 371 | |
---|
| 372 | ! |
---|
| 373 | !-- Calculate horizontal mean value of z0 used for logartihmic |
---|
[2312] | 374 | !-- interpolation. Note: this is not exact for heterogeneous z0. |
---|
| 375 | !-- However, sensitivity studies showed that the effect is |
---|
| 376 | !-- negligible. |
---|
[2232] | 377 | z0_av_local = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) + & |
---|
| 378 | SUM( surf_usm_h%z0 ) |
---|
[1359] | 379 | z0_av_global = 0.0_wp |
---|
[1314] | 380 | |
---|
[1320] | 381 | #if defined( __parallel ) |
---|
[1314] | 382 | CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, & |
---|
| 383 | comm2d, ierr ) |
---|
[1320] | 384 | #else |
---|
| 385 | z0_av_global = z0_av_local |
---|
| 386 | #endif |
---|
[1314] | 387 | |
---|
| 388 | z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) ) |
---|
| 389 | ! |
---|
| 390 | !-- Horizontal wind speed is zero below and at z0 |
---|
[1359] | 391 | log_z_z0(0) = 0.0_wp |
---|
[1314] | 392 | ! |
---|
| 393 | !-- Calculate vertical depth of the sublayers |
---|
[1322] | 394 | height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp ) |
---|
[1314] | 395 | ! |
---|
| 396 | !-- Precalculate LOG(z/z0) |
---|
[1929] | 397 | height_p = z0_av_global |
---|
[1314] | 398 | DO k = 1, number_of_sublayers |
---|
| 399 | |
---|
| 400 | height_p = height_p + height_int |
---|
| 401 | log_z_z0(k) = LOG( height_p / z0_av_global ) |
---|
| 402 | |
---|
| 403 | ENDDO |
---|
| 404 | |
---|
| 405 | ENDIF |
---|
| 406 | |
---|
| 407 | ! |
---|
[1359] | 408 | !-- Check boundary condition and set internal variables |
---|
| 409 | SELECT CASE ( bc_par_b ) |
---|
[2312] | 410 | |
---|
[1359] | 411 | CASE ( 'absorb' ) |
---|
| 412 | ibc_par_b = 1 |
---|
| 413 | |
---|
| 414 | CASE ( 'reflect' ) |
---|
| 415 | ibc_par_b = 2 |
---|
[2312] | 416 | |
---|
[1359] | 417 | CASE DEFAULT |
---|
| 418 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 419 | 'bc_par_b = "', TRIM( bc_par_b ), '"' |
---|
| 420 | CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 ) |
---|
[2312] | 421 | |
---|
[1359] | 422 | END SELECT |
---|
| 423 | SELECT CASE ( bc_par_t ) |
---|
[2312] | 424 | |
---|
[1359] | 425 | CASE ( 'absorb' ) |
---|
| 426 | ibc_par_t = 1 |
---|
| 427 | |
---|
| 428 | CASE ( 'reflect' ) |
---|
| 429 | ibc_par_t = 2 |
---|
[2312] | 430 | |
---|
[1359] | 431 | CASE DEFAULT |
---|
| 432 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 433 | 'bc_par_t = "', TRIM( bc_par_t ), '"' |
---|
| 434 | CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 ) |
---|
[2312] | 435 | |
---|
[1359] | 436 | END SELECT |
---|
| 437 | SELECT CASE ( bc_par_lr ) |
---|
| 438 | |
---|
| 439 | CASE ( 'cyclic' ) |
---|
| 440 | ibc_par_lr = 0 |
---|
| 441 | |
---|
| 442 | CASE ( 'absorb' ) |
---|
| 443 | ibc_par_lr = 1 |
---|
| 444 | |
---|
| 445 | CASE ( 'reflect' ) |
---|
| 446 | ibc_par_lr = 2 |
---|
[2312] | 447 | |
---|
[1359] | 448 | CASE DEFAULT |
---|
| 449 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 450 | 'bc_par_lr = "', TRIM( bc_par_lr ), '"' |
---|
| 451 | CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 ) |
---|
[2312] | 452 | |
---|
[1359] | 453 | END SELECT |
---|
| 454 | SELECT CASE ( bc_par_ns ) |
---|
| 455 | |
---|
| 456 | CASE ( 'cyclic' ) |
---|
| 457 | ibc_par_ns = 0 |
---|
| 458 | |
---|
| 459 | CASE ( 'absorb' ) |
---|
| 460 | ibc_par_ns = 1 |
---|
| 461 | |
---|
| 462 | CASE ( 'reflect' ) |
---|
| 463 | ibc_par_ns = 2 |
---|
[2312] | 464 | |
---|
[1359] | 465 | CASE DEFAULT |
---|
| 466 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 467 | 'bc_par_ns = "', TRIM( bc_par_ns ), '"' |
---|
| 468 | CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 ) |
---|
[2312] | 469 | |
---|
[1359] | 470 | END SELECT |
---|
[2263] | 471 | SELECT CASE ( splitting_mode ) |
---|
[2312] | 472 | |
---|
[2263] | 473 | CASE ( 'const' ) |
---|
| 474 | i_splitting_mode = 1 |
---|
[1359] | 475 | |
---|
[2263] | 476 | CASE ( 'cl_av' ) |
---|
| 477 | i_splitting_mode = 2 |
---|
| 478 | |
---|
| 479 | CASE ( 'gb_av' ) |
---|
| 480 | i_splitting_mode = 3 |
---|
[2312] | 481 | |
---|
[2263] | 482 | CASE DEFAULT |
---|
| 483 | WRITE( message_string, * ) 'unknown splitting condition ', & |
---|
| 484 | 'splitting_mode = "', TRIM( splitting_mode ), '"' |
---|
| 485 | CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 ) |
---|
[2312] | 486 | |
---|
[2263] | 487 | END SELECT |
---|
| 488 | SELECT CASE ( splitting_function ) |
---|
[2312] | 489 | |
---|
[2263] | 490 | CASE ( 'gamma' ) |
---|
| 491 | isf = 1 |
---|
| 492 | |
---|
| 493 | CASE ( 'log' ) |
---|
| 494 | isf = 2 |
---|
| 495 | |
---|
| 496 | CASE ( 'exp' ) |
---|
| 497 | isf = 3 |
---|
[2312] | 498 | |
---|
[2263] | 499 | CASE DEFAULT |
---|
| 500 | WRITE( message_string, * ) 'unknown splitting function ', & |
---|
| 501 | 'splitting_function = "', TRIM( splitting_function ), '"' |
---|
| 502 | CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 ) |
---|
[2312] | 503 | |
---|
[2263] | 504 | END SELECT |
---|
| 505 | |
---|
[2312] | 506 | |
---|
[1359] | 507 | ! |
---|
[828] | 508 | !-- Initialize collision kernels |
---|
| 509 | IF ( collision_kernel /= 'none' ) CALL init_kernels |
---|
| 510 | |
---|
| 511 | ! |
---|
[1] | 512 | !-- For the first model run of a possible job chain initialize the |
---|
[849] | 513 | !-- particles, otherwise read the particle data from restart file. |
---|
[1] | 514 | IF ( TRIM( initializing_actions ) == 'read_restart_data' & |
---|
| 515 | .AND. read_particles_from_restartfile ) THEN |
---|
| 516 | |
---|
[849] | 517 | CALL lpm_read_restart_file |
---|
[1] | 518 | |
---|
| 519 | ELSE |
---|
| 520 | |
---|
| 521 | ! |
---|
| 522 | !-- Allocate particle arrays and set attributes of the initial set of |
---|
| 523 | !-- particles, which can be also periodically released at later times. |
---|
[1359] | 524 | ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
| 525 | grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) ) |
---|
[1] | 526 | |
---|
[1359] | 527 | number_of_particles = 0 |
---|
[792] | 528 | |
---|
| 529 | sort_count = 0 |
---|
[1359] | 530 | prt_count = 0 |
---|
[792] | 531 | |
---|
[1] | 532 | ! |
---|
[2312] | 533 | !-- initialize counter for particle IDs |
---|
[2305] | 534 | grid_particles%id_counter = 1 |
---|
[2122] | 535 | |
---|
| 536 | ! |
---|
[1] | 537 | !-- Initialize all particles with dummy values (otherwise errors may |
---|
| 538 | !-- occur within restart runs). The reason for this is still not clear |
---|
| 539 | !-- and may be presumably caused by errors in the respective user-interface. |
---|
[1359] | 540 | zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & |
---|
| 541 | 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & |
---|
| 542 | 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & |
---|
[2312] | 543 | 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & |
---|
[2305] | 544 | 0, 0, 0_idp, .FALSE., -1 ) |
---|
[1822] | 545 | |
---|
[1359] | 546 | particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) |
---|
[1] | 547 | |
---|
| 548 | ! |
---|
| 549 | !-- Set values for the density ratio and radius for all particle |
---|
| 550 | !-- groups, if necessary |
---|
[1359] | 551 | IF ( density_ratio(1) == 9999999.9_wp ) density_ratio(1) = 0.0_wp |
---|
| 552 | IF ( radius(1) == 9999999.9_wp ) radius(1) = 0.0_wp |
---|
[1] | 553 | DO i = 2, number_of_particle_groups |
---|
[1359] | 554 | IF ( density_ratio(i) == 9999999.9_wp ) THEN |
---|
[1] | 555 | density_ratio(i) = density_ratio(i-1) |
---|
| 556 | ENDIF |
---|
[1359] | 557 | IF ( radius(i) == 9999999.9_wp ) radius(i) = radius(i-1) |
---|
[1] | 558 | ENDDO |
---|
| 559 | |
---|
| 560 | DO i = 1, number_of_particle_groups |
---|
[1359] | 561 | IF ( density_ratio(i) /= 0.0_wp .AND. radius(i) == 0 ) THEN |
---|
[254] | 562 | WRITE( message_string, * ) 'particle group #', i, 'has a', & |
---|
| 563 | 'density ratio /= 0 but radius = 0' |
---|
[849] | 564 | CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 ) |
---|
[1] | 565 | ENDIF |
---|
| 566 | particle_groups(i)%density_ratio = density_ratio(i) |
---|
| 567 | particle_groups(i)%radius = radius(i) |
---|
| 568 | ENDDO |
---|
| 569 | |
---|
| 570 | ! |
---|
[1359] | 571 | !-- Set a seed value for the random number generator to be exclusively |
---|
| 572 | !-- used for the particle code. The generated random numbers should be |
---|
| 573 | !-- different on the different PEs. |
---|
| 574 | iran_part = iran_part + myid |
---|
| 575 | |
---|
[1725] | 576 | CALL lpm_create_particle (PHASE_INIT) |
---|
[1359] | 577 | ! |
---|
| 578 | !-- User modification of initial particles |
---|
| 579 | CALL user_lpm_init |
---|
| 580 | |
---|
| 581 | ! |
---|
| 582 | !-- Open file for statistical informations about particle conditions |
---|
| 583 | IF ( write_particle_statistics ) THEN |
---|
| 584 | CALL check_open( 80 ) |
---|
| 585 | WRITE ( 80, 8000 ) current_timestep_number, simulated_time, & |
---|
[2312] | 586 | number_of_particles |
---|
[1359] | 587 | CALL close_file( 80 ) |
---|
| 588 | ENDIF |
---|
| 589 | |
---|
| 590 | ENDIF |
---|
| 591 | |
---|
| 592 | ! |
---|
[2312] | 593 | !-- To avoid programm abort, assign particles array to the local version of |
---|
[1359] | 594 | !-- first grid cell |
---|
| 595 | number_of_particles = prt_count(nzb+1,nys,nxl) |
---|
| 596 | particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles) |
---|
| 597 | ! |
---|
| 598 | !-- Formats |
---|
| 599 | 8000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10) |
---|
| 600 | |
---|
| 601 | END SUBROUTINE lpm_init |
---|
| 602 | |
---|
[1682] | 603 | !------------------------------------------------------------------------------! |
---|
| 604 | ! Description: |
---|
| 605 | ! ------------ |
---|
| 606 | !> @todo Missing subroutine description. |
---|
| 607 | !------------------------------------------------------------------------------! |
---|
[1359] | 608 | SUBROUTINE lpm_create_particle (phase) |
---|
| 609 | |
---|
| 610 | USE lpm_exchange_horiz_mod, & |
---|
| 611 | ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array |
---|
| 612 | |
---|
| 613 | USE lpm_pack_arrays_mod, & |
---|
| 614 | ONLY: lpm_pack_all_arrays |
---|
| 615 | |
---|
[1871] | 616 | USE particle_attributes, & |
---|
[2312] | 617 | ONLY: deleted_particles |
---|
[1871] | 618 | |
---|
[1359] | 619 | IMPLICIT NONE |
---|
| 620 | |
---|
[1929] | 621 | INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for particles |
---|
| 622 | INTEGER(iwp) :: i !< loop variable ( particle groups ) |
---|
| 623 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 624 | INTEGER(iwp) :: j !< loop variable ( particles per point ) |
---|
| 625 | INTEGER(iwp) :: jp !< index variable along y |
---|
[2232] | 626 | INTEGER(iwp) :: k !< index variable along z |
---|
| 627 | INTEGER(iwp) :: k_surf !< index of surface grid point |
---|
[1929] | 628 | INTEGER(iwp) :: kp !< index variable along z |
---|
| 629 | INTEGER(iwp) :: loop_stride !< loop variable for initialization |
---|
| 630 | INTEGER(iwp) :: n !< loop variable ( number of particles ) |
---|
| 631 | INTEGER(iwp) :: new_size !< new size of allocated memory for particles |
---|
[1359] | 632 | |
---|
[1929] | 633 | INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization |
---|
[1359] | 634 | |
---|
[1929] | 635 | INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< start address of new particle |
---|
| 636 | INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< start address of new particle |
---|
[1359] | 637 | |
---|
[1929] | 638 | LOGICAL :: first_stride !< flag for initialization |
---|
[1359] | 639 | |
---|
[2312] | 640 | REAL(wp) :: pos_x !< increment for particle position in x |
---|
| 641 | REAL(wp) :: pos_y !< increment for particle position in y |
---|
| 642 | REAL(wp) :: pos_z !< increment for particle position in z |
---|
[1929] | 643 | REAL(wp) :: rand_contr !< dummy argument for random position |
---|
[1359] | 644 | |
---|
[1929] | 645 | TYPE(particle_type),TARGET :: tmp_particle !< temporary particle used for initialization |
---|
[1359] | 646 | |
---|
| 647 | ! |
---|
| 648 | !-- Calculate particle positions and store particle attributes, if |
---|
| 649 | !-- particle is situated on this PE |
---|
| 650 | DO loop_stride = 1, 2 |
---|
| 651 | first_stride = (loop_stride == 1) |
---|
| 652 | IF ( first_stride ) THEN |
---|
| 653 | local_count = 0 ! count number of particles |
---|
| 654 | ELSE |
---|
| 655 | local_count = prt_count ! Start address of new particles |
---|
| 656 | ENDIF |
---|
| 657 | |
---|
[2182] | 658 | ! |
---|
| 659 | !-- Calculate initial_weighting_factor diagnostically |
---|
| 660 | IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN |
---|
[2312] | 661 | initial_weighting_factor = number_concentration * 1.0E6_wp * & |
---|
| 662 | pdx(1) * pdy(1) * pdz(1) |
---|
[2182] | 663 | END IF |
---|
| 664 | |
---|
[1] | 665 | n = 0 |
---|
| 666 | DO i = 1, number_of_particle_groups |
---|
| 667 | |
---|
| 668 | pos_z = psb(i) |
---|
| 669 | |
---|
| 670 | DO WHILE ( pos_z <= pst(i) ) |
---|
| 671 | |
---|
[2223] | 672 | IF ( pos_z >= 0.0_wp .AND. pos_z < zw(nzt) ) THEN |
---|
[1] | 673 | |
---|
| 674 | |
---|
[2223] | 675 | pos_y = pss(i) |
---|
[1] | 676 | |
---|
[2223] | 677 | DO WHILE ( pos_y <= psn(i) ) |
---|
[1] | 678 | |
---|
[2223] | 679 | IF ( pos_y >= ( nys - 0.5_wp ) * dy .AND. & |
---|
| 680 | pos_y < ( nyn + 0.5_wp ) * dy ) THEN |
---|
[1] | 681 | |
---|
[2223] | 682 | pos_x = psl(i) |
---|
[1] | 683 | |
---|
[2223] | 684 | xloop: DO WHILE ( pos_x <= psr(i) ) |
---|
[1] | 685 | |
---|
[2312] | 686 | IF ( pos_x >= ( nxl - 0.5_wp ) * dx .AND. & |
---|
[2223] | 687 | pos_x < ( nxr + 0.5_wp ) * dx ) THEN |
---|
| 688 | |
---|
| 689 | DO j = 1, particles_per_point |
---|
| 690 | |
---|
[2305] | 691 | |
---|
[2223] | 692 | n = n + 1 |
---|
| 693 | tmp_particle%x = pos_x |
---|
| 694 | tmp_particle%y = pos_y |
---|
| 695 | tmp_particle%z = pos_z |
---|
| 696 | tmp_particle%age = 0.0_wp |
---|
| 697 | tmp_particle%age_m = 0.0_wp |
---|
| 698 | tmp_particle%dt_sum = 0.0_wp |
---|
| 699 | tmp_particle%e_m = 0.0_wp |
---|
[2312] | 700 | tmp_particle%rvar1 = 0.0_wp |
---|
| 701 | tmp_particle%rvar2 = 0.0_wp |
---|
| 702 | tmp_particle%rvar3 = 0.0_wp |
---|
[2223] | 703 | tmp_particle%speed_x = 0.0_wp |
---|
| 704 | tmp_particle%speed_y = 0.0_wp |
---|
| 705 | tmp_particle%speed_z = 0.0_wp |
---|
| 706 | tmp_particle%origin_x = pos_x |
---|
| 707 | tmp_particle%origin_y = pos_y |
---|
| 708 | tmp_particle%origin_z = pos_z |
---|
[2312] | 709 | IF ( curvature_solution_effects ) THEN |
---|
| 710 | tmp_particle%aux1 = 0.0_wp ! dry aerosol radius |
---|
| 711 | tmp_particle%aux2 = dt_3d ! last Rosenbrock timestep |
---|
| 712 | ELSE |
---|
| 713 | tmp_particle%aux1 = 0.0_wp ! free to use |
---|
| 714 | tmp_particle%aux2 = 0.0_wp ! free to use |
---|
| 715 | ENDIF |
---|
[2223] | 716 | tmp_particle%radius = particle_groups(i)%radius |
---|
| 717 | tmp_particle%weight_factor = initial_weighting_factor |
---|
| 718 | tmp_particle%class = 1 |
---|
| 719 | tmp_particle%group = i |
---|
[2305] | 720 | tmp_particle%id = 0_idp |
---|
[2223] | 721 | tmp_particle%particle_mask = .TRUE. |
---|
| 722 | tmp_particle%block_nr = -1 |
---|
[1] | 723 | ! |
---|
[2223] | 724 | !-- Determine the grid indices of the particle position |
---|
[2232] | 725 | ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx |
---|
[2223] | 726 | jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy |
---|
| 727 | kp = tmp_particle%z / dz + 1 + offset_ocean_nzt |
---|
[2232] | 728 | ! |
---|
| 729 | !-- Determine surface level. Therefore, check for |
---|
[2317] | 730 | !-- upward-facing wall on w-grid. |
---|
| 731 | k_surf = get_topography_top_index( j, i, 'w' ) |
---|
[1] | 732 | |
---|
[2223] | 733 | IF ( seed_follows_topography ) THEN |
---|
[1575] | 734 | ! |
---|
[2223] | 735 | !-- Particle height is given relative to topography |
---|
[2232] | 736 | kp = kp + k_surf |
---|
| 737 | tmp_particle%z = tmp_particle%z + zw(k_surf) |
---|
| 738 | !-- Skip particle release if particle position is |
---|
| 739 | !-- above model top, or within topography in case |
---|
| 740 | !-- of overhanging structures. |
---|
| 741 | IF ( kp > nzt .OR. & |
---|
| 742 | .NOT. BTEST( wall_flags_0(kp,jp,ip), 0 ) ) THEN |
---|
[2223] | 743 | pos_x = pos_x + pdx(i) |
---|
| 744 | CYCLE xloop |
---|
| 745 | ENDIF |
---|
[2232] | 746 | ! |
---|
[2312] | 747 | !-- Skip particle release if particle position is |
---|
[2232] | 748 | !-- below surface, or within topography in case |
---|
| 749 | !-- of overhanging structures. |
---|
[2223] | 750 | ELSEIF ( .NOT. seed_follows_topography .AND. & |
---|
[2232] | 751 | tmp_particle%z <= zw(k_surf) .OR. & |
---|
| 752 | .NOT. BTEST( wall_flags_0(kp,jp,ip), 0 ) )& |
---|
| 753 | THEN |
---|
[1575] | 754 | pos_x = pos_x + pdx(i) |
---|
[2312] | 755 | CYCLE xloop |
---|
[1575] | 756 | ENDIF |
---|
| 757 | |
---|
[2223] | 758 | local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1 |
---|
[2182] | 759 | |
---|
[2223] | 760 | IF ( .NOT. first_stride ) THEN |
---|
| 761 | IF ( ip < nxl .OR. jp < nys .OR. kp < nzb+1 ) THEN |
---|
| 762 | write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1 |
---|
| 763 | ENDIF |
---|
| 764 | IF ( ip > nxr .OR. jp > nyn .OR. kp > nzt ) THEN |
---|
| 765 | write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt |
---|
| 766 | ENDIF |
---|
| 767 | grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle |
---|
| 768 | |
---|
[1359] | 769 | ENDIF |
---|
[2223] | 770 | ENDDO |
---|
[1929] | 771 | |
---|
[2223] | 772 | ENDIF |
---|
[1] | 773 | |
---|
[2223] | 774 | pos_x = pos_x + pdx(i) |
---|
[1] | 775 | |
---|
[2223] | 776 | ENDDO xloop |
---|
[1] | 777 | |
---|
[2223] | 778 | ENDIF |
---|
[1] | 779 | |
---|
[2223] | 780 | pos_y = pos_y + pdy(i) |
---|
[1] | 781 | |
---|
[2223] | 782 | ENDDO |
---|
[1] | 783 | |
---|
[2223] | 784 | ENDIF |
---|
[1] | 785 | |
---|
| 786 | pos_z = pos_z + pdz(i) |
---|
| 787 | |
---|
| 788 | ENDDO |
---|
| 789 | |
---|
| 790 | ENDDO |
---|
| 791 | |
---|
[1359] | 792 | IF ( first_stride ) THEN |
---|
| 793 | DO ip = nxl, nxr |
---|
| 794 | DO jp = nys, nyn |
---|
| 795 | DO kp = nzb+1, nzt |
---|
| 796 | IF ( phase == PHASE_INIT ) THEN |
---|
| 797 | IF ( local_count(kp,jp,ip) > 0 ) THEN |
---|
| 798 | alloc_size = MAX( INT( local_count(kp,jp,ip) * & |
---|
| 799 | ( 1.0_wp + alloc_factor / 100.0_wp ) ), & |
---|
| 800 | min_nr_particle ) |
---|
| 801 | ELSE |
---|
| 802 | alloc_size = min_nr_particle |
---|
| 803 | ENDIF |
---|
| 804 | ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size)) |
---|
| 805 | DO n = 1, alloc_size |
---|
| 806 | grid_particles(kp,jp,ip)%particles(n) = zero_particle |
---|
| 807 | ENDDO |
---|
| 808 | ELSEIF ( phase == PHASE_RELEASE ) THEN |
---|
| 809 | IF ( local_count(kp,jp,ip) > 0 ) THEN |
---|
| 810 | new_size = local_count(kp,jp,ip) + prt_count(kp,jp,ip) |
---|
| 811 | alloc_size = MAX( INT( new_size * ( 1.0_wp + & |
---|
| 812 | alloc_factor / 100.0_wp ) ), min_nr_particle ) |
---|
| 813 | IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN |
---|
[2182] | 814 | CALL realloc_particles_array(ip,jp,kp,alloc_size) |
---|
[1359] | 815 | ENDIF |
---|
| 816 | ENDIF |
---|
| 817 | ENDIF |
---|
| 818 | ENDDO |
---|
| 819 | ENDDO |
---|
| 820 | ENDDO |
---|
| 821 | ENDIF |
---|
[1929] | 822 | |
---|
[1359] | 823 | ENDDO |
---|
[1] | 824 | |
---|
[2182] | 825 | |
---|
| 826 | |
---|
[1359] | 827 | local_start = prt_count+1 |
---|
| 828 | prt_count = local_count |
---|
[1871] | 829 | |
---|
[1] | 830 | ! |
---|
[2122] | 831 | !-- Calculate particle IDs |
---|
| 832 | DO ip = nxl, nxr |
---|
| 833 | DO jp = nys, nyn |
---|
| 834 | DO kp = nzb+1, nzt |
---|
| 835 | number_of_particles = prt_count(kp,jp,ip) |
---|
| 836 | IF ( number_of_particles <= 0 ) CYCLE |
---|
| 837 | particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) |
---|
| 838 | |
---|
[2312] | 839 | DO n = local_start(kp,jp,ip), number_of_particles !only new particles |
---|
[2122] | 840 | |
---|
[2305] | 841 | particles(n)%id = 10000_idp**3 * grid_particles(kp,jp,ip)%id_counter + & |
---|
| 842 | 10000_idp**2 * kp + 10000_idp * jp + ip |
---|
| 843 | ! |
---|
| 844 | !-- Count the number of particles that have been released before |
---|
[2122] | 845 | grid_particles(kp,jp,ip)%id_counter = & |
---|
| 846 | grid_particles(kp,jp,ip)%id_counter + 1 |
---|
| 847 | |
---|
| 848 | ENDDO |
---|
| 849 | |
---|
| 850 | ENDDO |
---|
| 851 | ENDDO |
---|
| 852 | ENDDO |
---|
| 853 | |
---|
| 854 | ! |
---|
[1871] | 855 | !-- Initialize aerosol background spectrum |
---|
[2312] | 856 | IF ( curvature_solution_effects ) THEN |
---|
[1871] | 857 | CALL lpm_init_aerosols(local_start) |
---|
| 858 | ENDIF |
---|
| 859 | |
---|
| 860 | ! |
---|
[1929] | 861 | !-- Add random fluctuation to particle positions. |
---|
[1359] | 862 | IF ( random_start_position ) THEN |
---|
| 863 | DO ip = nxl, nxr |
---|
| 864 | DO jp = nys, nyn |
---|
| 865 | DO kp = nzb+1, nzt |
---|
| 866 | number_of_particles = prt_count(kp,jp,ip) |
---|
| 867 | IF ( number_of_particles <= 0 ) CYCLE |
---|
| 868 | particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) |
---|
[1929] | 869 | ! |
---|
[2312] | 870 | !-- Move only new particles. Moreover, limit random fluctuation |
---|
| 871 | !-- in order to prevent that particles move more than one grid box, |
---|
| 872 | !-- which would lead to problems concerning particle exchange |
---|
| 873 | !-- between processors in case pdx/pdy are larger than dx/dy, |
---|
| 874 | !-- respectively. |
---|
[1929] | 875 | DO n = local_start(kp,jp,ip), number_of_particles |
---|
[1359] | 876 | IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN |
---|
[1929] | 877 | rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & |
---|
| 878 | pdx(particles(n)%group) |
---|
[1359] | 879 | particles(n)%x = particles(n)%x + & |
---|
[2232] | 880 | MERGE( rand_contr, SIGN( dx, rand_contr ), & |
---|
[1929] | 881 | ABS( rand_contr ) < dx & |
---|
[2312] | 882 | ) |
---|
[1359] | 883 | ENDIF |
---|
| 884 | IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN |
---|
[1929] | 885 | rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & |
---|
| 886 | pdy(particles(n)%group) |
---|
[1359] | 887 | particles(n)%y = particles(n)%y + & |
---|
[2232] | 888 | MERGE( rand_contr, SIGN( dy, rand_contr ), & |
---|
[1929] | 889 | ABS( rand_contr ) < dy & |
---|
[2312] | 890 | ) |
---|
[1359] | 891 | ENDIF |
---|
| 892 | IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN |
---|
[1929] | 893 | rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & |
---|
| 894 | pdz(particles(n)%group) |
---|
[1359] | 895 | particles(n)%z = particles(n)%z + & |
---|
[2232] | 896 | MERGE( rand_contr, SIGN( dz, rand_contr ), & |
---|
[1929] | 897 | ABS( rand_contr ) < dz & |
---|
[2312] | 898 | ) |
---|
[1359] | 899 | ENDIF |
---|
| 900 | ENDDO |
---|
[1] | 901 | ! |
---|
[1929] | 902 | !-- Identify particles located outside the model domain and reflect |
---|
| 903 | !-- or absorb them if necessary. |
---|
[1359] | 904 | CALL lpm_boundary_conds( 'bottom/top' ) |
---|
[1929] | 905 | ! |
---|
| 906 | !-- Furthermore, remove particles located in topography. Note, as |
---|
[2312] | 907 | !-- the particle speed is still zero at this point, wall |
---|
[1929] | 908 | !-- reflection boundary conditions will not work in this case. |
---|
| 909 | particles => & |
---|
| 910 | grid_particles(kp,jp,ip)%particles(1:number_of_particles) |
---|
| 911 | DO n = local_start(kp,jp,ip), number_of_particles |
---|
| 912 | i = ( particles(n)%x + 0.5_wp * dx ) * ddx |
---|
| 913 | j = ( particles(n)%y + 0.5_wp * dy ) * ddy |
---|
[2232] | 914 | k = particles(n)%z / dz + 1 + offset_ocean_nzt |
---|
| 915 | ! |
---|
| 916 | !-- Check if particle is within topography |
---|
| 917 | IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) ) THEN |
---|
[1929] | 918 | particles(n)%particle_mask = .FALSE. |
---|
| 919 | deleted_particles = deleted_particles + 1 |
---|
| 920 | ENDIF |
---|
[2232] | 921 | |
---|
[1929] | 922 | ENDDO |
---|
[1359] | 923 | ENDDO |
---|
| 924 | ENDDO |
---|
| 925 | ENDDO |
---|
[1] | 926 | ! |
---|
[1359] | 927 | !-- Exchange particles between grid cells and processors |
---|
| 928 | CALL lpm_move_particle |
---|
| 929 | CALL lpm_exchange_horiz |
---|
[1] | 930 | |
---|
[1359] | 931 | ENDIF |
---|
[1] | 932 | ! |
---|
[2312] | 933 | !-- In case of random_start_position, delete particles identified by |
---|
| 934 | !-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks, |
---|
| 935 | !-- which is needed for a fast interpolation of the LES fields on the particle |
---|
[1359] | 936 | !-- position. |
---|
| 937 | CALL lpm_pack_all_arrays |
---|
[1] | 938 | |
---|
| 939 | ! |
---|
[2265] | 940 | !-- Determine the current number of particles |
---|
[1359] | 941 | DO ip = nxl, nxr |
---|
| 942 | DO jp = nys, nyn |
---|
| 943 | DO kp = nzb+1, nzt |
---|
| 944 | number_of_particles = number_of_particles & |
---|
| 945 | + prt_count(kp,jp,ip) |
---|
[1] | 946 | ENDDO |
---|
[1359] | 947 | ENDDO |
---|
| 948 | ENDDO |
---|
[1] | 949 | ! |
---|
[1822] | 950 | !-- Calculate the number of particles of the total domain |
---|
[1] | 951 | #if defined( __parallel ) |
---|
[1359] | 952 | IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 953 | CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, & |
---|
| 954 | MPI_INTEGER, MPI_SUM, comm2d, ierr ) |
---|
[1] | 955 | #else |
---|
[1359] | 956 | total_number_of_particles = number_of_particles |
---|
[1] | 957 | #endif |
---|
| 958 | |
---|
[1359] | 959 | RETURN |
---|
[1] | 960 | |
---|
[1359] | 961 | END SUBROUTINE lpm_create_particle |
---|
[336] | 962 | |
---|
[1871] | 963 | SUBROUTINE lpm_init_aerosols(local_start) |
---|
| 964 | |
---|
| 965 | USE arrays_3d, & |
---|
[2312] | 966 | ONLY: hyp, pt, q |
---|
[1871] | 967 | |
---|
| 968 | USE cloud_parameters, & |
---|
[2122] | 969 | ONLY: l_d_rv, rho_l, r_v |
---|
[1871] | 970 | |
---|
| 971 | USE constants, & |
---|
| 972 | ONLY: pi |
---|
| 973 | |
---|
| 974 | USE kinds |
---|
| 975 | |
---|
| 976 | USE particle_attributes, & |
---|
[2312] | 977 | ONLY: aero_type, aero_weight, log_sigma, molecular_weight_of_solute, & |
---|
| 978 | molecular_weight_of_water, na, rho_s, rm, vanthoff |
---|
[1871] | 979 | |
---|
| 980 | IMPLICIT NONE |
---|
| 981 | |
---|
[2122] | 982 | REAL(wp) :: afactor !< curvature effects |
---|
[1871] | 983 | REAL(wp) :: bfactor !< solute effects |
---|
[2312] | 984 | REAL(wp) :: dlogr !< logarithmic width of radius bin |
---|
[1871] | 985 | REAL(wp) :: e_a !< vapor pressure |
---|
| 986 | REAL(wp) :: e_s !< saturation vapor pressure |
---|
[2312] | 987 | REAL(wp) :: rmin = 0.005e-6_wp !< minimum aerosol radius |
---|
| 988 | REAL(wp) :: rmax = 10.0e-6_wp !< maximum aerosol radius |
---|
| 989 | REAL(wp) :: r_mid !< mean radius of bin |
---|
| 990 | REAL(wp) :: r_l !< left radius of bin |
---|
| 991 | REAL(wp) :: r_r !< right radius of bin |
---|
[2122] | 992 | REAL(wp) :: sigma !< surface tension |
---|
[1871] | 993 | REAL(wp) :: t_int !< temperature |
---|
| 994 | |
---|
[1890] | 995 | INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: local_start !< |
---|
[1871] | 996 | |
---|
| 997 | INTEGER(iwp) :: n !< |
---|
| 998 | INTEGER(iwp) :: ip !< |
---|
| 999 | INTEGER(iwp) :: jp !< |
---|
| 1000 | INTEGER(iwp) :: kp !< |
---|
| 1001 | |
---|
| 1002 | ! |
---|
[2312] | 1003 | !-- The following typical aerosol spectra are taken from Jaenicke (1993): |
---|
| 1004 | !-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions. |
---|
| 1005 | IF ( TRIM(aero_type) .EQ. 'polar' ) THEN |
---|
| 1006 | na = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6 |
---|
| 1007 | rm = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6 |
---|
| 1008 | log_sigma = (/ 0.245, 0.300, 0.291 /) |
---|
| 1009 | ELSEIF ( TRIM(aero_type) .EQ. 'background' ) THEN |
---|
| 1010 | na = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6 |
---|
| 1011 | rm = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6 |
---|
| 1012 | log_sigma = (/ 0.645, 0.253, 0.425 /) |
---|
| 1013 | ELSEIF ( TRIM(aero_type) .EQ. 'maritime' ) THEN |
---|
| 1014 | na = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6 |
---|
| 1015 | rm = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6 |
---|
| 1016 | log_sigma = (/ 0.657, 0.210, 0.396 /) |
---|
| 1017 | ELSEIF ( TRIM(aero_type) .EQ. 'continental' ) THEN |
---|
| 1018 | na = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6 |
---|
| 1019 | rm = (/ 0.01, 0.058, 0.9 /) * 1.0E-6 |
---|
| 1020 | log_sigma = (/ 0.161, 0.217, 0.380 /) |
---|
| 1021 | ELSEIF ( TRIM(aero_type) .EQ. 'desert' ) THEN |
---|
| 1022 | na = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6 |
---|
| 1023 | rm = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6 |
---|
| 1024 | log_sigma = (/ 0.247, 0.770, 0.438 /) |
---|
| 1025 | ELSEIF ( TRIM(aero_type) .EQ. 'rural' ) THEN |
---|
| 1026 | na = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6 |
---|
| 1027 | rm = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6 |
---|
| 1028 | log_sigma = (/ 0.225, 0.557, 0.266 /) |
---|
| 1029 | ELSEIF ( TRIM(aero_type) .EQ. 'urban' ) THEN |
---|
| 1030 | na = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6 |
---|
| 1031 | rm = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6 |
---|
| 1032 | log_sigma = (/ 0.245, 0.666, 0.337 /) |
---|
| 1033 | ELSEIF ( TRIM(aero_type) .EQ. 'user' ) THEN |
---|
| 1034 | CONTINUE |
---|
| 1035 | ELSE |
---|
| 1036 | WRITE( message_string, * ) 'unknown aerosol type ', & |
---|
| 1037 | 'aero_type = "', TRIM( aero_type ), '"' |
---|
| 1038 | CALL message( 'lpm_init', 'PA0459', 1, 2, 0, 6, 0 ) |
---|
[1871] | 1039 | ENDIF |
---|
| 1040 | |
---|
| 1041 | DO ip = nxl, nxr |
---|
| 1042 | DO jp = nys, nyn |
---|
| 1043 | DO kp = nzb+1, nzt |
---|
| 1044 | |
---|
| 1045 | number_of_particles = prt_count(kp,jp,ip) |
---|
| 1046 | IF ( number_of_particles <= 0 ) CYCLE |
---|
| 1047 | particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) |
---|
[2312] | 1048 | |
---|
| 1049 | dlogr = ( LOG10(rmax) - LOG10(rmin) ) / ( number_of_particles - local_start(kp,jp,ip) + 1 ) |
---|
[1871] | 1050 | ! |
---|
| 1051 | !-- Initialize the aerosols with a predefined spectral distribution |
---|
[2312] | 1052 | !-- of the dry radius (logarithmically increasing bins) and a varying |
---|
[1871] | 1053 | !-- weighting factor |
---|
[2312] | 1054 | DO n = local_start(kp,jp,ip), number_of_particles !only new particles |
---|
[1871] | 1055 | |
---|
[2312] | 1056 | r_l = 10.0**( LOG10( rmin ) + (n-1) * dlogr ) |
---|
| 1057 | r_r = 10.0**( LOG10( rmin ) + n * dlogr ) |
---|
| 1058 | r_mid = SQRT( r_l * r_r ) |
---|
[1871] | 1059 | |
---|
[2312] | 1060 | particles(n)%aux1 = r_mid |
---|
| 1061 | particles(n)%weight_factor = & |
---|
| 1062 | ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) * & |
---|
| 1063 | EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) + & |
---|
| 1064 | na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) * & |
---|
| 1065 | EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) + & |
---|
| 1066 | na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) * & |
---|
| 1067 | EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) ) & |
---|
| 1068 | ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dz ) |
---|
[1871] | 1069 | |
---|
[2312] | 1070 | ! |
---|
| 1071 | !-- Multiply weight_factor with the namelist parameter aero_weight |
---|
| 1072 | !-- to increase or decrease the number of simulated aerosols |
---|
| 1073 | particles(n)%weight_factor = particles(n)%weight_factor * aero_weight |
---|
[1871] | 1074 | |
---|
[2312] | 1075 | IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) & |
---|
| 1076 | .GT. random_function( iran_part ) ) THEN |
---|
| 1077 | particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0 |
---|
| 1078 | ELSE |
---|
| 1079 | particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) |
---|
[1871] | 1080 | ENDIF |
---|
| 1081 | ! |
---|
[2312] | 1082 | !-- Unnecessary particles will be deleted |
---|
| 1083 | IF ( particles(n)%weight_factor .LE. 0.0 ) particles(n)%particle_mask = .FALSE. |
---|
[1871] | 1084 | |
---|
[2312] | 1085 | ENDDO |
---|
[1871] | 1086 | ! |
---|
| 1087 | !-- Set particle radius to equilibrium radius based on the environmental |
---|
| 1088 | !-- supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids |
---|
[2312] | 1089 | !-- the sometimes lengthy growth toward their equilibrium radius within |
---|
[1871] | 1090 | !-- the simulation. |
---|
| 1091 | t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp |
---|
| 1092 | |
---|
| 1093 | e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) ) |
---|
| 1094 | e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp ) |
---|
| 1095 | |
---|
[2122] | 1096 | sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) |
---|
| 1097 | afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) |
---|
| 1098 | |
---|
| 1099 | bfactor = vanthoff * molecular_weight_of_water * & |
---|
| 1100 | rho_s / ( molecular_weight_of_solute * rho_l ) |
---|
[1871] | 1101 | ! |
---|
[2312] | 1102 | !-- The formula is only valid for subsaturated environments. For |
---|
[2122] | 1103 | !-- supersaturations higher than -5 %, the supersaturation is set to -5%. |
---|
| 1104 | IF ( e_a / e_s >= 0.95_wp ) e_a = 0.95_wp * e_s |
---|
[1871] | 1105 | |
---|
[2312] | 1106 | DO n = local_start(kp,jp,ip), number_of_particles !only new particles |
---|
[2122] | 1107 | ! |
---|
[2312] | 1108 | !-- For details on this equation, see Eq. (14) of Khvorostyanov and |
---|
| 1109 | !-- Curry (2007, JGR) |
---|
[2122] | 1110 | particles(n)%radius = bfactor**0.3333333_wp * & |
---|
[2312] | 1111 | particles(n)%aux1 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / & |
---|
[2122] | 1112 | ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp * & |
---|
[2312] | 1113 | particles(n)%aux1 ) ) / & |
---|
[2122] | 1114 | ( 1.0_wp - e_a / e_s )**0.6666666_wp & |
---|
| 1115 | ) |
---|
[1871] | 1116 | |
---|
[2122] | 1117 | ENDDO |
---|
[1871] | 1118 | |
---|
| 1119 | ENDDO |
---|
| 1120 | ENDDO |
---|
| 1121 | ENDDO |
---|
| 1122 | |
---|
| 1123 | END SUBROUTINE lpm_init_aerosols |
---|
| 1124 | |
---|
[1359] | 1125 | END MODULE lpm_init_mod |
---|