Changeset 2886 for palm/trunk/SOURCE/lpm_advec.f90
- Timestamp:
- Mar 14, 2018 11:51:53 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_advec.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r2718 r2886 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Bugfix in passive particle SGS Model: 28 ! Sometimes the added SGS velocities would lead to a violation of the CFL 29 ! criterion for single particles. For this a check was added after the 30 ! calculation of SGS velocities. 31 ! 32 ! 2718 2018-01-02 08:49:38Z maronga 27 33 ! Corrected "Former revisions" section 28 34 ! … … 218 224 REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level 219 225 REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level 220 REAL(wp) :: dt_gap !< remaining time until particle time integration reaches LES time221 226 REAL(wp) :: dt_particle_m !< previous particle time step 227 REAL(wp) :: dz_temp !< 222 228 REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level 223 229 REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level … … 227 233 REAL(wp) :: gg !< dummy argument for horizontal particle interpolation 228 234 REAL(wp) :: height_p !< dummy argument for logarithmic interpolation 229 REAL(wp) :: lagr_timescale !< Lagrangian timescale230 235 REAL(wp) :: location(1:30,1:3) !< wall locations 231 236 REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation … … 265 270 REAL(wp), DIMENSION(1:30) :: ei !< TKE at adjacent wall 266 271 267 REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not 268 REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !< 269 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position 270 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position 271 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position 272 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position 273 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step 274 REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position 275 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed 276 REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed 277 REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed 278 REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed 279 REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position 280 REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position 281 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 272 REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not 273 REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !< 274 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position 275 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position 276 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position 277 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position 278 REAL(wp), DIMENSION(number_of_particles) :: dt_gap !< remaining time until particle time integration reaches LES time 279 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step 280 REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position 281 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed 282 REAL(wp), DIMENSION(number_of_particles) :: lagr_timescale !< Lagrangian timescale 283 REAL(wp), DIMENSION(number_of_particles) :: rvar1_temp !< 284 REAL(wp), DIMENSION(number_of_particles) :: rvar2_temp !< 285 REAL(wp), DIMENSION(number_of_particles) :: rvar3_temp !< 286 REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed 287 REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed 288 REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed 289 REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position 290 REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position 291 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 282 292 283 293 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers … … 763 773 ! 764 774 !-- Calculate the Lagrangian timescale according to Weil et al. (2004). 765 lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &775 lagr_timescale(n) = ( 4.0_wp * e_int(n) + 1E-20_wp ) / & 766 776 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp ) 767 777 … … 769 779 !-- Calculate the next particle timestep. dt_gap is the time needed to 770 780 !-- complete the current LES timestep. 771 dt_gap = dt_3d - particles(n)%dt_sum 772 dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap ) 773 781 dt_gap(n) = dt_3d - particles(n)%dt_sum 782 dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale(n), dt_gap(n) ) 783 particles(n)%aux1 = lagr_timescale(n) 784 particles(n)%aux2 = dt_gap(n) 774 785 ! 775 786 !-- The particle timestep should not be too small in order to prevent 776 787 !-- the number of particle timesteps of getting too large 777 IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap ) THEN788 IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap(n) ) THEN 778 789 dt_particle(n) = dt_min_part 779 790 ENDIF 780 791 rvar1_temp(n) = particles(n)%rvar1 792 rvar2_temp(n) = particles(n)%rvar2 793 rvar3_temp(n) = particles(n)%rvar3 781 794 ! 782 795 !-- Calculate the SGS velocity components … … 787 800 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 788 801 !-- from becoming unrealistically large. 789 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 790 + 1E-20_wp ) * & 791 ( rg(n,1) - 1.0_wp ) 792 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 793 + 1E-20_wp ) * & 794 ( rg(n,2) - 1.0_wp ) 795 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 796 + 1E-20_wp ) * & 797 ( rg(n,3) - 1.0_wp ) 802 rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 803 + 1E-20_wp ) * ( rg(n,1) - 1.0_wp ) 804 rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 805 + 1E-20_wp ) * ( rg(n,2) - 1.0_wp ) 806 rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) & 807 + 1E-20_wp ) * ( rg(n,3) - 1.0_wp ) 798 808 799 809 ELSE … … 811 821 ENDIF 812 822 813 !814 823 !-- For old particles the SGS components are correlated with the 815 824 !-- values from the previous timestep. Random numbers have also to … … 827 836 ENDIF 828 837 829 CALL weil_stochastic_eq( particles(n)%rvar1, fs_int(n), e_int(n),&838 CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 830 839 de_dx_int(n), de_dt, diss_int(n), & 831 840 dt_particle(n), rg(n,1), term_1_2(n) ) 832 841 833 CALL weil_stochastic_eq( particles(n)%rvar2, fs_int(n), e_int(n),&842 CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 834 843 de_dy_int(n), de_dt, diss_int(n), & 835 844 dt_particle(n), rg(n,2), term_1_2(n) ) 836 845 837 CALL weil_stochastic_eq( particles(n)%rvar3, fs_int(n), e_int(n),&846 CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 838 847 de_dz_int(n), de_dt, diss_int(n), & 839 848 dt_particle(n), rg(n,3), term_1_2(n) ) … … 841 850 ENDIF 842 851 852 ENDDO 853 ENDDO 854 ! 855 !-- Check if the added SGS velocities result in a violation of the CFL- 856 !-- criterion. If yes choose a smaller timestep based on the new velocities 857 !-- and calculate SGS velocities again 858 dz_temp = zw(kp)-zw(kp-1) 859 860 DO nb = 0, 7 861 DO n = start_index(nb), end_index(nb) 862 IF ( .NOT. particles(n)%age == 0.0_wp .AND. & 863 (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n)) .OR. & 864 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n)) .OR. & 865 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN 866 867 dt_particle(n) = 0.9_wp * MIN( & 868 ( dx / ABS( u_int(n) + rvar1_temp(n) ) ), & 869 ( dy / ABS( v_int(n) + rvar2_temp(n) ) ), & 870 ( dz_temp / ABS( w_int(n) + rvar3_temp(n) ) ) ) 871 872 ! 873 !-- Reset temporary SGS velocites to "current" ones 874 rvar1_temp(n) = particles(n)%rvar1 875 rvar2_temp(n) = particles(n)%rvar2 876 rvar3_temp(n) = particles(n)%rvar3 877 878 de_dt_min = - e_int(n) / dt_particle(n) 879 880 de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m 881 882 IF ( de_dt < de_dt_min ) THEN 883 de_dt = de_dt_min 884 ENDIF 885 886 CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 887 de_dx_int(n), de_dt, diss_int(n), & 888 dt_particle(n), rg(n,1), term_1_2(n) ) 889 890 CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 891 de_dy_int(n), de_dt, diss_int(n), & 892 dt_particle(n), rg(n,2), term_1_2(n) ) 893 894 CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 895 de_dz_int(n), de_dt, diss_int(n), & 896 dt_particle(n), rg(n,3), term_1_2(n) ) 897 ENDIF 898 899 ! 900 !-- Update particle velocites 901 particles(n)%rvar1 = rvar1_temp(n) 902 particles(n)%rvar2 = rvar2_temp(n) 903 particles(n)%rvar3 = rvar3_temp(n) 843 904 u_int(n) = u_int(n) + particles(n)%rvar1 844 905 v_int(n) = v_int(n) + particles(n)%rvar2 … … 850 911 ENDDO 851 912 ENDDO 852 913 853 914 ELSE 854 915 ! … … 905 966 !-- (2010, Q. J. R. Meteorol. Soc.) 906 967 IF ( use_sgs_for_particles ) THEN 907 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )908 RL = EXP( -1.0_wp * dt_3d / lagr_timescale )968 lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 969 RL = EXP( -1.0_wp * dt_3d / lagr_timescale(n) ) 909 970 sigma = SQRT( e(kp,jp,ip) ) 910 971 … … 977 1038 !-- (2010, Q. J. R. Meteorol. Soc.) 978 1039 IF ( use_sgs_for_particles ) THEN 979 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )980 RL = EXP( -1.0_wp * dt_3d / lagr_timescale )1040 lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp ) 1041 RL = EXP( -1.0_wp * dt_3d / lagr_timescale(n) ) 981 1042 sigma = SQRT( e(kp,jp,ip) ) 982 1043
Note: See TracChangeset
for help on using the changeset viewer.
