- Timestamp:
- Apr 29, 2019 3:09:07 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
init_3d_model.f90 (modified) (7 diffs)
-
nesting_offl_mod.f90 (modified) (4 diffs)
-
synthetic_turbulence_generator_mod.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_3d_model.f90
r3900 r3937 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Move initialization of synthetic turbulence generator behind initialization 28 ! of offline nesting. Remove call for stg_adjust, as this is now already done 29 ! in stg_init. 30 ! 31 ! 3900 2019-04-16 15:17:43Z suehring 27 32 ! Fix problem with LOD = 2 initialization 28 33 ! … … 666 671 667 672 USE synthetic_turbulence_generator_mod, & 668 ONLY: parametrize_inflow_turbulence, stg_adjust, stg_init, & 669 use_syn_turb_gen 673 ONLY: parametrize_inflow_turbulence, stg_init, use_syn_turb_gen 670 674 671 675 USE surface_layer_fluxes_mod, & … … 1305 1309 !-- fluxes, etc. 1306 1310 CALL init_surfaces 1307 !1308 !-- Initialize turbulence generator1309 IF( use_syn_turb_gen ) CALL stg_init1310 1311 1311 1312 CALL location_message( 'initializing with INIFOR', 'finished' ) … … 1375 1376 !-- fluxes, etc. 1376 1377 CALL init_surfaces 1377 !1378 !-- Initialize synthetic turbulence generator if required1379 IF( use_syn_turb_gen ) CALL stg_init1380 1378 1381 1379 CALL location_message( 'initializing with 1D model profiles', 'finished' ) … … 1446 1444 !-- fluxes, etc. 1447 1445 CALL init_surfaces 1448 !1449 !-- Initialize synthetic turbulence generator if required1450 IF( use_syn_turb_gen ) CALL stg_init1451 1446 1452 1447 CALL location_message( 'initializing with constant profiles', 'finished' ) … … 1868 1863 ENDIF 1869 1864 IF ( passive_scalar ) ts_m = 0.0_wp 1870 !1871 !-- Initialize synthetic turbulence generator in case of restart.1872 IF ( TRIM( initializing_actions ) == 'read_restart_data' .AND. &1873 use_syn_turb_gen ) CALL stg_init1874 1865 1875 1866 CALL location_message( 'initializing in case of restart / cyclic_fill', 'finished' ) … … 2315 2306 IF ( agents_active ) CALL mas_init 2316 2307 ! 2317 !-- In case the synthetic turbulence generator does not have any information 2318 !-- about the inflow turbulence, these information will be parametrized 2319 !-- depending on the initial atmospheric conditions and surface properties. 2320 !-- Please note, within pre-determined time intervals these turbulence 2321 !-- information can be updated if desired. 2322 IF ( use_syn_turb_gen .AND. parametrize_inflow_turbulence ) THEN 2323 CALL stg_adjust 2324 ENDIF 2308 !-- Initialization of synthetic turbulence generator 2309 IF ( use_syn_turb_gen ) CALL stg_init 2325 2310 ! 2326 2311 !-- Initializing actions for all other modules -
palm/trunk/SOURCE/nesting_offl_mod.f90
r3891 r3937 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Set boundary conditon on upper-left and upper-south grid point for the u- and 28 ! v-component, respectively. 29 ! 30 ! 3891 2019-04-12 17:52:01Z suehring 27 31 ! Bugfix, do not overwrite lateral and top boundary data in case of restart 28 32 ! runs. … … 633 637 ENDDO 634 638 ENDDO 639 ! 640 !-- For left boundary set boundary condition for u-component also at top 641 !-- grid point. 642 !-- Note, this has no effect on the numeric solution, only for data output. 643 IF ( bc_dirichlet_l ) u(nzt+1,:,nxl) = u(nzt+1,:,nxlu) 635 644 636 645 DO i = nxl, nxr … … 644 653 ENDDO 645 654 ENDDO 655 ! 656 !-- For south boundary set boundary condition for v-component also at top 657 !-- grid point. 658 !-- Note, this has no effect on the numeric solution, only for data output. 659 IF ( bc_dirichlet_s ) v(nzt+1,nys,:) = v(nzt+1,nysv,:) 646 660 647 661 DO i = nxl, nxr … … 961 975 topo_max = topo_max_l 962 976 #endif 963 964 977 zi_ribulk = MAX( zi_ribulk, topo_max ) 965 978 966 979 END SUBROUTINE calc_zi 967 980 -
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 2019-04-17 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 - 10E-4_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 !< v-component 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 !-- pre-calculate 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 !-- near-surface 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 near-surface 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 pure-mean 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 sum-up. 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.
