Changeset 4588 for palm/trunk

Timestamp:

Jul 6, 2020 11:06:02 AM (4 years ago)

Author:

suehring

Message:

Simplify particle-speed interpolation in logarithmic layer

File:

• palm/trunk/SOURCE/lagrangian_particle_model_mod.f90 (modified) (6 diffs)

palm/trunk/SOURCE/lagrangian_particle_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Simplify particle-speed interpolation in logarithmic layer
+! 
+! 4585 2020-06-30 15:05:20Z suehring
 ! Limit logarithmically interpolated particle speed to the velocity component
 ! at the first prognostic grid point (since no stability corrected interpolation
@@ -345 +348 @@
     REAL(wp), DIMENSION(max_number_of_particle_groups) ::  radius = 9999999.9_wp         !< namelist parameter (see documentation)
 
-    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)  
+    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)
 
     INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !<
@@ -3488 +3491 @@
     REAL(wp) ::  unext              !< calculated particle u-velocity of corrector step
     REAL(wp) ::  us_int             !< friction velocity at particle grid box
-    REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
     REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
     REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
-    REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
     REAL(wp) ::  vnext              !< calculated particle v-velocity of corrector step
     REAL(wp) ::  vv_int             !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection 
@@ -3691 +3692 @@
                                       )
 !
-!--                Get friction velocity and momentum flux from new surface data
-!--                types.
-                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
-                        surf_def_h(0)%end_index(jp,ip) )  THEN
-                      surf_start = surf_def_h(0)%start_index(jp,ip)
-!--                   Limit friction velocity. In narrow canyons or holes the
-!--                   friction velocity can become very small, resulting in a too
-!--                   large particle speed.
-                      us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
-                      usws_int  = surf_def_h(0)%usws(surf_start)               &
-                                * drho_air_zw(k_wall)
-                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
-                            surf_lsm_h%end_index(jp,ip) )  THEN
-                      surf_start = surf_lsm_h%start_index(jp,ip)
-                      us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
-                      usws_int  = surf_lsm_h%usws(surf_start)                  &
-                                * drho_air_zw(k_wall)
-                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
-                            surf_usm_h%end_index(jp,ip) )  THEN
-                      surf_start = surf_usm_h%start_index(jp,ip)
-                      us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
-                      usws_int  = surf_usm_h%usws(surf_start)                  &
-                                * drho_air_zw(k_wall)
-                   ENDIF
-!
-!--                Neutral solution is applied for all situations, e.g. also for
+!--                Compute u*-portion for u-component based on mean roughness.
+!--                Note, neutral solution is applied for all situations, e.g. also for
 !--                unstable and stable situations. Even though this is not exact
 !--                this saves a lot of CPU time since several calls of intrinsic
 !--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
 !--                as sensitivity studies revealed no significant effect of
-!--                using the neutral solution also for un/stable situations.
-!--                However, as the stability correction is missing the particle
-!--                speed may be overestimated compared to the velocity component
-!--                at the first prognostic grid point. Hence, the particle speed
-!--                is limited.
-                   u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )        &
-                               * log_z_z0_int - u_gtrans
-
-                   u_int(n) = SIGN( MIN( ABS( u_int(n) ), ABS( u(k_wall+1,jp,ip) ) ), u_int(n) )
+!--                using the neutral solution also for un/stable situations. Based on the u*
+!--                recalculate the velocity at height z_particle. Since the analytical solution
+!--                only yields absolute values, include the sign using the intrinsic SIGN function.
+                   us_int   = kappa * 0.5_wp * ABS( u(k_wall+1,jp,ip) + u(k_wall+1,jp,ip+1) ) /    &
+                              log_z_z0(number_of_sublayers)
+                   u_int(n) = SIGN( 1.0_wp, u(k_wall+1,jp,ip) + u(k_wall+1,jp,ip+1) ) *            &
+                              log_z_z0_int * us_int / kappa - u_gtrans
 
                 ENDIF
@@ -3777 +3752 @@
                 ELSE
 !
-!--                Determine the sublayer. Further used as index. Please note,
-!--                logarithmus can not be reused from above, as in in case of
-!--                topography particle on u-grid can be above surface-layer height,
-!--                whereas it can be below on v-grid.
+!--                Determine the sublayer. Further used as index.
                    height_p = ( zv(n) - zw(k_wall) - z0_av_global )            &
-                                     * REAL( number_of_sublayers, KIND=wp )    &
-                                     * d_z_p_z0
+                                        * REAL( number_of_sublayers, KIND=wp ) &
+                                        * d_z_p_z0
 !
 !--                Calculate LOG(z/z0) for exact particle height. Therefore,
@@ -3793 +3765 @@
                                       )
 !
-!--                Get friction velocity and momentum flux from new surface data
-!--                types.
-                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
-                        surf_def_h(0)%end_index(jp,ip) )  THEN
-                      surf_start = surf_def_h(0)%start_index(jp,ip)
-!--                   Limit friction velocity. In narrow canyons or holes the
-!--                   friction velocity can become very small, resulting in a too
-!--                   large particle speed.
-                      us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_def_h(0)%vsws(surf_start)               &
-                                * drho_air_zw(k_wall)
-                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
-                            surf_lsm_h%end_index(jp,ip) )  THEN
-                      surf_start = surf_lsm_h%start_index(jp,ip)
-                      us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_lsm_h%vsws(surf_start)                  &
-                                * drho_air_zw(k_wall)
-                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
-                            surf_usm_h%end_index(jp,ip) )  THEN
-                      surf_start = surf_usm_h%start_index(jp,ip)
-                      us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_usm_h%vsws(surf_start)                  &
-                                * drho_air_zw(k_wall)
-                   ENDIF
-!
-!--                Neutral solution is applied for all situations, e.g. also for
+!--                Compute u*-portion for v-component based on mean roughness.
+!--                Note, neutral solution is applied for all situations, e.g. also for
 !--                unstable and stable situations. Even though this is not exact
 !--                this saves a lot of CPU time since several calls of intrinsic
 !--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
 !--                as sensitivity studies revealed no significant effect of
-!--                using the neutral solution also for un/stable situations.
-!--                However, as the stability correction is missing the particle
-!--                speed may be overestimated compared to the velocity component
-!--                at the first prognostic grid point. Hence, the particle speed
-!--                is limited.
-                   v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )        &
-                            * log_z_z0_int - v_gtrans
-
-                   v_int(n) = SIGN( MIN( ABS( v_int(n) ), ABS( v(k_wall+1,jp,ip) ) ), v_int(n) )
+!--                using the neutral solution also for un/stable situations. Based on the u*
+!--                recalculate the velocity at height z_particle. Since the analytical solution
+!--                only yields absolute values, include the sign using the intrinsic SIGN function.
+                   us_int   = kappa * 0.5_wp * ABS( v(k_wall+1,jp,ip) + v(k_wall+1,jp+1,ip) ) /    &
+                              log_z_z0(number_of_sublayers)
+                   v_int(n) = SIGN( 1.0_wp, v(k_wall+1,jp,ip) + v(k_wall+1,jp+1,ip) ) *            &
+                              log_z_z0_int * us_int / kappa - v_gtrans
 
                 ENDIF