Changeset 3937 for palm/trunk/SOURCE


Ignore:
Timestamp:
Apr 29, 2019 3:09:07 PM (6 years ago)
Author:
suehring
Message:

Move initialization of STG behind initialization of offline nesting; Bugfix in STG in case of very early restart; calculation of scaling parameters used for parametrization of synthetic turbulence profiles improved; in offline nesting, set boundary value at upper-left and upper-south grid points for u- and v-component, respectively

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r3900 r3937  
    2525! -----------------
    2626! $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
    2732! Fix problem with LOD = 2 initialization
    2833!
     
    666671
    667672    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
    670674
    671675    USE surface_layer_fluxes_mod,                                              &
     
    13051309!--       fluxes, etc.
    13061310          CALL init_surfaces
    1307 !
    1308 !--       Initialize turbulence generator
    1309           IF( use_syn_turb_gen )  CALL stg_init
    13101311
    13111312          CALL location_message( 'initializing with INIFOR', 'finished' )
     
    13751376!--       fluxes, etc.
    13761377          CALL init_surfaces
    1377 !
    1378 !--       Initialize synthetic turbulence generator if required
    1379           IF( use_syn_turb_gen )  CALL stg_init
    13801378
    13811379          CALL location_message( 'initializing with 1D model profiles', 'finished' )
     
    14461444!--       fluxes, etc.
    14471445          CALL init_surfaces
    1448 !
    1449 !--       Initialize synthetic turbulence generator if required
    1450           IF( use_syn_turb_gen )  CALL stg_init
    14511446         
    14521447          CALL location_message( 'initializing with constant profiles', 'finished' )
     
    18681863       ENDIF
    18691864       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_init
    18741865
    18751866       CALL location_message( 'initializing in case of restart / cyclic_fill', 'finished' )
     
    23152306    IF ( agents_active )  CALL mas_init
    23162307!
    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
    23252310!
    23262311!-- Initializing actions for all other modules
  • palm/trunk/SOURCE/nesting_offl_mod.f90

    r3891 r3937  
    2525! -----------------
    2626! $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
    2731! Bugfix, do not overwrite lateral and top boundary data in case of restart
    2832! runs.
     
    633637          ENDDO
    634638       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)
    635644
    636645       DO  i = nxl, nxr
     
    644653          ENDDO
    645654       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,:)
    646660
    647661       DO  i = nxl, nxr
     
    961975       topo_max     = topo_max_l
    962976#endif
    963 
    964977       zi_ribulk = MAX( zi_ribulk, topo_max )
    965 
     978       
    966979    END SUBROUTINE calc_zi
    967980   
  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r3909 r3937  
    2525! -----------------
    2626! $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
    2733! Minor bugfix for last commit
    2834!
     
    270276    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
    271277    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
    272     REAL(wp) ::  mc_factor                !< mass flux correction factor
     278    REAL(wp) ::  mc_factor = 1.0_wp       !< mass flux correction factor
    273279    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
    274280    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
     
    886892          DO  j = nysg, nyng
    887893             DO  k = nzb, nzt+1
    888 
    889894                IF  ( a11(k) .NE. 0._wp ) THEN
    890895                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
     
    17791784!--       zero.
    17801785          r11(k) = scale_us**2 * (                                             &
    1781                       MERGE( 0.35_wp * (                                       &
     1786                      MERGE( 0.35_wp * ABS(                                    &
    17821787                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
    1783                                        )**( 2.0_wp / 3.0_wp ),                 &
     1788                                          )**( 2.0_wp / 3.0_wp ),              &
    17841789                             0.0_wp,                                           &
    17851790                             scale_l < 0.0_wp )                                &
     
    20522057!------------------------------------------------------------------------------!
    20532058 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         
    20592062             
    20602063    USE indices,                                                               &
     
    20742077
    20752078    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
    20762082    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
    20772083    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
     
    20792085    REAL(wp) ::  v_int                  !< v-component
    20802086    REAL(wp) ::  w_convective           !< convective velocity scale
    2081     REAL(wp) ::  z0_mean                !< mean roughness length
    2082     REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
    20832087   
    20842088!
    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
    20932095    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)
    20962103    ENDDO   
    20972104    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)
    21002112    ENDDO
    21012113    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)
    21042121    ENDDO
    21052122   
    21062123#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,    &
    21082125                        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,    &
    21102131                        comm2d, ierr )
    21112132#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
    21142137#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 )   
    21992143!
    22002144!-- Compute mean convective velocity scale. Note, in case the mean heat flux
    22012145!-- is negative, set convective velocity scale to zero.
    22022146    IF ( shf_mean > 0.0_wp )  THEN
    2203        w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
     2147       w_convective = ( g * shf_mean * zi_ribulk / pt_surf_mean )**( 1.0_wp / 3.0_wp )
    22042148    ELSE
    22052149       w_convective = 0.0_wp
Note: See TracChangeset for help on using the changeset viewer.