Changeset 4588 for palm/trunk


Ignore:
Timestamp:
Jul 6, 2020 11:06:02 AM (4 years ago)
Author:
suehring
Message:

Simplify particle-speed interpolation in logarithmic layer

File:
1 edited

Legend:

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

    r4585 r4588  
    2525! -----------------
    2626! $Id$
     27! Simplify particle-speed interpolation in logarithmic layer
     28!
     29! 4585 2020-06-30 15:05:20Z suehring
    2730! Limit logarithmically interpolated particle speed to the velocity component
    2831! at the first prognostic grid point (since no stability corrected interpolation
     
    345348    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  radius = 9999999.9_wp         !< namelist parameter (see documentation)
    346349
    347     REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0) 
     350    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)
    348351
    349352    INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !<
     
    34883491    REAL(wp) ::  unext              !< calculated particle u-velocity of corrector step
    34893492    REAL(wp) ::  us_int             !< friction velocity at particle grid box
    3490     REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
    34913493    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
    34923494    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
    3493     REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
    34943495    REAL(wp) ::  vnext              !< calculated particle v-velocity of corrector step
    34953496    REAL(wp) ::  vv_int             !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection
     
    36913692                                      )
    36923693!
    3693 !--                Get friction velocity and momentum flux from new surface data
    3694 !--                types.
    3695                    IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
    3696                         surf_def_h(0)%end_index(jp,ip) )  THEN
    3697                       surf_start = surf_def_h(0)%start_index(jp,ip)
    3698 !--                   Limit friction velocity. In narrow canyons or holes the
    3699 !--                   friction velocity can become very small, resulting in a too
    3700 !--                   large particle speed.
    3701                       us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
    3702                       usws_int  = surf_def_h(0)%usws(surf_start)               &
    3703                                 * drho_air_zw(k_wall)
    3704                    ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
    3705                             surf_lsm_h%end_index(jp,ip) )  THEN
    3706                       surf_start = surf_lsm_h%start_index(jp,ip)
    3707                       us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
    3708                       usws_int  = surf_lsm_h%usws(surf_start)                  &
    3709                                 * drho_air_zw(k_wall)
    3710                    ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
    3711                             surf_usm_h%end_index(jp,ip) )  THEN
    3712                       surf_start = surf_usm_h%start_index(jp,ip)
    3713                       us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
    3714                       usws_int  = surf_usm_h%usws(surf_start)                  &
    3715                                 * drho_air_zw(k_wall)
    3716                    ENDIF
    3717 !
    3718 !--                Neutral solution is applied for all situations, e.g. also for
     3694!--                Compute u*-portion for u-component based on mean roughness.
     3695!--                Note, neutral solution is applied for all situations, e.g. also for
    37193696!--                unstable and stable situations. Even though this is not exact
    37203697!--                this saves a lot of CPU time since several calls of intrinsic
    37213698!--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    37223699!--                as sensitivity studies revealed no significant effect of
    3723 !--                using the neutral solution also for un/stable situations.
    3724 !--                However, as the stability correction is missing the particle
    3725 !--                speed may be overestimated compared to the velocity component
    3726 !--                at the first prognostic grid point. Hence, the particle speed
    3727 !--                is limited.
    3728                    u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )        &
    3729                                * log_z_z0_int - u_gtrans
    3730 
    3731                    u_int(n) = SIGN( MIN( ABS( u_int(n) ), ABS( u(k_wall+1,jp,ip) ) ), u_int(n) )
     3700!--                using the neutral solution also for un/stable situations. Based on the u*
     3701!--                recalculate the velocity at height z_particle. Since the analytical solution
     3702!--                only yields absolute values, include the sign using the intrinsic SIGN function.
     3703                   us_int   = kappa * 0.5_wp * ABS( u(k_wall+1,jp,ip) + u(k_wall+1,jp,ip+1) ) /    &
     3704                              log_z_z0(number_of_sublayers)
     3705                   u_int(n) = SIGN( 1.0_wp, u(k_wall+1,jp,ip) + u(k_wall+1,jp,ip+1) ) *            &
     3706                              log_z_z0_int * us_int / kappa - u_gtrans
    37323707
    37333708                ENDIF
     
    37773752                ELSE
    37783753!
    3779 !--                Determine the sublayer. Further used as index. Please note,
    3780 !--                logarithmus can not be reused from above, as in in case of
    3781 !--                topography particle on u-grid can be above surface-layer height,
    3782 !--                whereas it can be below on v-grid.
     3754!--                Determine the sublayer. Further used as index.
    37833755                   height_p = ( zv(n) - zw(k_wall) - z0_av_global )            &
    3784                                      * REAL( number_of_sublayers, KIND=wp )    &
    3785                                      * d_z_p_z0
     3756                                        * REAL( number_of_sublayers, KIND=wp ) &
     3757                                        * d_z_p_z0
    37863758!
    37873759!--                Calculate LOG(z/z0) for exact particle height. Therefore,
     
    37933765                                      )
    37943766!
    3795 !--                Get friction velocity and momentum flux from new surface data
    3796 !--                types.
    3797                    IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
    3798                         surf_def_h(0)%end_index(jp,ip) )  THEN
    3799                       surf_start = surf_def_h(0)%start_index(jp,ip)
    3800 !--                   Limit friction velocity. In narrow canyons or holes the
    3801 !--                   friction velocity can become very small, resulting in a too
    3802 !--                   large particle speed.
    3803                       us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
    3804                       vsws_int  = surf_def_h(0)%vsws(surf_start)               &
    3805                                 * drho_air_zw(k_wall)
    3806                    ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
    3807                             surf_lsm_h%end_index(jp,ip) )  THEN
    3808                       surf_start = surf_lsm_h%start_index(jp,ip)
    3809                       us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
    3810                       vsws_int  = surf_lsm_h%vsws(surf_start)                  &
    3811                                 * drho_air_zw(k_wall)
    3812                    ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
    3813                             surf_usm_h%end_index(jp,ip) )  THEN
    3814                       surf_start = surf_usm_h%start_index(jp,ip)
    3815                       us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
    3816                       vsws_int  = surf_usm_h%vsws(surf_start)                  &
    3817                                 * drho_air_zw(k_wall)
    3818                    ENDIF
    3819 !
    3820 !--                Neutral solution is applied for all situations, e.g. also for
     3767!--                Compute u*-portion for v-component based on mean roughness.
     3768!--                Note, neutral solution is applied for all situations, e.g. also for
    38213769!--                unstable and stable situations. Even though this is not exact
    38223770!--                this saves a lot of CPU time since several calls of intrinsic
    38233771!--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    38243772!--                as sensitivity studies revealed no significant effect of
    3825 !--                using the neutral solution also for un/stable situations.
    3826 !--                However, as the stability correction is missing the particle
    3827 !--                speed may be overestimated compared to the velocity component
    3828 !--                at the first prognostic grid point. Hence, the particle speed
    3829 !--                is limited.
    3830                    v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )        &
    3831                             * log_z_z0_int - v_gtrans
    3832 
    3833                    v_int(n) = SIGN( MIN( ABS( v_int(n) ), ABS( v(k_wall+1,jp,ip) ) ), v_int(n) )
     3773!--                using the neutral solution also for un/stable situations. Based on the u*
     3774!--                recalculate the velocity at height z_particle. Since the analytical solution
     3775!--                only yields absolute values, include the sign using the intrinsic SIGN function.
     3776                   us_int   = kappa * 0.5_wp * ABS( v(k_wall+1,jp,ip) + v(k_wall+1,jp+1,ip) ) /    &
     3777                              log_z_z0(number_of_sublayers)
     3778                   v_int(n) = SIGN( 1.0_wp, v(k_wall+1,jp,ip) + v(k_wall+1,jp+1,ip) ) *            &
     3779                              log_z_z0_int * us_int / kappa - v_gtrans
    38343780
    38353781                ENDIF
Note: See TracChangeset for help on using the changeset viewer.