Changeset 4430 for palm/trunk/SOURCE
- Timestamp:
- Feb 27, 2020 6:02:20 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4360 r4430 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Bugfix in logarithmic interpolation of near-ground particle speed (density 28 ! was not considered). 29 ! - Revise CFL-check when SGS particle speeds are considered. 30 ! - In nested case with SGS particle speeds in the child domain, do not give 31 ! warning that particles are on domain boundaries. At the end of the particle 32 ! time integration these will be transferred to the parent domain anyhow. 33 ! 34 ! 4360 2020-01-07 11:25:50Z suehring 27 35 ! Introduction of wall_flags_total_0, which currently sets bits based on static 28 36 ! topography information used in wall_flags_static_0 … … 126 134 127 135 USE arrays_3d, & 128 ONLY: de_dx, de_dy, de_dz, dzw, zu, zw, ql_c, ql_v, ql_vp, hyp, & 129 pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p, & 130 d_exner 136 ONLY: de_dx, de_dy, de_dz, & 137 d_exner, & 138 drho_air_zw, & 139 dzw, zu, zw, ql_c, ql_v, ql_vp, hyp, & 140 pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p 131 141 132 142 USE averaging, & … … 139 149 USE control_parameters, & 140 150 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, & 151 child_domain, & 141 152 cloud_droplets, constant_flux_layer, current_timestep_number, & 142 153 dt_3d, dt_3d_reached, first_call_lpm, humidity, & … … 830 841 ! 831 842 !-- Liquid water content, change in liquid water content 832 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &843 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 833 844 ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 834 !835 845 !-- Real volume of particles (with weighting), volume of particles 836 ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &837 846 ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 847 ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 838 848 ENDIF 839 849 840 850 841 ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &842 v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &851 ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 852 v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 843 853 w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 844 854 ! … … 2338 2348 !------------------------------------------------------------------------------! 2339 2349 SUBROUTINE particles_from_parent_to_child 2340 IMPLICIT NONE2341 2350 2342 2351 CALL pmcp_c_get_particle_from_parent ! Child actions … … 2344 2353 2345 2354 RETURN 2355 2346 2356 END SUBROUTINE particles_from_parent_to_child 2347 2357 … … 2353 2363 !------------------------------------------------------------------------------! 2354 2364 SUBROUTINE particles_from_child_to_parent 2355 IMPLICIT NONE2356 2365 2357 2366 CALL pmcp_c_send_particle_to_parent ! Child actions … … 2359 2368 2360 2369 RETURN 2370 2361 2371 END SUBROUTINE particles_from_child_to_parent 2362 2372 … … 3315 3325 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 3316 3326 3317 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers3327 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers 3318 3328 3319 3329 CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' ) … … 3461 3471 ! 3462 3472 !-- Determine the sublayer. Further used as index. 3463 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &3464 * REAL( number_of_sublayers, KIND=wp ) 3473 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3474 * REAL( number_of_sublayers, KIND=wp ) & 3465 3475 * d_z_p_z0 3466 3476 ! 3467 3477 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3468 3478 !-- interpolate linearly between precalculated logarithm. 3469 log_z_z0_int = log_z_z0(INT(height_p)) 3470 + ( height_p - INT(height_p) ) 3471 * ( log_z_z0(INT(height_p)+1) 3472 - log_z_z0(INT(height_p)) 3479 log_z_z0_int = log_z_z0(INT(height_p)) & 3480 + ( height_p - INT(height_p) ) & 3481 * ( log_z_z0(INT(height_p)+1) & 3482 - log_z_z0(INT(height_p)) & 3473 3483 ) 3474 3484 ! 3475 3485 !-- Get friction velocity and momentum flux from new surface data 3476 3486 !-- types. 3477 IF ( surf_def_h(0)%start_index(jp,ip) <= &3487 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3478 3488 surf_def_h(0)%end_index(jp,ip) ) THEN 3479 3489 surf_start = surf_def_h(0)%start_index(jp,ip) … … 3482 3492 !-- large particle speed. 3483 3493 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3484 usws_int = surf_def_h(0)%usws(surf_start) 3485 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3494 usws_int = surf_def_h(0)%usws(surf_start) & 3495 * drho_air_zw(k_wall) 3496 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3486 3497 surf_lsm_h%end_index(jp,ip) ) THEN 3487 3498 surf_start = surf_lsm_h%start_index(jp,ip) 3488 3499 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3489 usws_int = surf_lsm_h%usws(surf_start) 3490 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3500 usws_int = surf_lsm_h%usws(surf_start) & 3501 * drho_air_zw(k_wall) 3502 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3491 3503 surf_usm_h%end_index(jp,ip) ) THEN 3492 3504 surf_start = surf_usm_h%start_index(jp,ip) 3493 3505 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3494 usws_int = surf_usm_h%usws(surf_start) 3506 usws_int = surf_usm_h%usws(surf_start) & 3507 * drho_air_zw(k_wall) 3495 3508 ENDIF 3496 3509 ! … … 3501 3514 !-- as sensitivity studies revealed no significant effect of 3502 3515 !-- using the neutral solution also for un/stable situations. 3503 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) 3516 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 3504 3517 * log_z_z0_int - u_gtrans 3505 3518 ENDIF … … 3553 3566 !-- topography particle on u-grid can be above surface-layer height, 3554 3567 !-- whereas it can be below on v-grid. 3555 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &3556 * REAL( number_of_sublayers, KIND=wp ) 3568 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3569 * REAL( number_of_sublayers, KIND=wp ) & 3557 3570 * d_z_p_z0 3558 3571 ! 3559 3572 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3560 3573 !-- interpolate linearly between precalculated logarithm. 3561 log_z_z0_int = log_z_z0(INT(height_p)) 3562 + ( height_p - INT(height_p) ) 3563 * ( log_z_z0(INT(height_p)+1) 3564 - log_z_z0(INT(height_p)) 3574 log_z_z0_int = log_z_z0(INT(height_p)) & 3575 + ( height_p - INT(height_p) ) & 3576 * ( log_z_z0(INT(height_p)+1) & 3577 - log_z_z0(INT(height_p)) & 3565 3578 ) 3566 3579 ! 3567 3580 !-- Get friction velocity and momentum flux from new surface data 3568 3581 !-- types. 3569 IF ( surf_def_h(0)%start_index(jp,ip) <= &3582 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3570 3583 surf_def_h(0)%end_index(jp,ip) ) THEN 3571 3584 surf_start = surf_def_h(0)%start_index(jp,ip) … … 3574 3587 !-- large particle speed. 3575 3588 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3576 vsws_int = surf_def_h(0)%vsws(surf_start) 3577 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3589 vsws_int = surf_def_h(0)%vsws(surf_start) & 3590 * drho_air_zw(k_wall) 3591 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3578 3592 surf_lsm_h%end_index(jp,ip) ) THEN 3579 3593 surf_start = surf_lsm_h%start_index(jp,ip) 3580 3594 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3581 vsws_int = surf_lsm_h%vsws(surf_start) 3582 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3595 vsws_int = surf_lsm_h%vsws(surf_start) & 3596 * drho_air_zw(k_wall) 3597 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3583 3598 surf_usm_h%end_index(jp,ip) ) THEN 3584 3599 surf_start = surf_usm_h%start_index(jp,ip) 3585 3600 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3586 vsws_int = surf_usm_h%vsws(surf_start) 3601 vsws_int = surf_usm_h%vsws(surf_start) & 3602 * drho_air_zw(k_wall) 3587 3603 ENDIF 3588 3604 ! … … 3593 3609 !-- as sensitivity studies revealed no significant effect of 3594 3610 !-- using the neutral solution also for un/stable situations. 3595 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) 3611 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 3596 3612 * log_z_z0_int - v_gtrans 3597 3613 … … 3901 3917 ENDIF 3902 3918 ENDIF 3919 3903 3920 rvar1_temp(n) = particles(n)%rvar1 3904 3921 rvar2_temp(n) = particles(n)%rvar2 … … 3918 3935 rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 3919 3936 + 1E-20_wp ) * ( rg(n,3) - 1.0_wp ) 3920 3921 3937 ELSE 3922 3938 ! … … 3966 3982 ! 3967 3983 !-- Check if the added SGS velocities result in a violation of the CFL- 3968 !-- criterion. If yes choose a smaller timestep based on the new velocities 3969 !-- and calculate SGS velocities again 3984 !-- criterion. If yes, limt the SGS particle speed to match the 3985 !-- CFL criterion. Note, a re-calculation of the SGS particle speed with 3986 !-- smaller timestep does not necessarily fulfill the CFL criterion as the 3987 !-- new SGS speed can be even larger (due to the random term with scales with 3988 !-- the square-root of dt_particle, for small dt the random contribution increases). 3989 !-- Thus, we would need to re-calculate the SGS speeds as long as they would 3990 !-- fulfill the requirements, which could become computationally expensive, 3991 !-- Hence, we just limit them. 3970 3992 dz_temp = zw(kp)-zw(kp-1) 3971 3993 3972 3994 DO nb = 0, 7 3973 3995 DO n = start_index(nb), end_index(nb) 3974 IF ( .NOT. particles(n)%age == 0.0_wp .AND. & 3975 (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n)) .OR. & 3976 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n)) .OR. & 3977 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN 3978 3979 dt_particle(n) = 0.9_wp * MIN( & 3980 ( dx / ABS( u_int(n) + rvar1_temp(n) ) ), & 3981 ( dy / ABS( v_int(n) + rvar2_temp(n) ) ), & 3982 ( dz_temp / ABS( w_int(n) + rvar3_temp(n) ) ) ) 3983 3984 ! 3985 !-- Reset temporary SGS velocites to "current" ones 3986 rvar1_temp(n) = particles(n)%rvar1 3987 rvar2_temp(n) = particles(n)%rvar2 3988 rvar3_temp(n) = particles(n)%rvar3 3989 3990 de_dt_min = - e_int(n) / dt_particle(n) 3991 3992 de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m 3993 3994 IF ( de_dt < de_dt_min ) THEN 3995 de_dt = de_dt_min 3996 ENDIF 3997 3998 CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n), & 3999 de_dx_int(n), de_dt, diss_int(n), & 4000 dt_particle(n), rg(n,1), term_1_2(n) ) 4001 4002 CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n), & 4003 de_dy_int(n), de_dt, diss_int(n), & 4004 dt_particle(n), rg(n,2), term_1_2(n) ) 4005 4006 CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n), & 4007 de_dz_int(n), de_dt, diss_int(n), & 4008 dt_particle(n), rg(n,3), term_1_2(n) ) 3996 IF ( ABS( u_int(n) + rvar1_temp(n) ) > ( dx / dt_particle(n) ) .OR. & 3997 ABS( v_int(n) + rvar2_temp(n) ) > ( dy / dt_particle(n) ) .OR. & 3998 ABS( w_int(n) + rvar3_temp(n) ) > ( dz_temp / dt_particle(n) ) ) THEN 3999 ! 4000 !-- If total speed exceeds the allowed speed according to CFL 4001 !-- criterion, limit the SGS speed to 4002 !-- dx_i / dt_particle - u_resolved_i, considering a safty factor. 4003 rvar1_temp(n) = MERGE( rvar1_temp(n), & 4004 0.9_wp * & 4005 SIGN( dx / dt_particle(n) & 4006 - ABS( u_int(n) ), rvar1_temp(n) ), & 4007 ABS( u_int(n) + rvar1_temp(n) ) < & 4008 ( dx / dt_particle(n) ) ) 4009 rvar2_temp(n) = MERGE( rvar2_temp(n), & 4010 0.9_wp * & 4011 SIGN( dy / dt_particle(n) & 4012 - ABS( v_int(n) ), rvar2_temp(n) ), & 4013 ABS( v_int(n) + rvar2_temp(n) ) < & 4014 ( dy / dt_particle(n) ) ) 4015 rvar3_temp(n) = MERGE( rvar3_temp(n), & 4016 0.9_wp * & 4017 SIGN( zw(kp)-zw(kp-1) / dt_particle(n) & 4018 - ABS( w_int(n) ), rvar3_temp(n) ), & 4019 ABS( w_int(n) + rvar3_temp(n) ) < & 4020 ( zw(kp)-zw(kp-1) / dt_particle(n) ) ) 4009 4021 ENDIF 4010 4011 4022 ! 4012 4023 !-- Update particle velocites … … 4276 4287 ENDIF 4277 4288 ! 4278 !-- Loop from particle start to particle end 4289 !-- Loop from particle start to particle end and increment the particle 4290 !-- age and the total time that the particle has advanced within the 4291 !-- particle timestep procedure. 4279 4292 DO n = particle_start, particle_end 4280 !4281 !-- Increment the particle age and the total time that the particle4282 !-- has advanced within the particle timestep procedure4283 4293 particles(n)%age = particles(n)%age + dt_particle(n) 4284 4294 particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n) 4285 4286 ! 4287 !-- Check whether there is still a particle that has not yet completed 4288 !-- the total LES timestep 4289 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN 4295 ENDDO 4296 ! 4297 !-- Particles that leave the child domain during the SGS-timestep loop 4298 !-- must not continue timestepping until they are transferred to the 4299 !-- parent. Hence, set their dt_sum to dt. 4300 IF ( child_domain .AND. use_sgs_for_particles ) THEN 4301 DO n = particle_start, particle_end 4302 IF ( particles(n)%x < 0.0_wp .OR. & 4303 particles(n)%y < 0.0_wp .OR. & 4304 particles(n)%x > ( nx+1 ) * dx .OR. & 4305 particles(n)%y < ( ny+1 ) * dy ) THEN 4306 particles(n)%dt_sum = dt_3d 4307 ENDIF 4308 ENDDO 4309 ENDIF 4310 ! 4311 !-- Check whether there is still a particle that has not yet completed 4312 !-- the total LES timestep 4313 DO n = particle_start, particle_end 4314 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) & 4290 4315 dt_3d_reached_l = .FALSE. 4291 ENDIF4292 4293 4316 ENDDO 4294 4317 ENDDO … … 7060 7083 IF ( i < 0 ) THEN 7061 7084 ! 7062 !-- Apply boundary condition along x7085 !-- Apply boundary condition along x 7063 7086 IF ( ibc_par_lr == 0 ) THEN 7064 7087 ! … … 7440 7463 7441 7464 IF ( trsp_count_recv > 0 ) CALL lpm_add_particles_to_gridcell(rvsp(1:trsp_count_recv)) 7465 7442 7466 7443 7467 DEALLOCATE(rvsp) … … 7742 7766 ENDIF 7743 7767 ELSE 7744 WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn 7768 IF ( .NOT. child_domain ) THEN 7769 WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn 7770 ENDIF 7745 7771 ENDIF 7746 7772 ENDIF … … 7863 7889 !-- Note, check for CFL does not work at first particle timestep 7864 7890 !-- when both, age and age_m are zero. 7865 IF ( particles(n)%age - particles(n)%age_m > 0.0_wp ) THEN 7866 IF( ABS(particles(n)%speed_x) >&7867 ( dx/(particles(n)%age-particles(n)%age_m)) .OR.&7868 ABS(particles(n)%speed_y) >&7869 ( dy/(particles(n)%age-particles(n)%age_m)) .OR.&7870 ABS(particles(n)%speed_z) >&7871 ( (zw(k)-zw(k-1))/(particles(n)%age-particles(n)%age_m)))&7872 THEN7891 IF ( particles(n)%age - particles(n)%age_m > 0.0_wp ) THEN 7892 IF( ABS( particles(n)%speed_x ) > & 7893 ( dx / ( particles(n)%age - particles(n)%age_m) ) .OR. & 7894 ABS( particles(n)%speed_y ) > & 7895 ( dy / ( particles(n)%age - particles(n)%age_m) ) .OR. & 7896 ABS( particles(n)%speed_z ) > & 7897 ( ( zw(k)-zw(k-1) ) & 7898 / ( particles(n)%age - particles(n)%age_m) ) ) THEN 7873 7899 WRITE( message_string, * ) & 7874 7900 'Particle violated CFL-criterion: &particle with id ', & 7875 7901 particles(n)%id, ' will be deleted!' 7876 7902 CALL message( 'lpm_check_cfl', 'PA0475', 0, 1, -1, 6, 0 ) 7903 7877 7904 particles(n)%particle_mask= .FALSE. 7878 7905 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
