Changeset 3937 for palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
 Timestamp:
 Apr 29, 2019 3:09:07 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r3909 r3937 25 25 !  26 26 ! $Id$ 27 ! Minor bugfix in case of a very early restart where mc_factor is sill not 28 ! present. 29 ! Some modification and fixing of potential bugs in the calculation of scaling 30 ! parameters used for synthetic turbulence parametrization. 31 ! 32 ! 3909 20190417 09:13:25Z suehring 27 33 ! Minor bugfix for last commit 28 34 ! … … 270 276 REAL(wp) :: dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics 271 277 REAL(wp) :: dt_stg_call = 5.0_wp !< time interval for calling synthetic turbulence generator 272 REAL(wp) :: mc_factor 278 REAL(wp) :: mc_factor = 1.0_wp !< mass flux correction factor 273 279 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization  Obukhov length 274 280 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization  friction velocity … … 886 892 DO j = nysg, nyng 887 893 DO k = nzb, nzt+1 888 889 894 IF ( a11(k) .NE. 0._wp ) THEN 890 895 fu_yz(k,j) = ( u(k,j,i) / mc_factor  u_init(k) ) / a11(k) … … 1779 1784 ! zero. 1780 1785 r11(k) = scale_us**2 * ( & 1781 MERGE( 0.35_wp * (&1786 MERGE( 0.35_wp * ABS( & 1782 1787  zi_ribulk / ( kappa * scale_l  10E4_wp ) & 1783 )**( 2.0_wp / 3.0_wp ),&1788 )**( 2.0_wp / 3.0_wp ), & 1784 1789 0.0_wp, & 1785 1790 scale_l < 0.0_wp ) & … … 2052 2057 !! 2053 2058 SUBROUTINE calc_scaling_variables 2054 2055 2056 USE control_parameters, & 2057 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, & 2058 pt_surface 2059 2060 USE arrays_3d, & 2061 ONLY: drho_air 2059 2062 2060 2063 USE indices, & … … 2074 2077 2075 2078 REAL(wp) :: friction_vel_l !< mean friction veloctiy on subdomain 2079 REAL(wp) :: pt_surf_mean !< mean near surface temperature (at 1st grid point) 2080 REAL(wp) :: pt_surf_mean_l !< mean near surface temperature (at 1st grid point) on subdomain 2081 REAL(wp) :: scale_l_l !< mean Obukhov lenght on subdomain 2076 2082 REAL(wp) :: shf_mean !< mean surface sensible heat flux 2077 2083 REAL(wp) :: shf_mean_l !< mean surface sensible heat flux on subdomain … … 2079 2085 REAL(wp) :: v_int !< vcomponent 2080 2086 REAL(wp) :: w_convective !< convective velocity scale 2081 REAL(wp) :: z0_mean !< mean roughness length2082 REAL(wp) :: z0_mean_l !< mean roughness length on subdomain2083 2087 2084 2088 ! 2085 ! Mean friction velocity and velocity scale. Therefore, 2086 ! precalculate mean roughness length and surface sensible heat flux 2087 ! in the model domain, which are further used to estimate friction 2088 ! velocity and velocity scale. Note, for z0 linear averaging is applied, 2089 ! even though this is known to unestimate the effective roughness. 2090 ! This need to be revised later. 2091 z0_mean_l = 0.0_wp 2092 shf_mean_l = 0.0_wp 2089 ! Calculate mean friction velocity, velocity scale, heat flux and 2090 ! nearsurface temperature in the model domain. 2091 pt_surf_mean_l = 0.0_wp 2092 shf_mean_l = 0.0_wp 2093 scale_l_l = 0.0_wp 2094 friction_vel_l = 0.0_wp 2093 2095 DO m = 1, surf_def_h(0)%ns 2094 z0_mean_l = z0_mean_l + surf_def_h(0)%z0(m) 2095 shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m) 2096 i = surf_def_h(0)%i(m) 2097 j = surf_def_h(0)%j(m) 2098 k = surf_def_h(0)%k(m) 2099 friction_vel_l = friction_vel_l + surf_def_h(0)%us(m) 2100 shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m) * drho_air(k) 2101 scale_l_l = scale_l_l + surf_def_h(0)%ol(m) 2102 pt_surf_mean_l = pt_surf_mean_l + pt(k,j,i) 2096 2103 ENDDO 2097 2104 DO m = 1, surf_lsm_h%ns 2098 z0_mean_l = z0_mean_l + surf_lsm_h%z0(m) 2099 shf_mean_l = shf_mean_l + surf_lsm_h%shf(m) 2105 i = surf_lsm_h%i(m) 2106 j = surf_lsm_h%j(m) 2107 k = surf_lsm_h%k(m) 2108 friction_vel_l = friction_vel_l + surf_lsm_h%us(m) 2109 shf_mean_l = shf_mean_l + surf_lsm_h%shf(m) * drho_air(k) 2110 scale_l_l = scale_l_l + surf_lsm_h%ol(m) 2111 pt_surf_mean_l = pt_surf_mean_l + pt(k,j,i) 2100 2112 ENDDO 2101 2113 DO m = 1, surf_usm_h%ns 2102 z0_mean_l = z0_mean_l + surf_usm_h%z0(m) 2103 shf_mean_l = shf_mean_l + surf_usm_h%shf(m) 2114 i = surf_usm_h%i(m) 2115 j = surf_usm_h%j(m) 2116 k = surf_usm_h%k(m) 2117 friction_vel_l = friction_vel_l + surf_usm_h%us(m) 2118 shf_mean_l = shf_mean_l + surf_usm_h%shf(m) * drho_air(k) 2119 scale_l_l = scale_l_l + surf_usm_h%ol(m) 2120 pt_surf_mean_l = pt_surf_mean_l + pt(k,j,i) 2104 2121 ENDDO 2105 2122 2106 2123 #if defined( __parallel ) 2107 CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,&2124 CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM, & 2108 2125 comm2d, ierr ) 2109 CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM, & 2126 CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM, & 2127 comm2d, ierr ) 2128 CALL MPI_ALLREDUCE( scale_l_l, scale_l, 1, MPI_REAL, MPI_SUM, & 2129 comm2d, ierr ) 2130 CALL MPI_ALLREDUCE( pt_surf_mean_l, pt_surf_mean, 1, MPI_REAL, MPI_SUM, & 2110 2131 comm2d, ierr ) 2111 2132 #else 2112 z0_mean = z0_mean_l 2113 shf_mean = shf_mean_l 2133 scale_us = friction_vel_l 2134 shf_mean = shf_mean_l 2135 scale_l = scale_l_l 2136 pt_surf_mean = pt_surf_mean_l 2114 2137 #endif 2115 z0_mean = z0_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2116 shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2117 2118 ! 2119 ! Note, Inifor does not use logarithmic interpolation of the 2120 ! velocity components near the ground, so that nearsurface 2121 ! wind speed and thus the friction velocity is overestimated. 2122 ! However, friction velocity is used for turbulence 2123 ! parametrization, so that more physically meaningful values are important. 2124 ! Hence, derive friction velocity from wind speed at a reference height, 2125 ! which is 10 m, according to the height of the 1st vertical grid level 2126 ! in the COSMO level. However, in case of topography that is higher than 2127 ! the reference level, the k index is determined from the 1st vertical 2128 ! PALM grid level instead. 2129 ! For a first guess use 20 m, which is in the range of the first 2130 ! COSMO vertical level. 2131 k_ref = MINLOC( ABS( zu  cosmo_ref ), DIM = 1 )  1 2132 2133 friction_vel_l = 0.0_wp 2134 IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN 2135 2136 i = MERGE( 1, nxr + 1, bc_dirichlet_l ) 2137 2138 DO j = nys, nyn 2139 ! 2140 ! Determine the k index and topography top index 2141 k_topo = MAX( get_topography_top_index_ji( j, i, 'u' ), & 2142 get_topography_top_index_ji( j, i, 'v' ) ) 2143 k = MAX( k_ref, k_topo + 1 ) 2144 ! 2145 ! Note, in u and v component the imposed perturbations 2146 ! from the STG are already included. Check whether this 2147 ! makes any difference compared to using the puremean 2148 ! inflow profiles. 2149 u_int = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l ) 2150 v_int = v(k,j,i) 2151 ! 2152 ! Calculate friction velocity and sumup. Therefore, assume 2153 ! neutral condtions. 2154 friction_vel_l = friction_vel_l + kappa * & 2155 SQRT( u_int * u_int + v_int * v_int ) / & 2156 LOG( ( zu(k)  zu(k_topo) ) / z0_mean ) 2157 2158 ENDDO 2159 2160 ENDIF 2161 2162 IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN 2163 2164 j = MERGE( 1, nyn + 1, bc_dirichlet_s ) 2165 2166 DO i = nxl, nxr 2167 2168 k_topo = MAX( get_topography_top_index_ji( j, i, 'u' ), & 2169 get_topography_top_index_ji( j, i, 'v' ) ) 2170 k = MAX( k_ref, k_topo + 1 ) 2171 2172 u_int = u(k,j,i) 2173 v_int = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s ) 2174 2175 friction_vel_l = friction_vel_l + kappa * & 2176 SQRT( u_int * u_int + v_int * v_int ) / & 2177 LOG( ( zu(k)  zu(k_topo) ) / z0_mean ) 2178 2179 ENDDO 2180 2181 ENDIF 2182 2183 #if defined( __parallel ) 2184 CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM, & 2185 comm2d, ierr ) 2186 #else 2187 scale_us = friction_vel_l 2188 #endif 2189 scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp ) 2190 2191 ! 2192 ! Compute mean Obukhov length 2193 IF ( shf_mean > 0.0_wp ) THEN 2194 scale_l =  pt_surface / ( kappa * g ) * scale_us**3 / shf_mean 2195 ELSE 2196 scale_l = 0.0_wp 2197 ENDIF 2198 2138 2139 scale_us = scale_us / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2140 shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2141 scale_l = scale_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2142 pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2199 2143 ! 2200 2144 ! Compute mean convective velocity scale. Note, in case the mean heat flux 2201 2145 ! is negative, set convective velocity scale to zero. 2202 2146 IF ( shf_mean > 0.0_wp ) THEN 2203 w_convective = ( g * shf_mean * zi_ribulk / pt_surf ace)**( 1.0_wp / 3.0_wp )2147 w_convective = ( g * shf_mean * zi_ribulk / pt_surf_mean )**( 1.0_wp / 3.0_wp ) 2204 2148 ELSE 2205 2149 w_convective = 0.0_wp
Note: See TracChangeset
for help on using the changeset viewer.