Changeset 4430


Ignore:
Timestamp:
Feb 27, 2020 6:02:20 PM (15 months ago)
Author:
suehring
Message:

Lagrangian particle model: Bugfix in logarithmic interpolation of near-ground particle speed (density was not considered); Revise CFL-check when SGS particle speeds are considered; .palm.iofiles: missing output file connection for child particle statistics

Location:
palm/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/.palm.iofiles

    r4428 r4430  
    6767DATA_1D_SP_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _sp        nc
    6868DATA_1D_TS_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _ts        nc
    69 DATA_1D_PTS_NETCDF         out:tr   *         $base_data/$run_identifier/OUTPUT         _pts       nc
     69DATA_1D_PTS_NETCDF*        out:tr   *         $base_data/$run_identifier/OUTPUT         _pts       nc
    7070DATA_2D_XY_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _xy        nc
    7171DATA_2D_XY_AV_NETCDF*      out:tr   *         $base_data/$run_identifier/OUTPUT         _av_xy     nc
  • palm/trunk/SOURCE/lagrangian_particle_model_mod.f90

    r4360 r4430  
    2525! -----------------
    2626! $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
    2735! Introduction of wall_flags_total_0, which currently sets bits based on static
    2836! topography information used in wall_flags_static_0
     
    126134
    127135    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
    131141 
    132142    USE averaging,                                                             &
     
    139149    USE control_parameters,                                                    &
    140150        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
     151               child_domain,                                                   &
    141152               cloud_droplets, constant_flux_layer, current_timestep_number,   &
    142153               dt_3d, dt_3d_reached, first_call_lpm, humidity,                 &
     
    830841!
    831842!--    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),                         &
    833844                  ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    834 !
    835845!--    Real volume of particles (with weighting), volume of particles
    836        ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
    837                      ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     846       ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
     847                  ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    838848    ENDIF
    839849
    840850
    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),                              &
    843853              w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    844854!
     
    23382348!------------------------------------------------------------------------------!
    23392349 SUBROUTINE particles_from_parent_to_child
    2340     IMPLICIT NONE
    23412350
    23422351    CALL pmcp_c_get_particle_from_parent                         ! Child actions
     
    23442353
    23452354    RETURN
     2355
    23462356 END SUBROUTINE particles_from_parent_to_child
    23472357
     
    23532363!------------------------------------------------------------------------------!
    23542364 SUBROUTINE particles_from_child_to_parent
    2355     IMPLICIT NONE
    23562365
    23572366    CALL pmcp_c_send_particle_to_parent                         ! Child actions
     
    23592368
    23602369    RETURN
     2370
    23612371 END SUBROUTINE particles_from_child_to_parent
    23622372 
     
    33153325    REAL(wp), DIMENSION(number_of_particles) ::  zv             !< z-position
    33163326
    3317     REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
     3327    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg          !< vector of Gaussian distributed random numbers
    33183328
    33193329    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
     
    34613471!
    34623472!--                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 ) &
    34653475                                        * d_z_p_z0
    34663476!
    34673477!--                Calculate LOG(z/z0) for exact particle height. Therefore,
    34683478!--                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))             &
    34733483                                      )
    34743484!
    34753485!--                Get friction velocity and momentum flux from new surface data
    34763486!--                types.
    3477                    IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3487                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
    34783488                        surf_def_h(0)%end_index(jp,ip) )  THEN
    34793489                      surf_start = surf_def_h(0)%start_index(jp,ip)
     
    34823492!--                   large particle speed.
    34833493                      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) <=                   &
    34863497                            surf_lsm_h%end_index(jp,ip) )  THEN
    34873498                      surf_start = surf_lsm_h%start_index(jp,ip)
    34883499                      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) <=                   &
    34913503                            surf_usm_h%end_index(jp,ip) )  THEN
    34923504                      surf_start = surf_usm_h%start_index(jp,ip)
    34933505                      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)
    34953508                   ENDIF
    34963509!
     
    35013514!--                as sensitivity studies revealed no significant effect of
    35023515!--                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 )        &
    35043517                               * log_z_z0_int - u_gtrans
    35053518                ENDIF
     
    35533566!--                topography particle on u-grid can be above surface-layer height,
    35543567!--                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 )    &
    35573570                                     * d_z_p_z0
    35583571!
    35593572!--                Calculate LOG(z/z0) for exact particle height. Therefore,
    35603573!--                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))             &
    35653578                                      )
    35663579!
    35673580!--                Get friction velocity and momentum flux from new surface data
    35683581!--                types.
    3569                    IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3582                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
    35703583                        surf_def_h(0)%end_index(jp,ip) )  THEN
    35713584                      surf_start = surf_def_h(0)%start_index(jp,ip)
     
    35743587!--                   large particle speed.
    35753588                      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) <=                   &
    35783592                            surf_lsm_h%end_index(jp,ip) )  THEN
    35793593                      surf_start = surf_lsm_h%start_index(jp,ip)
    35803594                      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) <=                   &
    35833598                            surf_usm_h%end_index(jp,ip) )  THEN
    35843599                      surf_start = surf_usm_h%start_index(jp,ip)
    35853600                      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)
    35873603                   ENDIF
    35883604!
     
    35933609!--                as sensitivity studies revealed no significant effect of
    35943610!--                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 )        &
    35963612                            * log_z_z0_int - v_gtrans
    35973613
     
    39013917                ENDIF
    39023918             ENDIF
     3919
    39033920             rvar1_temp(n) = particles(n)%rvar1
    39043921             rvar2_temp(n) = particles(n)%rvar2
     
    39183935                rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
    39193936                                          + 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
    3920 
    39213937             ELSE
    39223938!
     
    39663982!
    39673983!--    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.
    39703992       dz_temp = zw(kp)-zw(kp-1)
    39713993
    39723994       DO  nb = 0, 7
    39733995          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) ) )
    40094021             ENDIF
    4010 
    40114022!
    40124023!--          Update particle velocites
     
    42764287       ENDIF
    42774288!
    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.
    42794292       DO  n = particle_start, particle_end
    4280 !
    4281 !--       Increment the particle age and the total time that the particle
    4282 !--       has advanced within the particle timestep procedure
    42834293          particles(n)%age    = particles(n)%age    + dt_particle(n)
    42844294          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 )                     &
    42904315             dt_3d_reached_l = .FALSE.
    4291           ENDIF
    4292 
    42934316       ENDDO
    42944317    ENDDO
     
    70607083                      IF ( i < 0 )  THEN
    70617084!
    7062 !--                   Apply boundary condition along x
     7085!--                      Apply boundary condition along x
    70637086                         IF ( ibc_par_lr == 0 )  THEN
    70647087!
     
    74407463
    74417464       IF ( trsp_count_recv > 0 )  CALL lpm_add_particles_to_gridcell(rvsp(1:trsp_count_recv))
     7465
    74427466
    74437467       DEALLOCATE(rvsp)
     
    77427766             ENDIF
    77437767          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
    77457771          ENDIF
    77467772       ENDIF
     
    78637889!--             Note, check for CFL does not work at first particle timestep
    78647890!--             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                    THEN
     7891                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
    78737899                      WRITE( message_string, * )                               &
    78747900                      'Particle violated CFL-criterion: &particle with id ',   &
    78757901                      particles(n)%id, ' will be deleted!'   
    78767902                      CALL message( 'lpm_check_cfl', 'PA0475', 0, 1, -1, 6, 0 )
     7903
    78777904                      particles(n)%particle_mask= .FALSE.
    78787905                   ENDIF
Note: See TracChangeset for help on using the changeset viewer.