Changeset 2886 for palm/trunk


Ignore:
Timestamp:
Mar 14, 2018 11:51:53 AM (4 years ago)
Author:
thiele
Message:

Bugfix in passive particle SGS Model

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_advec.f90

    r2718 r2886  
    2525! -----------------
    2626! $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
    2733! Corrected "Former revisions" section
    2834!
     
    218224    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
    219225    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 time
    221226    REAL(wp) ::  dt_particle_m      !< previous particle time step
     227    REAL(wp) ::  dz_temp            !<
    222228    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
    223229    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
     
    227233    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
    228234    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
    229     REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
    230235    REAL(wp) ::  location(1:30,1:3) !< wall locations
    231236    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
     
    265270    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
    266271
    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
    282292
    283293    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
     
    763773!
    764774!--          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 ) / &
    766776                              ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
    767777
     
    769779!--          Calculate the next particle timestep. dt_gap is the time needed to
    770780!--          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)
    774785!
    775786!--          The particle timestep should not be too small in order to prevent
    776787!--          the number of particle timesteps of getting too large
    777              IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
     788             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap(n) )  THEN
    778789                dt_particle(n) = dt_min_part
    779790             ENDIF
    780 
     791             rvar1_temp(n) = particles(n)%rvar1
     792             rvar2_temp(n) = particles(n)%rvar2
     793             rvar3_temp(n) = particles(n)%rvar3
    781794!
    782795!--          Calculate the SGS velocity components
     
    787800!--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
    788801!--             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 )
    798808
    799809             ELSE
     
    811821                ENDIF
    812822
    813 !
    814823!--             For old particles the SGS components are correlated with the
    815824!--             values from the previous timestep. Random numbers have also to
     
    827836                ENDIF
    828837
    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),&
    830839                                        de_dx_int(n), de_dt, diss_int(n),       &
    831840                                        dt_particle(n), rg(n,1), term_1_2(n) )
    832841
    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),&
    834843                                        de_dy_int(n), de_dt, diss_int(n),       &
    835844                                        dt_particle(n), rg(n,2), term_1_2(n) )
    836845
    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),&
    838847                                        de_dz_int(n), de_dt, diss_int(n),       &
    839848                                        dt_particle(n), rg(n,3), term_1_2(n) )
     
    841850             ENDIF
    842851
     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)
    843904             u_int(n) = u_int(n) + particles(n)%rvar1
    844905             v_int(n) = v_int(n) + particles(n)%rvar2
     
    850911          ENDDO
    851912       ENDDO
    852 
     913       
    853914    ELSE
    854915!
     
    905966!--                (2010, Q. J. R. Meteorol. Soc.)
    906967                   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) )
    909970                      sigma          = SQRT( e(kp,jp,ip) )
    910971
     
    9771038!--             (2010, Q. J. R. Meteorol. Soc.)
    9781039                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) )
    9811042                    sigma          = SQRT( e(kp,jp,ip) )
    9821043
Note: See TracChangeset for help on using the changeset viewer.