Changeset 3347


Ignore:
Timestamp:
Oct 15, 2018 2:21:08 PM (3 years ago)
Author:
suehring
Message:

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

Location:
palm/trunk/SOURCE
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r3337 r3347  
    2525# -----------------
    2626# $Id$
     27# Add module for offline nesting;
     28# Add surface_mod to synthetic-turbulence generator;
     29# Bugfix, missing dependency for turbulence generator in init_3d_model;
     30# Some formatting ajdustments
     31#
     32# 3343 2018-10-15 10:38:52Z suehring
    2733# (from branch resler)
    2834# Add biometeorology
     
    566572        netcdf_data_input_mod.f90 \
    567573        netcdf_interface_mod.f90 \
     574        nesting_offl_mod.f90 \
    568575        ocean_mod.f90 \
    569576        outflow_turbulence.f90 \
     
    10051012        modules.o \
    10061013        netcdf_interface_mod.o \
     1014        nesting_offl_mod.o \
    10071015        ocean_mod.o \
    10081016        plant_canopy_model_mod.o \
     
    10361044        netcdf_data_input_mod.o \
    10371045        netcdf_interface_mod.o \
     1046        nesting_offl_mod.o \
    10381047        ocean_mod.o \
    10391048        plant_canopy_model_mod.o \
     
    10451054        surface_layer_fluxes_mod.o \
    10461055        surface_mod.o \
     1056        synthetic_turbulence_generator_mod.o \
    10471057        turbulence_closure_mod.o \
    10481058        urban_surface_mod.o \
     
    11121122        mod_kinds.o \
    11131123        modules.o \
    1114         netcdf_data_input_mod.o \
    11151124        surface_mod.o
    11161125local_stop.o: \
     
    12881297        urban_surface_mod.o \
    12891298        uv_exposure_model_mod.o
     1299nesting_offl_mod.o: \
     1300        cpulog_mod.o \
     1301        mod_kinds.o \
     1302        modules.o \
     1303        netcdf_data_input_mod.o \
     1304        surface_mod.o
    12901305ocean_mod.o: \
    12911306        advec_s_pw.o \
     
    13351350        modules.o \
    13361351        multi_agent_system_mod.o \
     1352        nesting_offl_mod.o \
    13371353        netcdf_interface_mod.o \
    13381354        ocean_mod.o \
     
    15771593        mod_kinds.o \
    15781594        modules.o \
    1579         pmc_interface_mod.o
     1595        nesting_offl_mod.o \
     1596        pmc_interface_mod.o \
     1597        surface_mod.o
    15801598temperton_fft_mod.o: \
    15811599        mod_kinds.o \
     
    16001618        modules.o \
    16011619        multi_agent_system_mod.o \
     1620        nesting_offl_mod.o \
    16021621        ocean_mod.o \
    16031622        pmc_interface_mod.o \
  • palm/trunk/SOURCE/boundary_conds.f90

    r3294 r3347  
    2525! -----------------
    2626! $Id$
     27! Bugfix in setting boundary conditions in offline nesting
     28!
     29! 3341 2018-10-15 10:31:27Z suehring
    2730! changes concerning modularization of ocean option
    2831!
     
    293296!
    294297!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
    295     IF (  .NOT.  child_domain  .AND.                                           &
     298    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
    296299                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
    297300       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
     
    729732!-- in case of nest boundaries. This must not be done in case of vertical nesting
    730733!-- mode as in that case the lateral boundaries are actually cyclic.
    731 !-- @todo: Is this really needed? Boundary values will be overwritten in
    732 !--        coupler or by Inifor data.
    733734    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
    734735       IF ( bc_dirichlet_s )  THEN
  • palm/trunk/SOURCE/check_parameters.f90

    r3337 r3347  
    2525! -----------------
    2626! $Id$
     27! Remove offline nesting in if clause for pressure top boundary condition
     28!
     29! 3343 2018-10-15 10:38:52Z suehring
    2730! (from branch resler)
    2831! Add biometeorology,
     
    19541957       ibc_p_t = 0
    19551958!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
    1956     ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested'  .OR.                 &
    1957              bc_p_t == 'nesting_offline'  )  THEN
     1959    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
    19581960       ibc_p_t = 1
    19591961    ELSE
  • palm/trunk/SOURCE/header.f90

    r3302 r3347  
    2525! -----------------
    2626! $Id$
     27! Header output concerning offline nesting
     28!
     29! 3343 2018-10-15 10:38:52Z suehring
    2730! call of ocean_header
    2831!
     
    438441    USE netcdf_interface,                                                      &
    439442        ONLY:  netcdf_data_format, netcdf_data_format_string, netcdf_deflate
     443       
     444    USE nesting_offl_mod,                                                      &
     445        ONLY:  nesting_offl_header
    440446
    441447    USE ocean_mod,                                                             &
     
    903909!-- Large scale forcing and nudging
    904910    IF ( large_scale_forcing )  CALL lsf_nudging_header( io )
     911   
     912!
     913!-- Offline nesting
     914    IF ( nesting_offline )  CALL nesting_offl_header( io )
    905915
    906916!
  • palm/trunk/SOURCE/init_3d_model.f90

    r3302 r3347  
    2525! -----------------
    2626! $Id$
     27! - Separate offline nesting from large_scale_nudging_mod
     28! - Improve the synthetic turbulence generator
     29!
     30! 3298 2018-10-02 12:21:11Z kanani
     31! Minor formatting (kanani)
     32! Added CALL to the chem_emissions module (Russo)
     33!
     34! 3294 2018-10-01 02:37:10Z raasch
    2735! allocate and set stokes drift velocity profiles
    2836!
     
    563571        ONLY:  chem_emis, chem_emis_att, init_3d,                              &
    564572               netcdf_data_input_init_3d, netcdf_data_input_interpolate
     573       
     574    USE nesting_offl_mod,                                                      &
     575        ONLY:  nesting_offl_init
    565576
    566577    USE ocean_mod,                                                             &
     
    601612
    602613    USE synthetic_turbulence_generator_mod,                                    &
    603         ONLY:  stg_init, use_syn_turb_gen
    604 
     614        ONLY:  parametrize_inflow_turbulence, stg_adjust, stg_init,            &
     615               use_syn_turb_gen
     616               
    605617    USE surface_layer_fluxes_mod,                                              &
    606618        ONLY:  init_surface_layer_fluxes
     
    21522164    IF ( nudging )  CALL nudge_init
    21532165!
    2154 !-- Initialize 1D/3D offline-nesting with COSMO model and read data from
    2155 !-- external file.
    2156     IF ( large_scale_forcing  .OR.  nesting_offline )  CALL lsf_init
     2166!-- Initialize 1D large-scale forcing and nudging and read data from external
     2167!-- ASCII file
     2168    IF ( large_scale_forcing )  CALL lsf_init
    21572169!
    21582170!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
     
    21632175       ENDIF
    21642176    ENDIF
     2177!
     2178!-- Initializae 3D offline nesting in COSMO model and read data from
     2179!-- external NetCDF file.
     2180    IF ( nesting_offline )  CALL nesting_offl_init
    21652181!
    21662182!-- Initialize quantities for special advections schemes
     
    22792295       CALL location_message( 'finished', .TRUE. )
    22802296    ENDIF
    2281 
     2297!
     2298!-- In case the synthetic turbulence generator does not have any information
     2299!-- about the inflow turbulence, these information will be parametrized
     2300!-- depending on the initial atmospheric conditions and surface properties.
     2301!-- Please note, within pre-determined time intervals these turbulence
     2302!-- information can be updated if desired.
     2303    IF ( use_syn_turb_gen  .AND.  parametrize_inflow_turbulence )              &
     2304       CALL stg_adjust
    22822305!
    22832306!-- If required, set chemical emissions
  • palm/trunk/SOURCE/init_pegrid.f90

    r3241 r3347  
    2525! -----------------
    2626! $Id$
     27! Bugfix in setting number of ghost layers in neutral case
     28!
     29! 3341 2018-10-15 10:31:27Z suehring
    2730! unused variables removed
    2831!
     
    723726!
    724727!-- Determine the number of ghost point layers
    725     IF ( ( scalar_advec == 'ws-scheme' .AND. .NOT. neutral ) .OR.             &
     728    IF ( scalar_advec   == 'ws-scheme'  .OR.                                   &
    726729         momentum_advec == 'ws-scheme' )  THEN
    727730       nbgp = 3
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r3274 r3347  
    2525! -----------------
    2626! $Id$
     27! Assign real value instead of integer
     28!
     29! 3341 2018-10-15 10:31:27Z suehring
    2730! Modularization of all bulk cloud physics code components
    2831!
     
    27812784                         vegetation_type_f%var(j,i)                 = 1 ! bare soil
    27822785                         soil_type_f%var_2d(j,i)                    = 1
    2783                          surface_fraction_f%frac(ind_veg_wall,j,i)  = 1
    2784                          surface_fraction_f%frac(ind_pav_green,j,i) = 0
    2785                          surface_fraction_f%frac(ind_wat_win,j,i)   = 0
     2786                         surface_fraction_f%frac(ind_veg_wall,j,i)  = 1.0_wp
     2787                         surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
     2788                         surface_fraction_f%frac(ind_wat_win,j,i)   = 0.0_wp
    27862789                      ENDIF
    27872790                     
  • palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90

    r3294 r3347  
    2525! -----------------
    2626! $Id$
     27! Offline nesting parts are moved to an independent module nesting_offl_mod
     28!
     29! 3341 2018-10-15 10:31:27Z suehring
    2730! ocean renamed ocean_mode
    2831!
     
    8790
    8891    USE control_parameters,                                                    &
    89         ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
    90                bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
     92        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
    9193               constant_heatflux, constant_waterflux,                          &
    9294               data_output_pr, dt_3d, end_time,                                &
     
    9496               ibc_pt_b, ibc_q_b,                                              &
    9597               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
    96                lsf_exception, message_string, nesting_offline, neutral,        &
     98               lsf_exception, message_string, neutral,                         &
    9799               nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
    98100               surface_heatflux, surface_pressure, surface_waterflux,          &
    99101               topography, use_subsidence_tendencies
     102               
     103    USE cpulog,                                                                &
     104        ONLY:  cpu_log, log_point
    100105
    101106    USE grid_variables
     
    107112    USE kinds
    108113
    109     USE netcdf_data_input_mod,                                                 &
    110         ONLY:  nest_offl
    111 
    112114    USE pegrid
    113115
     
    120122    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
    121123    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
    122 
    123     REAL(wp) ::  d_area_t
    124124
    125125    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
     
    156156           lsf_nudging_check_parameters, nudge_init,                           &
    157157           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
    158            lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         &
    159158           nudge, nudge_ref
    160159           
     
    177176 CONTAINS
    178177
    179 
    180 !------------------------------------------------------------------------------!
    181 ! Description:
    182 ! ------------
    183 !> In this subroutine a constant mass within the model domain is guaranteed.
    184 !> Larger-scale models may be based on a compressible equation system, which is
    185 !> not consistent with PALMs incompressible equation system. In order to avoid
    186 !> a decrease or increase of mass during the simulation, non-divergent flow
    187 !> through the lateral and top boundaries is compensated by the vertical wind
    188 !> component at the top boundary.
    189 !------------------------------------------------------------------------------!
    190     SUBROUTINE lsf_nesting_offline_mass_conservation
    191 
    192        USE control_parameters,                                                 &
    193            ONLY:  volume_flow
    194 
    195        IMPLICIT NONE
    196 
    197        INTEGER(iwp) ::  i !< grid index in x-direction
    198        INTEGER(iwp) ::  j !< grid index in y-direction
    199        INTEGER(iwp) ::  k !< grid index in z-direction
    200 
    201        REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
    202        REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
    203 
    204        volume_flow   = 0.0_wp
    205        volume_flow_l = 0.0_wp
    206 
    207        d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
    208 
    209        IF ( bc_dirichlet_l )  THEN
    210           i = nxl
    211           DO  j = nys, nyn
    212              DO  k = nzb+1, nzt
    213                 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
    214                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    215                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    216              ENDDO
    217           ENDDO
    218        ENDIF
    219        IF ( bc_dirichlet_r )  THEN
    220           i = nxr+1
    221           DO  j = nys, nyn
    222              DO  k = nzb+1, nzt
    223                 volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
    224                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    225                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    226              ENDDO
    227           ENDDO
    228        ENDIF
    229        IF ( bc_dirichlet_s )  THEN
    230           j = nys
    231           DO  i = nxl, nxr
    232              DO  k = nzb+1, nzt
    233                 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
    234                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    235                                               BTEST( wall_flags_0(k,j,i), 2 ) )
    236              ENDDO
    237           ENDDO
    238        ENDIF
    239        IF ( bc_dirichlet_n )  THEN
    240           j = nyn+1
    241           DO  i = nxl, nxr
    242              DO  k = nzb+1, nzt
    243                 volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
    244                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    245                                               BTEST( wall_flags_0(k,j,i), 2 ) )
    246              ENDDO
    247           ENDDO
    248        ENDIF
    249 !
    250 !--    Top boundary
    251        k = nzt
    252        DO  i = nxl, nxr
    253           DO  j = nys, nyn
    254              volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
    255           ENDDO
    256        ENDDO
    257 
    258 #if defined( __parallel )
    259        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    260        CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
    261                            comm2d, ierr )
    262 #else
    263        volume_flow = volume_flow_l
    264 #endif
    265 
    266        w_correct = SUM( volume_flow ) * d_area_t
    267 
    268        DO  i = nxl, nxr
    269           DO  j = nys, nyn
    270              DO  k = nzt, nzt + 1
    271                 w(k,j,i) = w(k,j,i) + w_correct
    272              ENDDO
    273           ENDDO
    274        ENDDO
    275 
    276     END SUBROUTINE lsf_nesting_offline_mass_conservation
    277 
    278 
    279 !------------------------------------------------------------------------------!
    280 ! Description:
    281 ! ------------
    282 !> Set the lateral and top boundary conditions in case the PALM domain is
    283 !> nested offline in a mesoscale model.
    284 !------------------------------------------------------------------------------!
    285     SUBROUTINE lsf_nesting_offline
    286 
    287        USE control_parameters,                                                 &
    288            ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
    289                   bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
    290                   rans_tke_e, time_since_reference_point                     
    291 
    292        IMPLICIT NONE
    293 
    294        INTEGER(iwp) ::  i !< running index x-direction
    295        INTEGER(iwp) ::  j !< running index y-direction
    296        INTEGER(iwp) ::  k !< running index z-direction
    297 
    298        REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
    299        REAL(wp) ::  t_ref   !< time past since last reference step
    300      
    301 !
    302 !--    Calculate time interval of forcing data       
    303        ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
    304                             nest_offl%time(nest_offl%tind) )
    305 !
    306 !--    Calculate reziproke time past since last reference step. Please note,
    307 !--    the time coordinate is still not updated, so that the actual time need
    308 !--    to be incremented by dt_3d. Moreover, note that the simulation time
    309 !--    passed since simulation start is time_since_reference_point, not
    310 !--    simulated_time!
    311        t_ref = time_since_reference_point + dt_3d -                            &
    312                                             nest_offl%time(nest_offl%tind)
    313                                            
    314        IF ( bc_dirichlet_l )  THEN
    315 
    316           DO  j = nys, nyn
    317              DO  k = nzb+1, nzt
    318                 u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
    319                        ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
    320                          MERGE( 1.0_wp, 0.0_wp,                                &
    321                                 BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
    322              ENDDO
    323           ENDDO
    324 
    325           DO  j = nys, nyn
    326              DO  k = nzb+1, nzt-1
    327                 w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
    328                        ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
    329                          MERGE( 1.0_wp, 0.0_wp,                                &
    330                                 BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
    331              ENDDO
    332           ENDDO
    333 
    334           DO  j = nysv, nyn
    335              DO  k = nzb+1, nzt
    336                 v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
    337                        ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
    338                          MERGE( 1.0_wp, 0.0_wp,                                &
    339                                 BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
    340              ENDDO
    341           ENDDO
    342 
    343           IF ( .NOT. neutral )  THEN
    344              DO  j = nys, nyn
    345                 DO  k = nzb+1, nzt
    346                    pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
    347                                                                    t_ref   *   &
    348                        ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
    349  
    350                 ENDDO
    351              ENDDO
    352           ENDIF
    353 
    354           IF ( humidity )  THEN
    355              DO  j = nys, nyn
    356                 DO  k = nzb+1, nzt
    357                    q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
    358                                                                  t_ref     *   &
    359                        ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
    360  
    361                 ENDDO
    362              ENDDO
    363           ENDIF
    364 
    365        ENDIF
    366 
    367        IF ( bc_dirichlet_r )  THEN
    368 
    369           DO  j = nys, nyn
    370              DO  k = nzb+1, nzt
    371                 u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
    372                       ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
    373                          MERGE( 1.0_wp, 0.0_wp,                                 &
    374                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
    375 
    376              ENDDO
    377           ENDDO
    378           DO  j = nys, nyn
    379              DO  k = nzb+1, nzt-1
    380                 w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
    381                       ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
    382                          MERGE( 1.0_wp, 0.0_wp,                                 &
    383                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
    384              ENDDO
    385           ENDDO
    386 
    387           DO  j = nysv, nyn
    388              DO  k = nzb+1, nzt
    389                 v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
    390                       ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
    391                          MERGE( 1.0_wp, 0.0_wp,                                 &
    392                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
    393              ENDDO
    394           ENDDO
    395 
    396           IF ( .NOT. neutral )  THEN
    397              DO  j = nys, nyn
    398                 DO  k = nzb+1, nzt
    399                    pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
    400                                                                     t_ref   *  &
    401                      ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
    402  
    403                 ENDDO
    404              ENDDO
    405           ENDIF
    406 
    407           IF ( humidity )  THEN
    408              DO  j = nys, nyn
    409                 DO  k = nzb+1, nzt
    410                    q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
    411                                                                     t_ref   *  &
    412                        ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
    413  
    414                 ENDDO
    415              ENDDO
    416           ENDIF
    417 
    418        ENDIF
    419 
    420        IF ( bc_dirichlet_s )  THEN
    421 
    422           DO  i = nxl, nxr
    423              DO  k = nzb+1, nzt
    424                 v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
    425                       ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
    426                          MERGE( 1.0_wp, 0.0_wp,                                 &
    427                                 BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
    428              ENDDO
    429           ENDDO
    430 
    431           DO  i = nxl, nxr
    432              DO  k = nzb+1, nzt-1
    433                 w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
    434                         ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
    435                            MERGE( 1.0_wp, 0.0_wp,                                &
    436                                   BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
    437              ENDDO
    438           ENDDO
    439 
    440           DO  i = nxlu, nxr
    441              DO  k = nzb+1, nzt
    442                 u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
    443                         ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
    444                            MERGE( 1.0_wp, 0.0_wp,                                &
    445                                   BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
    446              ENDDO
    447           ENDDO
    448 
    449           IF ( .NOT. neutral )  THEN
    450              DO  i = nxl, nxr
    451                 DO  k = nzb+1, nzt
    452                    pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
    453                                                                     t_ref   *  &
    454                      ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
    455  
    456                 ENDDO
    457              ENDDO
    458           ENDIF
    459 
    460           IF ( humidity )  THEN
    461              DO  i = nxl, nxr
    462                 DO  k = nzb+1, nzt
    463                    q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
    464                                                                     t_ref   *  &
    465                        ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
    466  
    467                 ENDDO
    468              ENDDO
    469           ENDIF
    470 
    471        ENDIF
    472 
    473        IF ( bc_dirichlet_n )  THEN
    474 
    475           DO  i = nxl, nxr
    476              DO  k = nzb+1, nzt
    477                 v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
    478                         ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
    479                            MERGE( 1.0_wp, 0.0_wp,                                 &
    480                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
    481              ENDDO
    482           ENDDO
    483           DO  i = nxl, nxr
    484              DO  k = nzb+1, nzt-1
    485                 w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
    486                         ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
    487                            MERGE( 1.0_wp, 0.0_wp,                                &
    488                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
    489              ENDDO
    490           ENDDO
    491 
    492           DO  i = nxlu, nxr
    493              DO  k = nzb+1, nzt
    494                 u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
    495                         ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
    496                            MERGE( 1.0_wp, 0.0_wp,                                &
    497                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
    498 
    499              ENDDO
    500           ENDDO
    501 
    502           IF ( .NOT. neutral )  THEN
    503              DO  i = nxl, nxr
    504                 DO  k = nzb+1, nzt
    505                    pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
    506                                                                     t_ref   *  &
    507                      ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
    508  
    509                 ENDDO
    510              ENDDO
    511           ENDIF
    512 
    513           IF ( humidity )  THEN
    514              DO  i = nxl, nxr
    515                 DO  k = nzb+1, nzt
    516                    q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
    517                                                                     t_ref   *  &
    518                        ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
    519  
    520                 ENDDO
    521              ENDDO
    522           ENDIF
    523 
    524        ENDIF
    525 !
    526 !--    Top boundary.
    527        DO  i = nxlu, nxr
    528           DO  j = nys, nyn
    529              u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
    530                         ( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) *  &
    531                            MERGE( 1.0_wp, 0.0_wp,                              &
    532                                   BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
    533           ENDDO
    534        ENDDO
    535 
    536        DO  i = nxl, nxr
    537           DO  j = nysv, nyn
    538              v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
    539                         ( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) *  &
    540                            MERGE( 1.0_wp, 0.0_wp,                              &
    541                                   BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
    542           ENDDO
    543        ENDDO
    544 
    545        DO  i = nxl, nxr
    546           DO  j = nys, nyn
    547              w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
    548                         ( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) *  &
    549                            MERGE( 1.0_wp, 0.0_wp,                              &
    550                                   BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
    551           ENDDO
    552        ENDDO
    553 
    554 
    555        IF ( .NOT. neutral )  THEN
    556           DO  i = nxl, nxr
    557              DO  j = nys, nyn
    558                 pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
    559                         ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
    560              ENDDO
    561           ENDDO
    562        ENDIF
    563 
    564        IF ( humidity )  THEN
    565           DO  i = nxl, nxr
    566              DO  j = nys, nyn
    567                 q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
    568                         ( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
    569              ENDDO
    570           ENDDO
    571        ENDIF
    572 !
    573 !--    At the edges( left-south, left-north, right-south and right-north) set
    574 !--    data on ghost points.
    575        IF ( bc_dirichlet_l  .AND.  bc_dirichlet_s )  THEN
    576           DO  i = 1, nbgp
    577              u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
    578              w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
    579              IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
    580              IF ( humidity      )  q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
    581           ENDDO
    582           DO  i = 1, nbgp+1
    583              v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
    584           ENDDO
    585        ENDIF
    586        IF ( bc_dirichlet_l  .AND.  bc_dirichlet_n )  THEN
    587           DO  i = 1, nbgp
    588              u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
    589              v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
    590              w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
    591              IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
    592              IF ( humidity      )  q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
    593           ENDDO
    594        ENDIF
    595        IF ( bc_dirichlet_r  .AND.  bc_dirichlet_s )  THEN
    596           DO  i = 1, nbgp
    597              u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
    598              w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
    599              IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
    600              IF ( humidity      )  q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
    601           ENDDO
    602           DO  i = 1, nbgp+1
    603              v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
    604           ENDDO
    605        ENDIF
    606        IF ( bc_dirichlet_r  .AND.  bc_dirichlet_n )  THEN
    607           DO  i = 1, nbgp
    608              u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
    609              v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
    610              w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
    611              IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
    612              IF ( humidity      )  q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
    613           ENDDO
    614        ENDIF
    615 !
    616 !--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
    617 !--    passive scalar, dissipation, and chemical species if required
    618        IF ( rans_mode  .AND.  rans_tke_e )  THEN
    619           IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
    620           IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
    621           IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
    622           IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
    623        ENDIF
    624        IF ( .NOT. constant_diffusion )  THEN
    625           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
    626           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
    627           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
    628           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
    629           e(nzt+1,:,:) = e(nzt,:,:)
    630        ENDIF
    631        IF ( passive_scalar )  THEN
    632           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
    633           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
    634           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
    635           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
    636        ENDIF
    637 
    638 
    639 
    640        CALL exchange_horiz( u, nbgp )
    641        CALL exchange_horiz( v, nbgp )
    642        CALL exchange_horiz( w, nbgp )
    643        IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
    644        IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
    645 
    646 !
    647 !--    Set surface pressure. Please note, time-dependent surface
    648 !--    pressure would require changes in anelastic approximation and
    649 !--    treatment of fluxes.
    650 !--    For the moment, comment this out!
    651 !      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
    652 !                                                      ddt_lsf * t_ref *       &
    653 !                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
    654 !                                    - nest_offl%surface_pressure(nest_offl%tind) )
    655 
    656     END SUBROUTINE lsf_nesting_offline
    657178
    658179!------------------------------------------------------------------------------!
     
    935456           ONLY:  bc_lr_cyc, bc_ns_cyc
    936457
    937        USE netcdf_data_input_mod,                                              &
    938            ONLY:  netcdf_data_input_lsf 
    939 
    940458       IMPLICIT NONE
    941459
     
    967485       REAL(wp) ::  r_dummy           !<
    968486
    969        IF ( nesting_offline )  THEN
    970 !
    971 !--       Allocate arrays for geostrophic wind components. Arrays will
    972 !--       incorporate 2 time levels in order to interpolate in between. Please
    973 !--       note, forcing using geostrophic wind components is only required in
    974 !--       case of cyclic boundary conditions.
    975           IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
    976              ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
    977              ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
     487       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
     488                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
     489                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
     490                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
     491                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
     492                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
     493                 wsubs_vert(nzb:nzt+1,0:nlsf) )
     494
     495       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
     496       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
     497       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
     498       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
     499       wsubs_vert = 0.0_wp
     500
     501!
     502!--    Array for storing large scale forcing and nudging tendencies at each
     503!--    timestep for data output
     504       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
     505       sums_ls_l = 0.0_wp
     506
     507       ngp_sums_ls = (nz+2)*6
     508
     509       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
     510              FORM='FORMATTED', IOSTAT=ierrn )
     511
     512       IF ( ierrn /= 0 )  THEN
     513          message_string = 'file LSF_DATA does not exist'
     514          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
     515       ENDIF
     516
     517       ierrn = 0
     518!
     519!--    First three lines of LSF_DATA contain header
     520       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     521       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     522       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     523
     524       IF ( ierrn /= 0 )  THEN
     525          message_string = 'errors in file LSF_DATA'
     526          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     527       ENDIF
     528
     529!
     530!--    Surface values are read in
     531       nt     = 0
     532       ierrn = 0
     533
     534       DO WHILE ( time_surf(nt) < end_time )
     535          nt = nt + 1
     536          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
     537                                             qsws_surf(nt), pt_surf(nt),       &
     538                                             q_surf(nt), p_surf(nt)
     539
     540          IF ( ierrn /= 0 )  THEN
     541            WRITE ( message_string, * ) 'No time dependent surface ' //        &
     542                              'variables in & LSF_DATA for end of run found'
     543
     544             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
    978545          ENDIF
    979 !
    980 !--       Allocate arrays for reading boundary values. Arrays will incorporate 2
    981 !--       time levels in order to interpolate in between.
    982           IF ( bc_dirichlet_l )  THEN
    983              ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
    984              ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
    985              ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
    986              IF ( humidity )      ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
    987              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
    988           ENDIF
    989           IF ( bc_dirichlet_r )  THEN
    990              ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
    991              ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
    992              ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
    993              IF ( humidity )      ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
    994              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
    995           ENDIF
    996           IF ( bc_dirichlet_n )  THEN
    997              ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
    998              ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
    999              ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
    1000              IF ( humidity )      ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
    1001              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
    1002           ENDIF
    1003           IF ( bc_dirichlet_s )  THEN
    1004              ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
    1005              ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
    1006              ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
    1007              IF ( humidity )      ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
    1008              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
    1009           ENDIF
    1010          
    1011           ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
    1012           ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
    1013           ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
    1014           IF ( humidity )      ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
    1015           IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
    1016 
    1017 !
    1018 !--       Read COSMO data at lateral and top boundaries
    1019           CALL netcdf_data_input_lsf
    1020 !
    1021 !--       Write COSMO data at lateral and top boundaries
    1022           CALL lsf_nesting_offline
    1023 !
    1024 !--       After 3D data is initialized, ensure mass conservation
    1025           CALL lsf_nesting_offline_mass_conservation
    1026 !
    1027 !--       Initialize surface pressure. Please note, time-dependent surface
    1028 !--       pressure would require changes in anelastic approximation and
    1029 !--       treatment of fluxes.
    1030 !--       For the moment, comment this out!
    1031 !         surface_pressure = nest_offl%surface_pressure(0)
    1032 
    1033        ELSE
    1034 
    1035           ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
    1036                     qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
    1037                     td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
    1038                     td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
    1039                     time_vert(0:nlsf), time_surf(0:nlsf),                      &
    1040                     ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
    1041                     wsubs_vert(nzb:nzt+1,0:nlsf) )
    1042 
    1043           p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
    1044           shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
    1045           td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
    1046           time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
    1047           wsubs_vert = 0.0_wp
    1048 
    1049 !
    1050 !--       Array for storing large scale forcing and nudging tendencies at each
    1051 !--       timestep for data output
    1052           ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
    1053           sums_ls_l = 0.0_wp
    1054 
    1055           ngp_sums_ls = (nz+2)*6
    1056 
    1057           OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
    1058                  FORM='FORMATTED', IOSTAT=ierrn )
    1059 
    1060           IF ( ierrn /= 0 )  THEN
    1061              message_string = 'file LSF_DATA does not exist'
    1062              CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
    1063           ENDIF
    1064 
    1065           ierrn = 0
    1066 !
    1067 !--       First three lines of LSF_DATA contain header
    1068           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1069           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1070           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1071 
     546       ENDDO
     547
     548       IF ( time_surf(1) > end_time )  THEN
     549          WRITE ( message_string, * ) 'Time dependent surface variables in ' //&
     550                                      '&LSF_DATA set in after end of ' ,       &
     551                                      'simulation - lsf_surf is set to FALSE'
     552          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
     553          lsf_surf = .FALSE.
     554       ENDIF
     555
     556!
     557!--    Go to the end of the list with surface variables
     558       DO WHILE ( ierrn == 0 )
     559          READ ( finput, *, IOSTAT = ierrn ) r_dummy
     560       ENDDO
     561
     562!
     563!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
     564
     565       nt = 0
     566       DO WHILE ( time_vert(nt) < end_time )
     567          nt = nt + 1
     568          hash = "#"
     569          ierrn = 1 ! not zero
     570!
     571!--       Search for the next line consisting of "# time",
     572!--       from there onwards the profiles will be read
     573          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
     574             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
     575             IF ( ierrn < 0 )  THEN
     576                WRITE( message_string, * ) 'No time dependent vertical profiles',&
     577                                 ' in & LSF_DATA for end of run found'
     578                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
     579             ENDIF
     580          ENDDO
     581
     582          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
     583
     584          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
     585                                           lowwsubs_vert, low_td_lsa_lpt,      &
     586                                           low_td_lsa_q, low_td_sub_lpt,       &
     587                                           low_td_sub_q
    1072588          IF ( ierrn /= 0 )  THEN
    1073589             message_string = 'errors in file LSF_DATA'
     
    1075591          ENDIF
    1076592
    1077 !
    1078 !--       Surface values are read in
    1079           nt     = 0
    1080           ierrn = 0
    1081 
    1082           DO WHILE ( time_surf(nt) < end_time )
    1083              nt = nt + 1
    1084              READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
    1085                                                 qsws_surf(nt), pt_surf(nt),    &
    1086                                                 q_surf(nt), p_surf(nt)
    1087 
    1088              IF ( ierrn /= 0 )  THEN
    1089                WRITE ( message_string, * ) 'No time dependent surface ' //     &
    1090                                  'variables in & LSF_DATA for end of run found'
    1091 
    1092                 CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
     593          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
     594                                           highvg_vert, highwsubs_vert,        &
     595                                           high_td_lsa_lpt, high_td_lsa_q,     &
     596                                           high_td_sub_lpt, high_td_sub_q
     597       
     598          IF ( ierrn /= 0 )  THEN
     599             message_string = 'errors in file LSF_DATA'
     600             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     601          ENDIF
     602
     603
     604          DO  k = nzb, nzt+1
     605             IF ( highheight < zu(k) )  THEN
     606                lowheight      = highheight
     607                lowug_vert     = highug_vert
     608                lowvg_vert     = highvg_vert
     609                lowwsubs_vert  = highwsubs_vert
     610                low_td_lsa_lpt = high_td_lsa_lpt
     611                low_td_lsa_q   = high_td_lsa_q
     612                low_td_sub_lpt = high_td_sub_lpt
     613                low_td_sub_q   = high_td_sub_q
     614
     615                ierrn = 0
     616                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
     617                                                 highvg_vert, highwsubs_vert,  &
     618                                                 high_td_lsa_lpt,              &
     619                                                 high_td_lsa_q,                &
     620                                                 high_td_sub_lpt, high_td_sub_q
     621
     622                IF ( ierrn /= 0 )  THEN
     623                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',     &
     624                        'is higher than the maximum height in LSF_DATA ',      &
     625                        'which is ', lowheight, 'm. Interpolation on PALM ',   &
     626                        'grid is not possible.'
     627                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
     628                ENDIF
     629
    1093630             ENDIF
     631
     632!
     633!--          Interpolation of prescribed profiles in space
     634             fac = (highheight-zu(k))/(highheight - lowheight)
     635
     636             ug_vert(k,nt)    = fac * lowug_vert                               &
     637                                + ( 1.0_wp - fac ) * highug_vert
     638             vg_vert(k,nt)    = fac * lowvg_vert                               &
     639                                + ( 1.0_wp - fac ) * highvg_vert
     640             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
     641                                + ( 1.0_wp - fac ) * highwsubs_vert
     642
     643             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
     644                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
     645             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
     646                                + ( 1.0_wp - fac ) * high_td_lsa_q
     647             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
     648                                + ( 1.0_wp - fac ) * high_td_sub_lpt
     649             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
     650                                + ( 1.0_wp - fac ) * high_td_sub_q
     651
    1094652          ENDDO
    1095653
    1096           IF ( time_surf(1) > end_time )  THEN
    1097              WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
    1098                                          '&LSF_DATA set in after end of ' ,        &
    1099                                          'simulation - lsf_surf is set to FALSE'
    1100              CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
    1101              lsf_surf = .FALSE.
    1102           ENDIF
    1103 
    1104 !
    1105 !--       Go to the end of the list with surface variables
    1106           DO WHILE ( ierrn == 0 )
    1107              READ ( finput, *, IOSTAT = ierrn ) r_dummy
    1108           ENDDO
    1109 
    1110 !
    1111 !--       Profiles of ug, vg and w_subs are read in (large scale forcing)
    1112 
    1113           nt = 0
    1114           DO WHILE ( time_vert(nt) < end_time )
    1115              nt = nt + 1
    1116              hash = "#"
    1117              ierrn = 1 ! not zero
    1118 !
    1119 !--          Search for the next line consisting of "# time",
    1120 !--          from there onwards the profiles will be read
    1121              DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
    1122                 READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
    1123                 IF ( ierrn < 0 )  THEN
    1124                    WRITE( message_string, * ) 'No time dependent vertical profiles',&
    1125                                     ' in & LSF_DATA for end of run found'
    1126                    CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
    1127                 ENDIF
    1128              ENDDO
    1129 
    1130              IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
    1131 
    1132              READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
    1133                                               lowwsubs_vert, low_td_lsa_lpt,    &
    1134                                               low_td_lsa_q, low_td_sub_lpt,     &
    1135                                               low_td_sub_q
    1136              IF ( ierrn /= 0 )  THEN
    1137                 message_string = 'errors in file LSF_DATA'
    1138                 CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    1139              ENDIF
    1140 
    1141              READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
    1142                                               highvg_vert, highwsubs_vert,     &
    1143                                               high_td_lsa_lpt, high_td_lsa_q,  &
    1144                                               high_td_sub_lpt, high_td_sub_q
    1145          
    1146              IF ( ierrn /= 0 )  THEN
    1147                 message_string = 'errors in file LSF_DATA'
    1148                 CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    1149              ENDIF
    1150 
    1151 
    1152              DO  k = nzb, nzt+1
    1153                 IF ( highheight < zu(k) )  THEN
    1154                    lowheight      = highheight
    1155                    lowug_vert     = highug_vert
    1156                    lowvg_vert     = highvg_vert
    1157                    lowwsubs_vert  = highwsubs_vert
    1158                    low_td_lsa_lpt = high_td_lsa_lpt
    1159                    low_td_lsa_q   = high_td_lsa_q
    1160                    low_td_sub_lpt = high_td_sub_lpt
    1161                    low_td_sub_q   = high_td_sub_q
    1162 
    1163                    ierrn = 0
    1164                    READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
    1165                                                     highvg_vert, highwsubs_vert,&
    1166                                                     high_td_lsa_lpt,            &
    1167                                                     high_td_lsa_q,              &
    1168                                                     high_td_sub_lpt, high_td_sub_q
    1169 
    1170                    IF ( ierrn /= 0 )  THEN
    1171                       WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',  &
    1172                            'is higher than the maximum height in LSF_DATA ',   &
    1173                            'which is ', lowheight, 'm. Interpolation on PALM ',&
    1174                            'grid is not possible.'
    1175                       CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
    1176                    ENDIF
    1177 
    1178                 ENDIF
    1179 
    1180 !
    1181 !--             Interpolation of prescribed profiles in space
    1182                 fac = (highheight-zu(k))/(highheight - lowheight)
    1183 
    1184                 ug_vert(k,nt)    = fac * lowug_vert                            &
    1185                                    + ( 1.0_wp - fac ) * highug_vert
    1186                 vg_vert(k,nt)    = fac * lowvg_vert                            &
    1187                                    + ( 1.0_wp - fac ) * highvg_vert
    1188                 wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
    1189                                    + ( 1.0_wp - fac ) * highwsubs_vert
    1190 
    1191                 td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
    1192                                    + ( 1.0_wp - fac ) * high_td_lsa_lpt
    1193                 td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
    1194                                    + ( 1.0_wp - fac ) * high_td_lsa_q
    1195                 td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
    1196                                    + ( 1.0_wp - fac ) * high_td_sub_lpt
    1197                 td_sub_q(k,nt)   = fac * low_td_sub_q                          &
    1198                                    + ( 1.0_wp - fac ) * high_td_sub_q
    1199 
    1200              ENDDO
    1201 
    1202           ENDDO
    1203 
    1204 !
    1205 !--       Large scale vertical velocity has to be zero at the surface
    1206           wsubs_vert(nzb,:) = 0.0_wp
     654       ENDDO
     655
     656!
     657!--    Large scale vertical velocity has to be zero at the surface
     658       wsubs_vert(nzb,:) = 0.0_wp
    1207659   
    1208           IF ( time_vert(1) > end_time )  THEN
    1209              WRITE ( message_string, * ) 'Time dependent large scale profile ',&
    1210                                 'forcing from&LSF_DATA sets in after end of ' ,&
    1211                                 'simulation - lsf_vert is set to FALSE'
    1212              CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
    1213              lsf_vert = .FALSE.
    1214           ENDIF
    1215 
    1216           CLOSE( finput )
    1217 
    1218        ENDIF
     660       IF ( time_vert(1) > end_time )  THEN
     661          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
     662                             'forcing from&LSF_DATA sets in after end of ' ,   &
     663                             'simulation - lsf_vert is set to FALSE'
     664          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
     665          lsf_vert = .FALSE.
     666       ENDIF
     667
     668       CLOSE( finput )
    1219669
    1220670    END SUBROUTINE lsf_init
  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r3337 r3347  
    2525! -----------------
    2626! $Id$
     27! Subroutine renamed
     28!
     29! 3257 2018-09-17 17:11:46Z suehring
    2730! (from branch resler)
    2831! Formatting
     
    678681    END INTERFACE netcdf_data_input_init_lsm
    679682
    680     INTERFACE netcdf_data_input_lsf
    681        MODULE PROCEDURE netcdf_data_input_lsf
    682     END INTERFACE netcdf_data_input_lsf
     683    INTERFACE netcdf_data_input_offline_nesting
     684       MODULE PROCEDURE netcdf_data_input_offline_nesting
     685    END INTERFACE netcdf_data_input_offline_nesting
    683686
    684687    INTERFACE netcdf_data_input_surface_data
     
    735738           netcdf_data_input_init, netcdf_data_input_init_lsm,                 &
    736739           netcdf_data_input_init_3d,                                          &
    737            netcdf_data_input_interpolate, netcdf_data_input_lsf,               &
     740           netcdf_data_input_interpolate, netcdf_data_input_offline_nesting,   &
    738741           netcdf_data_input_surface_data, netcdf_data_input_topo
    739742
     
    30783081!> (COSMO) by Inifor.
    30793082!------------------------------------------------------------------------------!
    3080     SUBROUTINE netcdf_data_input_lsf
     3083    SUBROUTINE netcdf_data_input_offline_nesting
    30813084
    30823085       USE control_parameters,                                                 &
     
    31573160!
    31583161!--    Obtain time index for current input starting at 0.
    3159 !--    @todo: At the moment time, in INIFOR and simulated time correspond
     3162!--    @todo: At the moment INIFOR and simulated time correspond
    31603163!--           to each other. If required, adjust to daytime.
    31613164       nest_offl%tind = MINLOC( ABS( nest_offl%time -                          &
     
    33483351       CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
    33493352
    3350     END SUBROUTINE netcdf_data_input_lsf
     3353    END SUBROUTINE netcdf_data_input_offline_nesting
    33513354
    33523355
  • palm/trunk/SOURCE/parin.f90

    r3313 r3347  
    2525! -----------------
    2626! $Id$
     27! - offline nesting separated from large-scale forcing module
     28! - top boundary condition for pressure in offline nesting changed
     29!
     30! 3343 2018-10-15 10:38:52Z suehring
    2731! Introduced reading of date_init to inipar.(Russo)
    2832! Add extra profiles for chemistry (basit)
     
    477481        ONLY:  mas_parin
    478482
     483    USE nesting_offl_mod,                                                      &
     484        ONLY:  nesting_offl_parin
     485       
    479486    USE netcdf_interface,                                                      &
    480487        ONLY:  netcdf_data_format, netcdf_deflate, netcdf_precision
     
    570577             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
    571578             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    572              most_method, nesting_offline,                  &
     579             most_method,                                                      &
    573580             netcdf_precision, neutral, ngsrb,                                 &
    574581             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
     
    643650             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
    644651             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    645              most_method, nesting_offline,                                     &
     652             most_method,                                                      &
    646653             netcdf_precision, neutral, ngsrb,                                 &
    647654             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
     
    878885          CALL gust_parin
    879886          CALL mas_parin
     887          CALL nesting_offl_parin
    880888          CALL ocean_parin
    881889          CALL pcm_parin
     
    979987           !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear
    980988           !  bc_cs_t  = 'nesting_offline'  ! same for chemical species
    981              bc_p_t   = 'neumann'
     989!
     990!--          For the pressure set Dirichlet conditions, in contrast to the
     991!--          self nesting. This gives less oscilations within the
     992!--          free atmosphere since the pressure solver has more degrees of
     993!--          freedom. In constrast to the self nesting, this might be
     994!--          justified since the top boundary is far away from the domain
     995!--          of interest.
     996             bc_p_t   = 'dirichlet' !'neumann'
    982997          ENDIF
    983998
  • palm/trunk/SOURCE/pres.f90

    r3183 r3347  
    2525! -----------------
    2626! $Id$
     27! Bugfixes in offline nesting.
     28! Add comment.
     29!
     30! 3341 2018-10-15 10:31:27Z suehring
    2731! Rename variables for boundary flags and nesting
    2832!
     
    169173               dt_3d, gathered_size, ibc_p_b, ibc_p_t,                         &
    170174               intermediate_timestep_count, intermediate_timestep_count_max,   &
    171                mg_switch_to_pe0_level, psolver, subdomain_size,                &
     175               mg_switch_to_pe0_level, nesting_offline,                        &
     176               psolver, subdomain_size,                                        &
    172177               topography, volume_flow, volume_flow_area, volume_flow_initial
    173178
     
    386391    ENDIF
    387392
    388     IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.                              &
     393    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline  .AND. &
    389394         .NOT. child_domain_nvn  .AND. intermediate_timestep_count /= 0 )       &
    390395    THEN
     
    719724       volume_flow_l(2) = 0.0_wp
    720725    ENDIF
    721 
     726!
     727!-- Add pressure gradients to the velocity components. Note, the loops are
     728!-- running over the entire model domain, even though, in case of non-cyclic
     729!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
     730!-- respectiveley. However, in case of Dirichlet boundary conditions for the
     731!-- velocities, zero-gradient conditions for the pressure are set, so that
     732!-- no modification is imposed at the boundaries.
    722733    !$OMP PARALLEL PRIVATE (i,j,k)
    723734    !$OMP DO
  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r3289 r3347  
    2525! -----------------
    2626! $Id$
     27! - Revise structure of init routine
     28! - introduce new parameters to skip STG for some timesteps
     29! - introduce time-dependent parametrization of Reynolds-stress tensor
     30! - Bugfix in allocation of mean_inflow_profiles
     31!
     32! 3341 2018-10-15 10:31:27Z suehring
    2733! Introduce module parameter for number of inflow profiles
    2834!
     
    122128! Authors:
    123129! --------
    124 ! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
     130! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
     131!
    125132!
    126133!
     
    133140!> a time scale for each velocity component. The profiles of length and time
    134141!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
    135 !> tensor are read from file STG_PROFILES.
     142!> tensor can be either read from file STG_PROFILES, or will be parametrized
     143!> within the boundary layer.
    136144!>
    137145!> @todo test restart
     
    149157
    150158    USE arrays_3d,                                                             &
    151         ONLY:  mean_inflow_profiles, u, v, w
     159        ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
    152160
    153161    USE basic_constants_and_equations_mod,                                     &
    154         ONLY:  pi
     162        ONLY:  g, kappa, pi
    155163
    156164    USE control_parameters,                                                    &
     
    159167
    160168    USE cpulog,                                                                &
    161         ONLY:  cpu_log, log_point, log_point_s
     169        ONLY:  cpu_log, log_point
    162170
    163171    USE indices,                                                               &
     
    170178#endif
    171179
     180    USE nesting_offl_mod,                                                      &
     181        ONLY:  zi_ribulk
     182
    172183    USE pegrid,                                                                &
    173184        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
     
    184195
    185196
    186     LOGICAL :: velocity_seed_initialized = .FALSE.  !< true after first call of stg_main
    187     LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
     197    LOGICAL ::  velocity_seed_initialized = .FALSE.     !< true after first call of stg_main
     198    LOGICAL ::  parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized
     199    LOGICAL ::  use_syn_turb_gen = .FALSE.              !< switch to use synthetic turbulence generator
    188200
    189201    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
     
    202214    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
    203215
    204     REAL(wp) :: mc_factor    !< mass flux correction factor
    205 
    206     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz      !< displacement for MPI_GATHERV
    207     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz  !< receive count for MPI_GATHERV
    208     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz      !< displacement for MPI_GATHERV
    209     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz  !< receive count for MPI_GATHERV
    210     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux            !< length scale of u in x direction (in gp)
    211     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy            !< length scale of u in y direction (in gp)
    212     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz            !< length scale of u in z direction (in gp)
    213     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx            !< length scale of v in x direction (in gp)
    214     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy            !< length scale of v in y direction (in gp)
    215     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz            !< length scale of v in z direction (in gp)
    216     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx            !< length scale of w in x direction (in gp)
    217     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy            !< length scale of w in y direction (in gp)
    218     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz            !< length scale of w in z direction (in gp)
    219 
    220     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed        !< seed of random number for rn-generator
    221 
    222     REAL(wp), DIMENSION(:), ALLOCATABLE :: a11             !< coefficient for Lund rotation
    223     REAL(wp), DIMENSION(:), ALLOCATABLE :: a21             !< coefficient for Lund rotation
    224     REAL(wp), DIMENSION(:), ALLOCATABLE :: a22             !< coefficient for Lund rotation
    225     REAL(wp), DIMENSION(:), ALLOCATABLE :: a31             !< coefficient for Lund rotation
    226     REAL(wp), DIMENSION(:), ALLOCATABLE :: a32             !< coefficient for Lund rotation
    227     REAL(wp), DIMENSION(:), ALLOCATABLE :: a33             !< coefficient for Lund rotation
    228     REAL(wp), DIMENSION(:), ALLOCATABLE :: tu              !< Lagrangian time scale of u
    229     REAL(wp), DIMENSION(:), ALLOCATABLE :: tv              !< Lagrangian time scale of v
    230     REAL(wp), DIMENSION(:), ALLOCATABLE :: tw              !< Lagrangian time scale of w
    231 
    232     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux           !< filter function for u in x direction
    233     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy           !< filter function for u in y direction
    234     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz           !< filter function for u in z direction
    235     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx           !< filter function for v in x direction
    236     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy           !< filter function for v in y direction
    237     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz           !< filter function for v in z direction
    238     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx           !< filter function for w in y direction
    239     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy           !< filter function for w in y direction
    240     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz           !< filter function for w in z direction
    241     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz         !< velocity seed for u at xz plane
    242     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz        !< velocity seed for u at xz plane with new random number
    243     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz         !< velocity seed for u at yz plane
    244     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz        !< velocity seed for u at yz plane with new random number
    245     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz         !< velocity seed for v at xz plane
    246     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz        !< velocity seed for v at xz plane with new random number
    247     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz         !< velocity seed for v at yz plane
    248     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz        !< velocity seed for v at yz plane with new random number
    249     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz         !< velocity seed for w at xz plane
    250     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz        !< velocity seed for w at xz plane with new random number
    251     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz         !< velocity seed for w at yz plane
    252     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz        !< velocity seed for w at yz plane with new random number
    253 
     216    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_xz      !< displacement for MPI_GATHERV
     217    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_xz  !< receive count for MPI_GATHERV
     218    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_yz      !< displacement for MPI_GATHERV
     219    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_yz  !< receive count for MPI_GATHERV
     220    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nux            !< length scale of u in x direction (in gp)
     221    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuy            !< length scale of u in y direction (in gp)
     222    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuz            !< length scale of u in z direction (in gp)
     223    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvx            !< length scale of v in x direction (in gp)
     224    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvy            !< length scale of v in y direction (in gp)
     225    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvz            !< length scale of v in z direction (in gp)
     226    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwx            !< length scale of w in x direction (in gp)
     227    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwy            !< length scale of w in y direction (in gp)
     228    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
     229    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
     230
     231    REAL(wp) ::  cosmo_ref = 10.0_wp      !< height of first vertical grid level in mesoscale model, used for calculation of scaling parameters
     232    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
     233    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
     234    REAL(wp) ::  mc_factor                !< mass flux correction factor
     235    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
     236    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
     237    REAL(wp) ::  scale_wm                 !< scaling parameter used for turbulence parametrization - momentum scale 
     238    REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
     239    REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
     240   
     241   
     242    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r11              !< Reynolds parameter
     243    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r21              !< Reynolds parameter
     244    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r22              !< Reynolds parameter
     245    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r31              !< Reynolds parameter
     246    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r32              !< Reynolds parameter
     247    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r33              !< Reynolds parameter
     248   
     249    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a11             !< coefficient for Lund rotation
     250    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a21             !< coefficient for Lund rotation
     251    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a22             !< coefficient for Lund rotation
     252    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a31             !< coefficient for Lund rotation
     253    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a32             !< coefficient for Lund rotation
     254    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a33             !< coefficient for Lund rotation
     255    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tu              !< Lagrangian time scale of u
     256    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tv              !< Lagrangian time scale of v
     257    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tw              !< Lagrangian time scale of w
     258
     259    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bux           !< filter function for u in x direction
     260    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buy           !< filter function for u in y direction
     261    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buz           !< filter function for u in z direction
     262    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvx           !< filter function for v in x direction
     263    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvy           !< filter function for v in y direction
     264    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvz           !< filter function for v in z direction
     265    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwx           !< filter function for w in y direction
     266    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwy           !< filter function for w in y direction
     267    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwz           !< filter function for w in z direction
     268    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_xz         !< velocity seed for u at xz plane
     269    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_xz        !< velocity seed for u at xz plane with new random number
     270    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_yz         !< velocity seed for u at yz plane
     271    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_yz        !< velocity seed for u at yz plane with new random number
     272    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_xz         !< velocity seed for v at xz plane
     273    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_xz        !< velocity seed for v at xz plane with new random number
     274    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_yz         !< velocity seed for v at yz plane
     275    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_yz        !< velocity seed for v at yz plane with new random number
     276    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_xz         !< velocity seed for w at xz plane
     277    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_xz        !< velocity seed for w at xz plane with new random number
     278    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_yz         !< velocity seed for w at yz plane
     279    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_yz        !< velocity seed for w at yz plane with new random number
     280   
     281    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz     !< imposed disturbances at north/south boundary
     282    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz     !< imposed disturbances at north/south boundary
    254283
    255284!
    256285!-- PALM interfaces:
     286!-- Adjust time and lenght scales, Reynolds stress, and filter functions
     287    INTERFACE stg_adjust
     288       MODULE PROCEDURE stg_adjust
     289    END INTERFACE stg_adjust
     290!
    257291!-- Input parameter checks to be done in check_parameters
    258292    INTERFACE stg_check_parameters
     
    319353!
    320354!-- Public interfaces
    321     PUBLIC  stg_check_parameters, stg_header, stg_init, stg_main, stg_parin,   &
    322             stg_wrd_global, stg_rrd_global
     355    PUBLIC  stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main,  &
     356            stg_parin, stg_rrd_global, stg_wrd_global
    323357
    324358!
    325359!-- Public variables
    326     PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
    327             use_syn_turb_gen
     360    PUBLIC  dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north,             &
     361            id_stg_right, id_stg_south, parametrize_inflow_turbulence,         &
     362            time_stg_adjust, time_stg_call, use_syn_turb_gen
    328363
    329364
     
    376411             message_string = 'Using synthetic turbulence generator ' //       &
    377412                              'requires %initializing_actions = '         //   &
    378                               '"set_constant_profiles" or "read_restart_data"'
     413                              '"set_constant_profiles" or "read_restart_data"' //&
     414                              ', if not offline nesting is applied.'
    379415             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
    380416          ENDIF
     
    382418          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
    383419             message_string = 'Using synthetic turbulence generator ' //       &
    384                               'requires &bc_lr = "dirichlet/radiation"'
     420                              'requires &bc_lr = "dirichlet/radiation", ' //   &
     421                              'if not offline nesting is applied.'
    385422             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
    386423          ENDIF
    387424          IF ( bc_ns /= 'cyclic' )  THEN
    388425             message_string = 'Using synthetic turbulence generator ' //       &
    389                               'requires &bc_ns = "cyclic"'
     426                              'requires &bc_ns = "cyclic", ' //               &
     427                              'if not offline nesting is applied.'
    390428             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
    391429          ENDIF
     
    425463       WRITE( io, 3 )
    426464    ENDIF
     465   
     466    IF ( parametrize_inflow_turbulence )  THEN
     467       WRITE( io, 4 ) dt_stg_adjust
     468    ELSE
     469       WRITE( io, 5 )
     470    ENDIF
    427471
    4284721   FORMAT (//' Synthetic turbulence generator information:'/                  &
    429473              ' ------------------------------------------'/)
    430 2   FORMAT ('    synthetic turbulence generator switched on')
    431 3   FORMAT ('    synthetic turbulence generator switched off')
     4742   FORMAT ('    synthetic turbulence generator is switched on')
     4753   FORMAT ('    synthetic turbulence generator is switched off')
     4764   FORMAT ('    imposed turbulence statistics are parametrized and ' \\       &
     477            'ajdusted to boundary-layer development each ', F8.2, ' s' )
     4785   FORMAT ('    imposed turbulence is read from file' )
    432479
    433480 END SUBROUTINE stg_header
     
    443490
    444491    USE arrays_3d,                                                             &
    445         ONLY:  dzw, ddzw, u_init, v_init, zu, zw
     492        ONLY:  dzw, ddzw, u_init, v_init, zu
    446493
    447494    USE control_parameters,                                                    &
     
    460507    IMPLICIT NONE
    461508
    462     LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
     509    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
    463510
    464511#if defined( __parallel )
     
    484531    REAL(wp) :: d21     !< vertical interpolation for a21
    485532    REAL(wp) :: d22     !< vertical interpolation for a22
    486     REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
    487533    REAL(wp) :: luy     !< length scale for u in y direction
    488534    REAL(wp) :: luz     !< length scale for u in z direction
     
    494540    REAL(wp) :: zz      !< height
    495541
    496     REAL(wp) :: length_scale_surface !< typical length scale
    497     REAL(wp) :: r_ii_0               !< correction factor
    498     REAL(wp) :: time_scale           !< typical time scale
    499     REAL(wp) :: length_scale_z       !< typical length scale
    500 
    501     REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
    502     REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
    503     REAL(wp),DIMENSION(nzb:nzt+1) :: r22  !< Reynolds parameter
    504     REAL(wp),DIMENSION(nzb:nzt+1) :: r31  !< Reynolds parameter
    505     REAL(wp),DIMENSION(nzb:nzt+1) :: r32  !< Reynolds parameter
    506     REAL(wp),DIMENSION(nzb:nzt+1) :: r33  !< Reynolds parameter
    507 
    508542
    509543#if defined( __parallel )
     
    512546
    513547    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
    514 
    515     IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                             &
    516        ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
    517 
    518     ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
    519                a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
    520                nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
    521                nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
    522                nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
    523                tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
    524548
    525549#if defined( __parallel )
     
    615639    ENDDO
    616640    CALL RANDOM_SEED(put=seed)
    617 
     641!
     642!-- Allocate required arrays
     643!-- mean_inflow profiles must not be allocated in offline nesting
     644    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
     645       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                          &
     646          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
     647    ENDIF
     648
     649    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
     650               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
     651               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
     652               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
     653               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
     654               r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1),                 &
     655               r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1),                 &
     656               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
     657               
     658    ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
     659    ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
     660    dist_xz = 0.0_wp
     661    dist_yz = 0.0_wp
     662!
    618663!-- Read inflow profile
    619664!-- Height levels of profiles in input profile is as follows:
     
    669714       
    670715       CLOSE( 90 )
    671 
     716!
     717!--    Calculate coefficient matrix from Reynolds stress tensor 
     718!--    (Lund rotation)
     719       CALL calc_coeff_matrix
     720!
     721!-- No information about turbulence and its length scales are available.
     722!-- Instead, parametrize turbulence which is imposed at the boundaries.
     723!-- Set flag which indicates that turbulence is parametrized, which is done
     724!-- later when energy-balance models are already initialized. This is
     725!-- because the STG needs information about surface properties such as
     726!-- roughness to build 'realistic' turbulence profiles.
    672727    ELSE
    673728!
    674 !--    Set-up default turbulent length scales. From the numerical point of
    675 !--    view the imposed perturbations should not be immediately dissipated
    676 !--    by the numerics. The numerical dissipation, however, acts on scales
    677 !--    up to 8 x the grid spacing. For this reason, set the turbulence
    678 !--    length scale to 8 time the grid spacing. Further, above the boundary
    679 !--    layer height, set turbulence lenght scales to zero (equivalent to not
    680 !--    imposing any perturbations) in order to save computational costs.
    681 !--    Typical length (time) scales of 100 m (s) should be a good compromise
    682 !--    between all stratrifications. Near-surface variances are fixed to
    683 !--    0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
    684 !--    Vertical fluxes
    685        length_scale_surface = 100.0_wp
    686        r_ii_0               = 0.1_wp
    687        time_scale           = 100.0_wp
    688        DO  k = nzb+1, nzt+1
    689           dum_exp        = MERGE( -zu(k) / ( 0.3* zu(nzt) ),                   &
    690                                   0.0_wp,                                      &
    691                                   zu(k) > 0.3 * zu(nzt)                        &
    692                                 )
    693           length_scale_z = 8.0_wp * MIN( dx, dy, dzw(k) )
    694          
    695 !           IF ( zu(k) > zi_richardson ) 
    696 
    697           nux(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    698           nuy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    699           nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    700           nvx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    701           nvy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    702           nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    703           nwx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    704           nwy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    705           nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    706 
    707           r11(k) = r_ii_0 * EXP( dum_exp )
    708           r22(k) = r_ii_0 * EXP( dum_exp )
    709           r33(k) = r_ii_0 * EXP( dum_exp )
    710 
    711           r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    712           r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    713           r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    714 
    715           tu(k)  = time_scale
    716           tv(k)  = time_scale
    717           tw(k)  = time_scale
    718 
    719        ENDDO
    720        
    721        nux(nzb) = nux(nzb+1)
    722        nuy(nzb) = nuy(nzb+1)
    723        nuz(nzb) = nuz(nzb+1)
    724        nvx(nzb) = nvx(nzb+1)
    725        nvy(nzb) = nvy(nzb+1)
    726        nvz(nzb) = nvz(nzb+1)
    727        nwx(nzb) = nwx(nzb+1)
    728        nwy(nzb) = nwy(nzb+1)
    729        nwz(nzb) = nwz(nzb+1)
    730 
    731        r11(nzb) = r11(nzb+1)
    732        r22(nzb) = r22(nzb+1)
    733        r33(nzb) = r33(nzb+1)
    734 
    735        r21(nzb) = r11(nzb+1)
    736        r31(nzb) = r31(nzb+1)
    737        r32(nzb) = r32(nzb+1)
    738 
    739        tu(nzb)  = time_scale
    740        tv(nzb)  = time_scale
    741        tw(nzb)  = time_scale
    742 
     729!--    Set flag indicating that turbulence is parametrized
     730       parametrize_inflow_turbulence = .TRUE.
     731!
     732!--    Determine boundary-layer depth, which is used to initialize lenght
     733!--    scales
     734       CALL calc_scaling_variables
     735!
     736!--    Initialize lenght and time scales, which in turn are used
     737!--    to initialize the filter functions.
     738       CALL calc_length_and_time_scale
     739           
    743740    ENDIF
    744741
    745742!
    746 !-- Assign initial profiles
     743!-- Assign initial profiles. Note, this is only required if turbulent
     744!-- inflow from the left is desired, not in case of any of the
     745!-- nesting (offline or self nesting) approaches.
    747746    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
    748747       u_init = mean_inflow_profiles(:,1)
     
    751750       e_init = MAXVAL( mean_inflow_profiles(:,5) )
    752751    ENDIF
    753 !
    754 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
    755     DO  k = nzb, nzt+1
    756        IF ( r11(k) > 0.0_wp )  THEN
    757           a11(k) = SQRT( r11(k) )
    758           a21(k) = r21(k) / a11(k)
    759        ELSE
    760           a11(k) = 0.0_wp
    761           a21(k) = 0.0_wp
    762        ENDIF
    763 
    764        a22(k) = r22(k) - a21(k)**2
    765        IF ( a22(k) > 0.0_wp )  THEN
    766           a22(k) = SQRT( a22(k) )
    767        ELSE
    768           a22(k) = 0.0_wp
    769        ENDIF
    770 
    771 !
    772 !--    a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
    773 !--    d21, d22) because of different vertical grid
    774        IF ( k .le. nzt )  THEN
    775           d11 = 0.5_wp * ( r11(k) + r11(k+1) )
    776           IF ( d11 > 0.0_wp )  THEN
    777              d11 = SQRT( d11 )
    778              d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
    779              a31(k) = r31(k) / d11
    780           ELSE
    781              d21 = 0.0_wp
    782              a31(k) = 0.0_wp
    783           ENDIF
    784 
    785           d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
    786           IF ( d22 > 0.0_wp )  THEN
    787              a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
    788           ELSE
    789              a32(k) = 0.0_wp
    790           ENDIF
    791 
    792           a33(k) = r33(k) - a31(k) ** 2 - a32(k) ** 2
    793           IF ( a33(k) > 0.0_wp )  THEN
    794              a33(k) = SQRT( a33(k) )
    795           ELSE
    796              a33(k) = 0.0_wp
    797           ENDIF
    798        ELSE
    799           a31(k) = a31(k-1)
    800           a32(k) = a32(k-1)
    801           a33(k) = a33(k-1)
    802        ENDIF
    803 
    804     ENDDO
     752   
    805753!
    806754!-- Define the size of the filter functions and allocate them.
     
    961909    ENDDO
    962910
    963 END SUBROUTINE stg_filter_func
     911 END SUBROUTINE stg_filter_func
    964912
    965913
     
    977925
    978926
    979     NAMELIST /stg_par/  use_syn_turb_gen
     927    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
    980928
    981929    line = ' '
     
    11081056    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
    11091057
    1110     REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
    1111     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
    1112 
    1113 
    11141058    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
    11151059
    11161060!
    11171061!-- Calculate time step which is needed for filter functions
    1118     dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
    1119 
     1062    dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count)
    11201063!
    11211064!-- Initial value of fu, fv, fw
     
    11451088       velocity_seed_initialized = .TRUE.
    11461089    ENDIF
     1090   
    11471091!
    11481092!-- New set of fu, fv, fw
     
    11501094    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
    11511095    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
    1152 
     1096   
    11531097    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
    1154                              .AND.  .NOT.  rans_mode ) )  THEN
    1155 !
    1156 !--       Generate turbulence at right boundary
    1157           CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
    1158           CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
    1159           CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
    1160 !
    1161 !--       Generate turbulence at north boundary
    1162           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
    1163           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
    1164           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
    1165 !
    1166 !--       Generate turbulence at south boundary
    1167           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
    1168           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
    1169           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1098                          .AND.  .NOT.  rans_mode ) )  THEN
     1099!
     1100!--    Generate turbulence at right boundary
     1101       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
     1102       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
     1103       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
     1104!
     1105!--    Generate turbulence at north boundary
     1106       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
     1107       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
     1108       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
     1109!
     1110!--    Generate turbulence at south boundary
     1111       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
     1112       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
     1113       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
    11701114    ENDIF
     1115   
    11711116!
    11721117!-- Turbulence generation at left and or right boundary
    11731118    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
    1174 
     1119   
    11751120       DO  j = nysg, nyng
    11761121          DO  k = nzb, nzt + 1
     
    11801125                fu_yz(k,j) = fuo_yz(k,j)
    11811126             ELSE
    1182                 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
     1127                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
    11831128                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
    11841129             ENDIF
     
    11871132                fv_yz(k,j) = fvo_yz(k,j)
    11881133             ELSE
    1189                 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
     1134                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
    11901135                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    11911136             ENDIF
     
    11941139                fw_yz(k,j) = fwo_yz(k,j)
    11951140             ELSE
    1196                 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
     1141                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
    11971142                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    11981143             ENDIF
     
    12041149                dist_yz(k,j,2) = 0.0_wp
    12051150                dist_yz(k,j,3) = 0.0_wp
    1206 !                 dist_yz(k,j,4) = 0.0_wp
    1207 !                 dist_yz(k,j,5) = 0.0_wp
    12081151             ELSE
    1209                 dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
     1152                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
    12101153                !experimental test of 1.2
    1211                 dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
    1212                                          * 1.2_wp )                            &
    1213                                        * (   a21(k) * fu_yz(k,j)               &
    1214                                            + a22(k) * fv_yz(k,j) )
    1215                 dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
    1216                                          * 1.3_wp )                            &
    1217                                        * (   a31(k) * fu_yz(k,j)               &
    1218                                            + a32(k) * fv_yz(k,j)               &
    1219                                            + a33(k) * fw_yz(k,j) )
    1220                 ! Calculation for pt and e not yet implemented
    1221 !                 dist_yz(k,j,4) = 0.0_wp
    1222 !                 dist_yz(k,j,5) = 0.0_wp
     1154                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
     1155                                      * 1.2_wp )                               &
     1156                                     * (   a21(k) * fu_yz(k,j)                 &
     1157                                         + a22(k) * fv_yz(k,j) ), 3.0_wp )
     1158                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
     1159                                      * 1.3_wp )                               &
     1160                                    * (   a31(k) * fu_yz(k,j)                  &
     1161                                        + a32(k) * fv_yz(k,j)                  &
     1162                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )
    12231163             ENDIF
    12241164
     
    12411181
    12421182#if defined( __parallel )
    1243           CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
    1244                               1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1183          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
     1184                              comm1dy, ierr )
    12451185#else
    12461186          mc_factor = mc_factor_l
     
    12981238!--       Add disturbances
    12991239          IF ( myidx == id_stg_left  )  THEN
    1300 
    13011240             DO  j = nys, nyn
    1302                 DO  k = nzb, nzt
     1241                DO  k = nzb+1, nzt
    13031242                   u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
    13041243                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
    13051244                                            BTEST( wall_flags_0(k,j,0), 1 ) )
     1245                   u(k,j,-1) = u(k,j,0)
    13061246                   v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
    13071247                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
     
    13141254          ENDIF
    13151255          IF ( myidx == id_stg_right  )  THEN
    1316 
    13171256             DO  j = nys, nyn
    1318                 DO  k = nzb, nzt
     1257                DO  k = nzb+1, nzt
    13191258                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
    13201259                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
     
    13351274!-- Turbulence generation at north and south boundary
    13361275    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
    1337 
     1276   
    13381277       DO  i = nxlg, nxrg
    13391278          DO  k = nzb, nzt + 1
     
    13511290             ELSE
    13521291                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
    1353                          fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1292                      fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    13541293             ENDIF
    1355 
     1294           
    13561295             IF ( tw(k) == 0.0_wp )  THEN
    13571296                fw_xz(k,i) = fwo_xz(k,i)
    13581297             ELSE
    13591298                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
    1360                          fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1299                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    13611300             ENDIF
    13621301!
     
    13691308
    13701309             ELSE
    1371                 dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
     1310                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
    13721311                !experimental test of 1.2
    1373                 dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
    1374                                          * 1.2_wp )                            &
    1375                                        * (   a21(k) * fu_xz(k,i)               &
    1376                                            + a22(k) * fv_xz(k,i) )
    1377                 dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
    1378                                          * 1.3_wp )                            &
    1379                                        * (   a31(k) * fu_xz(k,i)               &
    1380                                            + a32(k) * fv_xz(k,i)               &
    1381                                            + a33(k) * fw_xz(k,i) )
     1312                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
     1313                                      * 1.2_wp )                               &
     1314                                     * (   a21(k) * fu_xz(k,i)                 &
     1315                                         + a22(k) * fv_xz(k,i) ), 3.0_wp )
     1316                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
     1317                                      * 1.3_wp )                               &
     1318                                    * (   a31(k) * fu_xz(k,i)                  &
     1319                                        + a32(k) * fv_xz(k,i)                  &
     1320                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )
    13821321             ENDIF
    13831322
    13841323          ENDDO
    13851324       ENDDO
     1325
    13861326!
    13871327!--    Mass flux correction following Kim et al. (2013)
     
    14251365
    14261366          DO  i = nxl, nxr
    1427              DO  k = nzb, nzt
     1367             DO  k = nzb+1, nzt
    14281368                u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
    14291369                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
     
    14321372                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
    14331373                                           BTEST( wall_flags_0(k,0,i), 2 ) )
     1374                v(k,-1,i) = v(k,0,i)
    14341375                w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
    14351376                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
     
    14411382
    14421383          DO  i = nxl, nxr
    1443              DO  k = nzb, nzt
    1444                 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1445                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1384             DO  k = nzb+1, nzt
     1385                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
     1386                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14461387                                          BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
    1447                 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1448                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1388                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
     1389                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14491390                                          BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
    1450                 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1451                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1391                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
     1392                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14521393                                          BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
    14531394             ENDDO
     
    14551396       ENDIF
    14561397    ENDIF
     1398!
     1399!-- Finally, set time counter for calling STG to zero
     1400    time_stg_call = 0.0_wp
    14571401
    14581402    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
     
    14931437    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    14941438
    1495     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
    1496     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
     1439    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
     1440    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    14971441    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
    14981442    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
     
    16131557    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
    16141558    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
    1615     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
     1559    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_x*2
    16161560    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
    16171561
     
    16201564    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    16211565
    1622     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
    1623     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
     1566    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
     1567    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    16241568    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
    16251569    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
     
    17121656 END SUBROUTINE stg_generate_seed_xz
    17131657
     1658!------------------------------------------------------------------------------!
     1659! Description:
     1660! ------------
     1661!> Parametrization of the Reynolds stress tensor, following the parametrization
     1662!> described in Rotach et al. (1996), which is applied in state-of-the-art
     1663!> dispserion modelling. Please note, the parametrization does not distinguish
     1664!> between along-wind and cross-wind turbulence.
     1665!------------------------------------------------------------------------------!
     1666 SUBROUTINE parametrize_reynolds_stress
     1667
     1668    USE arrays_3d,                                                             &
     1669        ONLY:  zu
     1670
     1671    IMPLICIT NONE
     1672
     1673    INTEGER(iwp) :: k            !< loop index in z-direction
     1674   
     1675    REAL(wp)     ::  zzi         !< ratio of z/zi
     1676   
     1677!
     1678!--
     1679    DO  k = nzb+1, nzt+1
     1680         
     1681       IF ( zu(k) <= zi_ribulk )  THEN
     1682!
     1683!--       Determine normalized height coordinate
     1684          zzi = zu(k) / zi_ribulk
     1685!
     1686!--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
     1687!--       to the denominator, else the merge-function can crash if scale_l is
     1688!--       zero.
     1689          r11(k) = scale_us**2 * (                                             &
     1690                      MERGE( 0.35_wp * (                                       &
     1691                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
     1692                                       )**( 2.0_wp / 3.0_wp ),                 &
     1693                             0.0_wp,                                           &
     1694                             scale_l < 0.0_wp )                                &
     1695                    + 5.0_wp - 4.0_wp * zzi                                    &
     1696                                 )
     1697                                 
     1698          r22(k) = r11(k)
     1699!
     1700!--       w'w'
     1701          r33(k) = scale_wm**2 * (                                             &
     1702                      1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
     1703                    + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
     1704                                 )
     1705!
     1706!--       u'w' and v'w'. Assume isotropy.
     1707          r31(k) = - scale_us**2 * (                                           &
     1708                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
     1709                                   )
     1710                             
     1711          r32(k) = r31(k)
     1712!
     1713!--       For u'v' no parametrization exist so far - ?. For simplicity assume
     1714!--       a similar profile as for u'w'.
     1715          r21(k) = r31(k)
     1716!
     1717!--    Above the boundary layer, assmume laminar flow conditions.
     1718       ELSE
     1719          r11(k) = 10E-8_wp
     1720          r22(k) = 10E-8_wp
     1721          r33(k) = 10E-8_wp
     1722          r21(k) = 10E-8_wp
     1723          r31(k) = 10E-8_wp
     1724          r32(k) = 10E-8_wp
     1725       ENDIF
     1726!        write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l
     1727    ENDDO
     1728
     1729!
     1730!-- Set bottom boundary condition   
     1731    r11(nzb) = r11(nzb+1)
     1732    r22(nzb) = r22(nzb+1)
     1733    r33(nzb) = r33(nzb+1)
     1734
     1735    r21(nzb) = r11(nzb+1)
     1736    r31(nzb) = r31(nzb+1)
     1737    r32(nzb) = r32(nzb+1)
     1738   
     1739
     1740 END SUBROUTINE parametrize_reynolds_stress
     1741 
     1742!------------------------------------------------------------------------------!
     1743! Description:
     1744! ------------
     1745!> Calculate the coefficient matrix from the Lund rotation.
     1746!------------------------------------------------------------------------------!
     1747 SUBROUTINE calc_coeff_matrix
     1748
     1749    IMPLICIT NONE
     1750
     1751    INTEGER(iwp) :: k   !< loop index in z-direction
     1752   
     1753!
     1754!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
     1755    DO  k = nzb+1, nzt+1
     1756       IF ( r11(k) > 0.0_wp )  THEN
     1757          a11(k) = SQRT( r11(k) )
     1758          a21(k) = r21(k) / a11(k)
     1759          a31(k) = r31(k) / a11(k)
     1760       ELSE
     1761          a11(k) = 10E-8_wp
     1762          a21(k) = 10E-8_wp
     1763          a31(k) = 10E-8_wp
     1764       ENDIF
     1765    ENDDO
     1766    DO  k = nzb+1, nzt+1
     1767       a22(k) = r22(k) - a21(k)**2
     1768       IF ( a22(k) > 0.0_wp )  THEN
     1769          a22(k) = SQRT( a22(k) )
     1770          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
     1771       ELSE
     1772          a22(k) = 10E-8_wp
     1773          a32(k) = 10E-8_wp
     1774       ENDIF
     1775    ENDDO
     1776    DO  k = nzb+1, nzt+1
     1777       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
     1778       IF ( a33(k) > 0.0_wp )  THEN
     1779          a33(k) =  SQRT( a33(k) )
     1780       ELSE
     1781          a33(k) = 10E-8_wp
     1782       ENDIF
     1783    ENDDO   
     1784!
     1785!-- Set bottom boundary condition
     1786    a11(nzb) = a11(nzb+1)
     1787    a22(nzb) = a22(nzb+1)
     1788    a21(nzb) = a21(nzb+1)
     1789    a33(nzb) = a33(nzb+1)   
     1790    a31(nzb) = a31(nzb+1)
     1791    a32(nzb) = a32(nzb+1)   
     1792
     1793 END SUBROUTINE calc_coeff_matrix
     1794 
     1795!------------------------------------------------------------------------------!
     1796! Description:
     1797! ------------
     1798!> This routine controls the re-adjustment of the turbulence statistics used
     1799!> for generating turbulence at the lateral boundaries. 
     1800!------------------------------------------------------------------------------!
     1801 SUBROUTINE stg_adjust
     1802
     1803    IMPLICIT NONE
     1804   
     1805    INTEGER(iwp) ::  k
     1806
     1807    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
     1808!
     1809!-- Compute mean boundary layer height according to Richardson-Bulk
     1810!-- criterion using the inflow profiles. Further velocity scale as well as
     1811!-- mean friction velocity are calculated.
     1812    CALL calc_scaling_variables
     1813!
     1814!-- Set length and time scales depending on boundary-layer height
     1815    CALL calc_length_and_time_scale
     1816!
     1817!-- Parametrize Reynolds-stress tensor, diagonal elements as well
     1818!-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows
     1819!-- Rotach et al. (1996) and is based on boundary-layer depth,
     1820!-- friction velocity and velocity scale.
     1821    CALL parametrize_reynolds_stress
     1822!
     1823!-- Calculate coefficient matrix from Reynolds stress tensor 
     1824!-- (Lund rotation)
     1825    CALL calc_coeff_matrix
     1826!
     1827!-- Determine filter functions on basis of updated length scales
     1828    CALL stg_filter_func( nux, bux ) !filter ux
     1829    CALL stg_filter_func( nuy, buy ) !filter uy
     1830    CALL stg_filter_func( nuz, buz ) !filter uz
     1831    CALL stg_filter_func( nvx, bvx ) !filter vx
     1832    CALL stg_filter_func( nvy, bvy ) !filter vy
     1833    CALL stg_filter_func( nvz, bvz ) !filter vz
     1834    CALL stg_filter_func( nwx, bwx ) !filter wx
     1835    CALL stg_filter_func( nwy, bwy ) !filter wy
     1836    CALL stg_filter_func( nwz, bwz ) !filter wz
     1837!
     1838!-- Reset time counter for controlling next adjustment to zero
     1839    time_stg_adjust = 0.0_wp
     1840   
     1841    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
     1842   
     1843 END SUBROUTINE stg_adjust
     1844 
     1845 
     1846!------------------------------------------------------------------------------!
     1847! Description:
     1848! ------------
     1849!> Calculates turbuluent length and time scales if these are not available
     1850!> from measurements.
     1851!------------------------------------------------------------------------------!
     1852 SUBROUTINE calc_length_and_time_scale
     1853
     1854    USE arrays_3d,                                                             &
     1855        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
     1856
     1857    USE grid_variables,                                                        &
     1858        ONLY:  ddx, ddy, dx, dy
     1859       
     1860    USE statistics,                                                            &
     1861        ONLY:  hom
     1862
     1863    IMPLICIT NONE
     1864
     1865
     1866    INTEGER(iwp) :: k            !< loop index in z-direction
     1867
     1868    REAL(wp) ::  length_scale    !< typical length scale
     1869   
     1870!
     1871!-- In initial call the boundary-layer depth can be zero. This case, set
     1872!-- minimum value for boundary-layer depth, to setup length scales correctly.
     1873    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
     1874!
     1875!-- Set-up default turbulent length scales. From the numerical point of
     1876!-- view the imposed perturbations should not be immediately dissipated
     1877!-- by the numerics. The numerical dissipation, however, acts on scales
     1878!-- up to 8 x the grid spacing. For this reason, set the turbulence
     1879!-- length scale to 8 time the grid spacing. Further, above the boundary
     1880!-- layer height, set turbulence lenght scales to zero (equivalent to not
     1881!-- imposing any perturbations) in order to save computational costs.
     1882!-- Typical time scales are derived by assuming Taylors's hypothesis,
     1883!-- using the length scales and the mean profiles of u- and v-component.
     1884    DO  k = nzb+1, nzt+1
     1885   
     1886       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
     1887         
     1888       IF ( zu(k) <= zi_ribulk )  THEN 
     1889!
     1890!--       Assume isotropic turbulence length scales
     1891          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
     1892          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1893          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1894          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
     1895          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1896          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1897          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
     1898          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1899          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1900!
     1901!--       Limit time scales, else they become very larger for low wind speed,
     1902!--       imposing long-living inflow perturbations which in turn propagate
     1903!--       further into the model domain. Use u_init and v_init to calculate
     1904!--       the time scales, which will be equal to the inflow profiles, both,
     1905!--       in offline nesting mode or in dirichlet/radiation mode.
     1906          tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
     1907                                  ( ABS( u_init(k) ) + 0.1_wp ) )
     1908          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
     1909                                  ( ABS( v_init(k) ) + 0.1_wp ) )
     1910!
     1911!--       Time scale of w-component is a mixture from u- and v-component.
     1912          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
     1913!
     1914!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
     1915!--    is not correlated in space and time, just white noise. This saves
     1916!--    computations power.
     1917       ELSE
     1918          nux(k) = 0.0_wp
     1919          nuy(k) = 0.0_wp
     1920          nuz(k) = 0.0_wp
     1921          nvx(k) = 0.0_wp
     1922          nvy(k) = 0.0_wp
     1923          nvz(k) = 0.0_wp
     1924          nwx(k) = 0.0_wp
     1925          nwy(k) = 0.0_wp
     1926          nwz(k) = 0.0_wp
     1927         
     1928          tu(k)  = 0.0_wp
     1929          tv(k)  = 0.0_wp
     1930          tw(k)  = 0.0_wp
     1931       ENDIF
     1932    ENDDO
     1933!
     1934!-- Set bottom boundary condition for the length and time scales
     1935    nux(nzb) = nux(nzb+1)
     1936    nuy(nzb) = nuy(nzb+1)
     1937    nuz(nzb) = nuz(nzb+1)
     1938    nvx(nzb) = nvx(nzb+1)
     1939    nvy(nzb) = nvy(nzb+1)
     1940    nvz(nzb) = nvz(nzb+1)
     1941    nwx(nzb) = nwx(nzb+1)
     1942    nwy(nzb) = nwy(nzb+1)
     1943    nwz(nzb) = nwz(nzb+1)
     1944
     1945    tu(nzb)  = tu(nzb+1)
     1946    tv(nzb)  = tv(nzb+1)
     1947    tw(nzb)  = tw(nzb+1)
     1948
     1949
     1950 END SUBROUTINE calc_length_and_time_scale
     1951
     1952!------------------------------------------------------------------------------!
     1953! Description:
     1954! ------------
     1955!> Calculate scaling variables which are used for turbulence parametrization
     1956!> according to Rotach et al. (1996). Scaling variables are: friction velocity,
     1957!> boundary-layer depth, momentum velocity scale, and Obukhov length.
     1958!------------------------------------------------------------------------------!
     1959 SUBROUTINE calc_scaling_variables
     1960
     1961
     1962    USE control_parameters,                                                    &
     1963        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
     1964               humidity, pt_surface
     1965             
     1966    USE indices,                                                               &
     1967        ONLY:  nx, ny 
     1968       
     1969    USE surface_mod,                                                           &
     1970        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
     1971
     1972    IMPLICIT NONE
     1973
     1974    INTEGER(iwp) :: i            !< loop index in x-direction
     1975    INTEGER(iwp) :: j            !< loop index in y-direction
     1976    INTEGER(iwp) :: k            !< loop index in z-direction
     1977    INTEGER(iwp) :: k_ref        !< index in z-direction for reference height
     1978    INTEGER(iwp) :: m            !< surface element index
     1979
     1980    REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
     1981    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
     1982    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
     1983    REAL(wp) ::  u_int                  !< u-component
     1984    REAL(wp) ::  v_int                  !< v-component
     1985    REAL(wp) ::  vpt_surface            !< near-surface virtual potential temperature
     1986    REAL(wp) ::  w_convective           !< convective velocity scale
     1987    REAL(wp) ::  z0_mean                !< mean roughness length
     1988    REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
     1989   
     1990!
     1991!-- Mean friction velocity and velocity scale. Therefore,
     1992!-- pre-calculate mean roughness length and surface sensible heat flux
     1993!-- in the model domain, which are further used to estimate friction
     1994!-- velocity and velocity scale. Note, for z0 linear averaging is applied,
     1995!-- even though this is known to unestimate the effective roughness.
     1996!-- This need to be revised later.
     1997    z0_mean_l  = 0.0_wp
     1998    shf_mean_l = 0.0_wp
     1999    DO  m = 1, surf_def_h(0)%ns
     2000       z0_mean_l  = z0_mean_l  + surf_def_h(0)%z0(m)
     2001       shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m)
     2002    ENDDO   
     2003    DO  m = 1, surf_lsm_h%ns
     2004       z0_mean_l  = z0_mean_l  + surf_lsm_h%z0(m)
     2005       shf_mean_l = shf_mean_l + surf_lsm_h%shf(m)
     2006    ENDDO
     2007    DO  m = 1, surf_usm_h%ns
     2008       z0_mean_l  = z0_mean_l  + surf_usm_h%z0(m)
     2009       shf_mean_l = shf_mean_l + surf_usm_h%shf(m)
     2010    ENDDO
     2011   
     2012#if defined( __parallel )
     2013    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
     2014                        comm2d, ierr )
     2015    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
     2016                        comm2d, ierr )
     2017#else
     2018    z0_mean  = z0_mean_l
     2019    shf_mean = shf_mean_l
     2020#endif
     2021    z0_mean  = z0_mean  / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
     2022    shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
     2023       
     2024!
     2025!-- Note, Inifor does not use logarithmic interpolation of the
     2026!-- velocity components near the ground, so that near-surface
     2027!-- wind speed and thus the friction velocity is overestimated.
     2028!-- However, friction velocity is used for turbulence
     2029!-- parametrization, so that more physically meaningful values are important.
     2030!-- Hence, derive friction velocity from wind speed at a reference height.
     2031!-- For a first guess use 20 m, which is in the range of the first
     2032!-- COSMO vertical level.
     2033    k_ref = MINLOC( ABS( zu - cosmo_ref ), DIM = 1 ) - 1
     2034   
     2035    friction_vel_l = 0.0_wp
     2036    IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
     2037       
     2038       i = MERGE( -1, nxr + 1, bc_dirichlet_l )
     2039
     2040       DO  j = nys, nyn
     2041!
     2042!--       Note, in u- and v- component the imposed perturbations
     2043!--       from the STG are already included. Check whether this
     2044!--       makes any difference compared to using the pure-mean
     2045!--       inflow profiles.
     2046          u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l )
     2047          v_int = v(k_ref,j,i)
     2048!
     2049!--       Calculate friction velocity and sum-up. Therefore, assume
     2050!--       neutral condtions.
     2051          friction_vel_l = friction_vel_l + kappa *                            &
     2052                            SQRT( u_int * u_int + v_int * v_int )  /           &
     2053                            LOG( zu(k_ref) / z0_mean )
     2054         
     2055       ENDDO
     2056       
     2057    ENDIF
     2058   
     2059    IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
     2060       
     2061       j = MERGE( -1, nyn + 1, bc_dirichlet_s )
     2062   
     2063       DO  i = nxl, nxr
     2064         
     2065          u_int = u(k_ref,j,i)
     2066          v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s )
     2067 
     2068          friction_vel_l = friction_vel_l + kappa *                            &
     2069                            SQRT( u_int * u_int + v_int * v_int )  /           &
     2070                            LOG( zu(k_ref) / z0_mean )                       
     2071         
     2072       ENDDO
     2073       
     2074    ENDIF
     2075   
     2076#if defined( __parallel )
     2077    CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM,        &
     2078                        comm2d, ierr )
     2079#else
     2080    scale_us = friction_vel_l
     2081#endif
     2082    scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp )
     2083   
     2084!
     2085!-- Compute mean Obukhov length
     2086    IF ( shf_mean > 0.0_wp )  THEN
     2087       scale_l = - pt_surface / ( kappa * g ) * scale_us**3 / shf_mean
     2088    ELSE
     2089       scale_l = 0.0_wp
     2090    ENDIF
     2091   
     2092!
     2093!-- Compute mean convective velocity scale. Note, in case the mean heat flux
     2094!-- is negative, set convective velocity scale to zero.
     2095    IF ( shf_mean > 0.0_wp )  THEN
     2096       w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
     2097    ELSE
     2098       w_convective = 0.0_wp
     2099    ENDIF
     2100!
     2101!-- Finally, in order to consider also neutral or stable stratification,
     2102!-- compute momentum velocity scale from u* and convective velocity scale,
     2103!-- according to Rotach et al. (1996).
     2104    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
     2105   
     2106 END SUBROUTINE calc_scaling_variables
     2107
    17142108 END MODULE synthetic_turbulence_generator_mod
  • palm/trunk/SOURCE/time_integration.f90

    r3298 r3347  
    2525! -----------------
    2626! $Id$
     27! - offline nesting separated from large-scale forcing module
     28! - changes for synthetic turbulence generator
     29!
     30! 3343 2018-10-15 10:38:52Z suehring
    2731! - Formatting, clean-up, comments (kanani)
    2832! - Added CALL to chem_emissions_setup (Russo)
     
    459463
    460464    USE lsf_nudging_mod,                                                       &
    461         ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref,       &
    462                lsf_nesting_offline, lsf_nesting_offline_mass_conservation
     465        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref
    463466
    464467    USE multi_agent_system_mod,                                                &
    465468        ONLY:  agents_active, multi_agent_system
    466 
     469       
     470    USE nesting_offl_mod,                                                      &
     471        ONLY:  nesting_offl_bc, nesting_offl_mass_conservation
     472       
    467473    USE netcdf_data_input_mod,                                                 &
    468         ONLY:  chem_emis, chem_emis_att, nest_offl, netcdf_data_input_lsf
     474        ONLY:  chem_emis, chem_emis_att, nest_offl,                            &
     475               netcdf_data_input_offline_nesting
    469476
    470477    USE ocean_mod,                                                             &
     
    515522
    516523    USE synthetic_turbulence_generator_mod,                                    &
    517         ONLY:  stg_main, use_syn_turb_gen
     524        ONLY:  dt_stg_call, dt_stg_adjust, parametrize_inflow_turbulence,      &
     525               stg_adjust, stg_main, time_stg_adjust, time_stg_call,           &
     526               use_syn_turb_gen
    518527
    519528    USE user_actions_mod,                                                      &
     
    638647       IF ( nesting_offline )  THEN
    639648          IF ( nest_offl%time(nest_offl%tind_p) <= time_since_reference_point )&
    640              CALL netcdf_data_input_lsf
     649             CALL netcdf_data_input_offline_nesting
    641650       ENDIF
    642651
     
    918927!--       Impose a turbulent inflow using the recycling method
    919928          IF ( turbulent_inflow )  CALL  inflow_turbulence
    920 
    921 !
    922 !--       Impose a turbulent inflow using synthetic generated turbulence
    923           IF ( use_syn_turb_gen )  CALL  stg_main
    924929
    925930!
     
    960965!--       boundaries.
    961966          IF ( nesting_offline  .AND.  intermediate_timestep_count ==          &
    962                                        intermediate_timestep_count_max  )  THEN
    963              CALL lsf_nesting_offline
    964 !
    965 !--          Moreover, ensure mass conservation
    966              CALL lsf_nesting_offline_mass_conservation
    967           ENDIF
    968 
     967                                       intermediate_timestep_count_max  )      &
     968             CALL nesting_offl_bc
     969!
     970!--       Impose a turbulent inflow using synthetic generated turbulence,
     971!--       only once per time step.
     972          IF ( use_syn_turb_gen  .AND.  time_stg_call >= dt_stg_call  .AND.    &
     973             intermediate_timestep_count == intermediate_timestep_count_max )  THEN! &
     974             CALL stg_main
     975          ENDIF
     976!
     977!--       Ensure mass conservation. This need to be done after imposing
     978!--       synthetic turbulence and top boundary condition for pressure is set to
     979!--       Neumann conditions.
     980!--       Is this also required in case of Dirichlet?
     981          IF ( nesting_offline )  CALL nesting_offl_mass_conservation
    969982!
    970983!--       Reduce the velocity divergence via the equation for perturbation
     
    12131226       time_dopr_listing          = time_dopr_listing        + dt_3d
    12141227       time_run_control   = time_run_control + dt_3d
     1228!
     1229!--    In case of synthetic turbulence generation and parametrized turbulence
     1230!--    information, update the time counter and if required, adjust the
     1231!--    STG to new atmospheric conditions.
     1232       IF ( use_syn_turb_gen  )  THEN
     1233          IF ( parametrize_inflow_turbulence )  THEN
     1234             time_stg_adjust = time_stg_adjust + dt_3d
     1235             IF ( time_stg_adjust >= dt_stg_adjust )  CALL stg_adjust
     1236          ENDIF
     1237          time_stg_call = time_stg_call + dt_3d
     1238       ENDIF
    12151239
    12161240!
     
    14721496
    14731497    ENDDO   ! time loop
    1474 
     1498   
     1499!
    14751500!-- Vertical nesting: Deallocate variables initialized for vertical nesting   
    14761501    IF ( vnest_init )  CALL vnest_deallocate
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3337 r3347  
    2828! -----------------
    2929! $Id$
     30! Enable USM initialization with default building parameters in case no static
     31! input file exist.
     32!
     33! 3343 2018-10-15 10:38:52Z suehring
    3034! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
    3135!
     
    69026906        INTEGER(iwp)                                          :: wtc
    69036907        REAL(wp), DIMENSION(n_surface_params)                 :: wtp
     6908       
     6909        LOGICAL                                               :: ascii_file = .FALSE.
    69046910   
    69056911        INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg)   :: usm_par
     
    69226928        IF ( building_type_f%from_file  .OR.  building_pars_f%from_file )      &
    69236929           RETURN
     6930!
     6931!--     Check if ASCII input file exists. If not, return and initialize USM
     6932!--     with default settings.
     6933        INQUIRE( FILE = 'SURFACE_PARAMETERS' // coupling_char,                 &
     6934                 EXIST = ascii_file )
     6935                 
     6936        IF ( .NOT. ascii_file )  RETURN
    69246937
    69256938!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracChangeset for help on using the changeset viewer.