Changeset 1822 for palm/trunk/SOURCE/lpm_advec.f90
- Timestamp:
- Apr 7, 2016 7:49:42 AM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_advec.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r1818 r1822 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Random velocity fluctuations for particles added. Terminal fall velocity 22 ! for droplets is calculated from a parameterization (which is better than 23 ! the previous, physically correct calculation, which demands a very short 24 ! time step that is not used in the model). 25 ! 26 ! Unused variables deleted. 22 27 ! 23 28 ! Former revisions: … … 78 83 79 84 USE arrays_3d, & 80 ONLY: de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu, & 81 zw 85 ONLY: de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw 82 86 83 87 USE cpulog … … 87 91 USE control_parameters, & 88 92 ONLY: atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d, & 89 dt_3d_reached_l, dz, g, kappa, molecular_viscosity, topography, & 90 u_gtrans, v_gtrans, simulated_time 93 dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans 91 94 92 95 USE grid_variables, & … … 99 102 100 103 USE particle_attributes, & 101 ONLY: block_offset, c_0, d ensity_ratio, dt_min_part, grid_particles,&104 ONLY: block_offset, c_0, dt_min_part, grid_particles, & 102 105 iran_part, log_z_z0, number_of_particles, number_of_sublayers, & 103 106 particles, particle_groups, offset_ocean_nzt, sgs_wfu_part, & … … 140 143 REAL(wp) :: de_dz_int_l !< 141 144 REAL(wp) :: de_dz_int_u !< 145 REAL(wp) :: diameter !< diamter of droplet 142 146 REAL(wp) :: diss_int_l !< 143 147 REAL(wp) :: diss_int_u !< … … 150 154 REAL(wp) :: exp_term !< 151 155 REAL(wp) :: gg !< 152 REAL(wp) :: height_int !<153 156 REAL(wp) :: height_p !< 154 REAL(wp) :: lagr_timescale !< 157 REAL(wp) :: lagr_timescale !< Lagrangian timescale 155 158 REAL(wp) :: location(1:30,1:3) !< 156 159 REAL(wp) :: random_gauss !< 160 REAL(wp) :: RL !< Lagrangian autocorrelation coefficient 161 REAL(wp) :: rg1 !< Gaussian distributed random number 162 REAL(wp) :: rg2 !< Gaussian distributed random number 163 REAL(wp) :: rg3 !< Gaussian distributed random number 164 REAL(wp) :: sigma !< velocity standard deviation 157 165 REAL(wp) :: u_int_l !< 158 166 REAL(wp) :: u_int_u !< … … 163 171 REAL(wp) :: w_int_l !< 164 172 REAL(wp) :: w_int_u !< 173 REAL(wp) :: w_s !< terminal velocity of droplets 165 174 REAL(wp) :: x !< 166 175 REAL(wp) :: y !< 167 REAL(wp) :: z_p !< 176 REAL(wp) :: z_p !< 177 178 REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity 179 REAL(wp), PARAMETER :: b_rog = 10.43_wp !< parameter for fall velocity 180 REAL(wp), PARAMETER :: c_rog = 0.6_wp !< parameter for fall velocity 181 REAL(wp), PARAMETER :: k_cap_rog = 4.0_wp !< parameter for fall velocity 182 REAL(wp), PARAMETER :: k_low_rog = 12.0_wp !< parameter for fall velocity 183 REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter 168 184 169 185 REAL(wp), DIMENSION(1:30) :: d_gp_pl !< … … 388 404 !-- Interpolate and calculate quantities needed for calculating the SGS 389 405 !-- velocities 390 IF ( use_sgs_for_particles ) THEN406 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 391 407 392 408 IF ( topography == 'flat' ) THEN … … 1336 1352 !-- Update of the particle velocity 1337 1353 IF ( cloud_droplets ) THEN 1338 exp_arg = 4.5_wp * dens_ratio(n) * molecular_viscosity / & 1339 ( particles(n)%radius )**2 * & 1340 ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius & 1341 * SQRT( ( u_int(n) - particles(n)%speed_x )**2 + & 1342 ( v_int(n) - particles(n)%speed_y )**2 + & 1343 ( w_int(n) - particles(n)%speed_z )**2 ) & 1344 / molecular_viscosity )**0.687_wp & 1345 ) 1346 1347 exp_term = EXP( -exp_arg * dt_particle(n) ) 1348 ELSEIF ( use_sgs_for_particles ) THEN 1354 ! 1355 !-- Terminal velocity is computed for vertical direction (Rogers et 1356 !-- al., 1993, J. Appl. Meteorol.) 1357 diameter = particles(n)%radius * 2000.0_wp !diameter in mm 1358 IF ( diameter <= d0_rog ) THEN 1359 w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) ) 1360 ELSE 1361 w_s = a_rog - b_rog * EXP( -c_rog * diameter ) 1362 ENDIF 1363 1364 ! 1365 !-- If selected, add random velocities following Soelch and Kaercher 1366 !-- (2010, Q. J. R. Meteorol. Soc.) 1367 IF ( use_sgs_for_particles ) THEN 1368 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 1369 RL = EXP( -1.0_wp * dt_3d / lagr_timescale ) 1370 sigma = SQRT( e(kp,jp,ip) ) 1371 1372 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1373 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1374 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1375 1376 particles(n)%rvar1 = RL * particles(n)%rvar1 + & 1377 SQRT( 1.0_wp - RL**2 ) * sigma * rg1 1378 particles(n)%rvar2 = RL * particles(n)%rvar2 + & 1379 SQRT( 1.0_wp - RL**2 ) * sigma * rg2 1380 particles(n)%rvar3 = RL * particles(n)%rvar3 + & 1381 SQRT( 1.0_wp - RL**2 ) * sigma * rg3 1382 1383 particles(n)%speed_x = u_int(n) + particles(n)%rvar1 1384 particles(n)%speed_y = v_int(n) + particles(n)%rvar2 1385 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s 1386 ELSE 1387 particles(n)%speed_x = u_int(n) 1388 particles(n)%speed_y = v_int(n) 1389 particles(n)%speed_z = w_int(n) - w_s 1390 ENDIF 1391 1392 ELSE 1393 1394 IF ( use_sgs_for_particles ) THEN 1395 exp_arg = particle_groups(particles(n)%group)%exp_arg 1396 exp_term = EXP( -exp_arg * dt_particle(n) ) 1397 ELSE 1398 exp_arg = particle_groups(particles(n)%group)%exp_arg 1399 exp_term = particle_groups(particles(n)%group)%exp_term 1400 ENDIF 1401 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 1402 u_int(n) * ( 1.0_wp - exp_term ) 1403 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 1404 v_int(n) * ( 1.0_wp - exp_term ) 1405 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 1406 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * & 1407 g / exp_arg ) * ( 1.0_wp - exp_term ) 1408 ENDIF 1409 1410 ENDIF 1411 1412 ENDDO 1413 1414 ELSE 1415 1416 DO n = 1, number_of_particles 1417 1418 !-- Transport of particles with inertia 1419 particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n) 1420 particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n) 1421 particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n) 1422 ! 1423 !-- Update of the particle velocity 1424 IF ( cloud_droplets ) THEN 1425 ! 1426 !-- Terminal velocity is computed for vertical direction (Rogers et al., 1427 !-- 1993, J. Appl. Meteorol.) 1428 diameter = particles(n)%radius * 2000.0_wp !diameter in mm 1429 IF ( diameter <= d0_rog ) THEN 1430 w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) ) 1431 ELSE 1432 w_s = a_rog - b_rog * EXP( -c_rog * diameter ) 1433 ENDIF 1434 1435 ! 1436 !-- If selected, add random velocities following Soelch and Kaercher 1437 !-- (2010, Q. J. R. Meteorol. Soc.) 1438 IF ( use_sgs_for_particles ) THEN 1439 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 1440 RL = EXP( -1.0_wp * dt_3d / lagr_timescale ) 1441 sigma = SQRT( e(kp,jp,ip) ) 1442 1443 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1444 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1445 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp 1446 1447 particles(n)%rvar1 = RL * particles(n)%rvar1 + & 1448 SQRT( 1.0_wp - RL**2 ) * sigma * rg1 1449 particles(n)%rvar2 = RL * particles(n)%rvar2 + & 1450 SQRT( 1.0_wp - RL**2 ) * sigma * rg2 1451 particles(n)%rvar3 = RL * particles(n)%rvar3 + & 1452 SQRT( 1.0_wp - RL**2 ) * sigma * rg3 1453 1454 particles(n)%speed_x = u_int(n) + particles(n)%rvar1 1455 particles(n)%speed_y = v_int(n) + particles(n)%rvar2 1456 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s 1457 ELSE 1458 particles(n)%speed_x = u_int(n) 1459 particles(n)%speed_y = v_int(n) 1460 particles(n)%speed_z = w_int(n) - w_s 1461 ENDIF 1462 1463 ELSE 1464 1465 IF ( use_sgs_for_particles ) THEN 1349 1466 exp_arg = particle_groups(particles(n)%group)%exp_arg 1350 1467 exp_term = EXP( -exp_arg * dt_particle(n) ) … … 1353 1470 exp_term = particle_groups(particles(n)%group)%exp_term 1354 1471 ENDIF 1355 particles(n)%speed_x = particles(n)%speed_x * exp_term + &1472 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 1356 1473 u_int(n) * ( 1.0_wp - exp_term ) 1357 particles(n)%speed_y = particles(n)%speed_y * exp_term + &1474 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 1358 1475 v_int(n) * ( 1.0_wp - exp_term ) 1359 particles(n)%speed_z = particles(n)%speed_z * exp_term + &1360 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &1361 g /exp_arg ) * ( 1.0_wp - exp_term )1476 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 1477 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / & 1478 exp_arg ) * ( 1.0_wp - exp_term ) 1362 1479 ENDIF 1363 1480 1364 ENDDO1365 1366 ELSE1367 1368 DO n = 1, number_of_particles1369 1370 !-- Transport of particles with inertia1371 particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)1372 particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)1373 particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)1374 !1375 !-- Update of the particle velocity1376 IF ( cloud_droplets ) THEN1377 1378 exp_arg = 4.5_wp * dens_ratio(n) * molecular_viscosity / &1379 ( particles(n)%radius )**2 * &1380 ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius * &1381 SQRT( ( u_int(n) - particles(n)%speed_x )**2 + &1382 ( v_int(n) - particles(n)%speed_y )**2 + &1383 ( w_int(n) - particles(n)%speed_z )**2 ) / &1384 molecular_viscosity )**0.687_wp &1385 )1386 1387 exp_term = EXP( -exp_arg * dt_particle(n) )1388 ELSEIF ( use_sgs_for_particles ) THEN1389 exp_arg = particle_groups(particles(n)%group)%exp_arg1390 exp_term = EXP( -exp_arg * dt_particle(n) )1391 ELSE1392 exp_arg = particle_groups(particles(n)%group)%exp_arg1393 exp_term = particle_groups(particles(n)%group)%exp_term1394 ENDIF1395 particles(n)%speed_x = particles(n)%speed_x * exp_term + &1396 u_int(n) * ( 1.0_wp - exp_term )1397 particles(n)%speed_y = particles(n)%speed_y * exp_term + &1398 v_int(n) * ( 1.0_wp - exp_term )1399 particles(n)%speed_z = particles(n)%speed_z * exp_term + &1400 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &1401 exp_arg ) * ( 1.0_wp - exp_term )1402 1481 ENDDO 1403 1482
Note: See TracChangeset
for help on using the changeset viewer.
