Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 2938 for palm


Ignore:
Timestamp:
Mar 27, 2018 3:52:42 PM (7 years ago)
Author:
suehring
Message:

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

Location:
palm/trunk/SOURCE
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r2934 r2938  
    2525# -----------------
    2626# $Id$
     27# No initialization of child domains via dynamic input file, except for soil
     28# moisture and temperature
     29# Apply turbulence generator at non-cyclic lateral boundary in nesting case
     30#
     31# 2936 2018-03-27 14:49:27Z suehring
    2732# Added dependencies for parent and child synchronization
    2833#
    … …  
    10081013        modules.o \
    10091014        netcdf_data_input_mod.o \
     1015        pmc_interface_mod.o \
    10101016        radiation_model_mod.o \
    10111017        surface_mod.o
    … …  
    14341440        cpulog_mod.o \
    14351441        mod_kinds.o \
    1436         modules.o
     1442        modules.o \
     1443        pmc_interface_mod.o
    14371444temperton_fft_mod.o: \
    14381445        mod_kinds.o \
    … …  
    15101517        modules.o \
    15111518        plant_canopy_model_mod.o \
     1519        pmc_interface_mod.o \
    15121520        surface_mod.o \
    15131521        user_actions.o

  • palm/trunk/SOURCE/boundary_conds.f90

    r2766 r2938  
    2525! -----------------
    2626! $Id$
     27! Set boundary condition for TKE and TKE dissipation rate in case of nesting
     28! and if parent model operates in RANS mode but child model in LES mode.
     29! mode
     30!
     31! 2793 2018-02-07 10:54:33Z suehring
    2732! Removed preprocessor directive __chem
    2833!
    … …  
    193198               inflow_s, intermediate_timestep_count,                          &
    194199               microphysics_morrison, microphysics_seifert, nest_domain,       &
    195                nest_bound_l, nest_bound_s, nudging, ocean, outflow_l,          &
    196                outflow_n, outflow_r, outflow_s, passive_scalar, rans_tke_e,    &
    197                tsc, use_cmax
     200               nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, nudging,&
     201               ocean, outflow_l, outflow_n, outflow_r, outflow_s,              &
     202               passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
    198203
    199204    USE grid_variables,                                                        &
    … …  
    209214
    210215    USE pmc_interface,                                                         &
    211         ONLY : nesting_mode
     216        ONLY : nesting_mode, rans_mode_parent
    212217
    213218    USE surface_mod,                                                           &
    … …  
    321326
    322327!
    323 !-- Boundary conditions for TKE
    324 !-- Generally Neumann conditions with de/dz=0 are assumed
     328!-- Boundary conditions for TKE.
     329!-- Generally Neumann conditions with de/dz=0 are assumed.
    325330    IF ( .NOT. constant_diffusion )  THEN
    326331
    … …  
    343348       IF ( .NOT. nest_domain )  THEN
    344349          e_p(nzt+1,:,:) = e_p(nzt,:,:)
    345        ENDIF
    346     ENDIF
    347 
    348 !
    349 !-- Boundary conditions for TKE dissipation rate
     350!
     351!--    Nesting case: if parent operates in RANS mode and child in LES mode,
     352!--    no TKE is transfered. This case, set Neumann conditions at lateral and
     353!--    top child boundaries.
     354!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
     355!--    is treated in the nesting.
     356       ELSE
     357
     358          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
     359
     360
     361
     362             e_p(nzt+1,:,:) = e_p(nzt,:,:)
     363             IF ( nest_bound_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
     364             IF ( nest_bound_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
     365             IF ( nest_bound_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
     366             IF ( nest_bound_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
     367
     368          ENDIF
     369       ENDIF
     370    ENDIF
     371
     372!
     373!-- Boundary conditions for TKE dissipation rate.
    350374    IF ( rans_tke_e .AND. .NOT. nest_domain )  THEN
    351375       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)

  • palm/trunk/SOURCE/check_parameters.f90

    r2934 r2938  
    2525! -----------------
    2626! $Id$
     27! - Revise start and end indices for imposing random disturbances in case of
     28!   nesting or large-scale forcing.
     29! - Remove check for inifor initialization in case of nested runs.
     30! - Adapt call to check for synthetic turbulence geneartor settings.
     31!
     32! 2936 2018-03-27 14:49:27Z suehring
    2733! Check spinup in case of nested runs, spinup time and timestep must be
    2834! identical to assure synchronuous simulations.
    … …  
    13821388    ENDIF
    13831389!
    1384 !-- In case of nested run assure that all domains are initialized the same
    1385 !-- way, i.e. if at least at one domain is initialized with soil and
    1386 !-- atmospheric data provided by COSMO, all domains must be initialized the
    1387 !-- same way, to assure that soil and atmospheric quantities are
    1388 !-- consistent.
    1389     IF ( nested_run )  THEN
    1390        check_nest = .TRUE.
    1391 #if defined( __parallel )
    1392        CALL MPI_ALLREDUCE( TRIM( initializing_actions ) == 'inifor',           &
    1393                            check_nest, 1, MPI_LOGICAL,                         &
    1394                            MPI_LAND, MPI_COMM_WORLD, ierr )
    1395 
    1396        IF ( TRIM( initializing_actions ) == 'inifor'  .AND.                    &
    1397             .NOT.  check_nest )  THEN
    1398           message_string = 'In case of nesting, if at least in one ' //        &
    1399                            'domain initializing_actions = inifor, '  //        &
    1400                            'all domains need to be initialized that way.'
    1401           CALL message( 'netcdf_data_input_mod', 'PA0430', 3, 2, 0, 6, 0 )
    1402        ENDIF
    1403 #endif
    1404     ENDIF
    1405 !
    14061390!-- In case of spinup and nested run, spinup end time must be identical
    14071391!-- in order to have synchronously running simulations.
    … …  
    14361420
    14371421!
    1438 !-- When synthetic turbulence generator is used, peform addtional checks
    1439     IF ( syn_turb_gen )  CALL stg_check_parameters
     1422!-- Check for synthetic turbulence generator settings
     1423    CALL stg_check_parameters
    14401424!
    14411425!-- When plant canopy model is used, peform addtional checks
    … …  
    38983882          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
    38993883       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
    3900           dist_nxl = MAX( inflow_disturbance_begin, nxl )
     3884          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
     3885          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
    39013886       ENDIF
    39023887       IF ( bc_ns == 'dirichlet/radiation' )  THEN
    … …  
    39073892          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
    39083893       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
    3909           dist_nys = MAX( inflow_disturbance_begin, nys )
     3894          dist_nys    = MAX( inflow_disturbance_begin, nys )
     3895          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
    39103896       ENDIF
    39113897    ELSE
    … …  
    39183904          dist_nxl    = inflow_disturbance_begin
    39193905          dist_nxr(1) = inflow_disturbance_end
     3906       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
     3907          dist_nxr    = nx - inflow_disturbance_begin
     3908          dist_nxl    = inflow_disturbance_begin
    39203909       ENDIF
    39213910       IF ( bc_ns == 'dirichlet/radiation' )  THEN
    … …  
    39253914          dist_nys    = inflow_disturbance_begin
    39263915          dist_nyn(1) = inflow_disturbance_end
     3916       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
     3917          dist_nyn    = ny - inflow_disturbance_begin
     3918          dist_nys    = inflow_disturbance_begin
    39273919       ENDIF
    39283920    ENDIF

  • palm/trunk/SOURCE/init_3d_model.f90

    r2934 r2938  
    2525! -----------------
    2626! $Id$
     27! - Revise Inifor initialization for geostrophic wind components
     28! - Initialize synthetic turbulence generator in case of Inifor initialization 
     29!
     30! 2936 2018-03-27 14:49:27Z suehring
    2731! Synchronize parent and child models after initialization.
    2832! Remove obsolete masking of topography grid points for Runge-Kutta weighted
    … …  
    11891193!--           are not identical to the initial wind profiles in this case.
    11901194!--           This need to be further revised.
    1191 !           ug(:) = u_init(:)
    1192 !           vg(:) = v_init(:)
     1195          IF ( init_3d%from_file_ug )  THEN
     1196             ug(:) = init_3d%ug_init(:)
     1197          ENDIF
     1198          IF ( init_3d%from_file_vg )  THEN
     1199             vg(:) = init_3d%vg_init(:)
     1200          ENDIF
     1201
     1202          ug(nzt+1) = ug(nzt)
     1203          vg(nzt+1) = vg(nzt)
     1204
    11931205!
    11941206!--       Set inital w to 0
    … …  
    12431255!--       fluxes, etc.
    12441256          CALL init_surfaces
     1257!
     1258!--       Initialize turbulence generator
     1259          IF( use_syn_turb_gen )  CALL stg_init
    12451260
    12461261          CALL location_message( 'finished', .TRUE. )
    … …  
    13351350!
    13361351!--       Overwrite initial profiles in case of synthetic turbulence generator
    1337           IF( use_syn_turb_gen ) THEN
    1338              CALL stg_init
    1339           ENDIF
     1352          IF( use_syn_turb_gen )  CALL stg_init
    13401353
    13411354!

  • palm/trunk/SOURCE/init_pegrid.f90

    r2776 r2938  
    2525! -----------------
    2626! $Id$
     27! - No checks for domain decomposition in case of turbulence generator
     28!  (is done in stg module)
     29! - Introduce ids to indicate lateral processors for turbulence generator
     30!
     31! 2936 2018-03-27 14:49:27Z suehring
    2732! Variable use_synthetic_turbulence_generator has been abbreviated
    2833!
    … …  
    251256
    252257    USE synthetic_turbulence_generator_mod,                                    &
    253         ONLY:  use_syn_turb_gen
     258        ONLY:  id_stg_left, id_stg_north, id_stg_right, id_stg_south,          &
     259               use_syn_turb_gen
    254260
    255261    USE transpose_indices,                                                     &
    … …  
    267273    INTEGER(iwp) ::  id_outflow_source_l      !< local value of id_outflow_source
    268274    INTEGER(iwp) ::  id_recycling_l           !<
     275    INTEGER(iwp) ::  id_stg_left_l            !< left lateral boundary local core id in case of turbulence generator 
     276    INTEGER(iwp) ::  id_stg_north_l           !< north lateral boundary local core id in case of turbulence generator 
     277    INTEGER(iwp) ::  id_stg_right_l           !< right lateral boundary local core id in case of turbulence generator 
     278    INTEGER(iwp) ::  id_stg_south_l           !< south lateral boundary local core id in case of turbulence generator 
    269279    INTEGER(iwp) ::  ind(5)                   !<
    270280    INTEGER(iwp) ::  j                        !<
    … …  
    518528!-- 1. transposition  z --> x
    519529!-- This transposition is not neccessary in case of a 1d-decomposition along x
    520     IF ( psolver == 'poisfft'  .OR.  calculate_spectra  .OR.                   &
    521          use_syn_turb_gen )  THEN
     530    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
    522531
    523532       IF ( pdims(2) /= 1 )  THEN
    … …  
    12631272       ENDIF
    12641273    ENDIF
    1265 
     1274!
     1275!-- In case of synthetic turbulence geneartor determine ids.
     1276!-- Please note, if no forcing or nesting is applied, the generator is applied
     1277!-- only at the left lateral boundary.
     1278    IF ( use_syn_turb_gen )  THEN
     1279       IF ( force_bound_l  .OR.  nest_bound_l  .OR.  inflow_l )  THEN
     1280          id_stg_left_l = myidx
     1281       ELSE
     1282          id_stg_left_l = 0
     1283       ENDIF
     1284       IF ( force_bound_r  .OR.  nest_bound_r )  THEN
     1285          id_stg_right_l = myidx
     1286       ELSE
     1287          id_stg_right_l = 0
     1288       ENDIF
     1289       IF ( force_bound_s  .OR.  nest_bound_s )  THEN
     1290          id_stg_south_l = myidy
     1291       ELSE
     1292          id_stg_south_l = 0
     1293       ENDIF
     1294       IF ( force_bound_n  .OR.  nest_bound_n )  THEN
     1295          id_stg_north_l = myidy
     1296       ELSE
     1297          id_stg_north_l = 0
     1298       ENDIF
     1299
     1300       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1301       CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left,   1, MPI_INTEGER,       &
     1302                           MPI_SUM, comm1dx, ierr )
     1303
     1304       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1305       CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER,       &
     1306                           MPI_SUM, comm1dx, ierr )
     1307
     1308       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1309       CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER,       &
     1310                           MPI_SUM, comm1dy, ierr )
     1311
     1312       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1313       CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER,       &
     1314                           MPI_SUM, comm1dy, ierr )
     1315
     1316    ENDIF
     1317 
    12661318!
    12671319!-- Broadcast the id of the inflow PE

  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2932 r2938  
    2525! -----------------
    2626! $Id$
     27! Initialization of soil moisture and temperature via Inifor-provided data also
     28! in nested child domains, even if no dynamic input file is available for
     29! child domain. 1D soil profiles are received from parent model. 
     30!
     31! 2936 2018-03-27 14:49:27Z suehring
    2732! renamed lsm_par to land_surface_parameters. Bugfix in message calls
    2833!
    … …  
    23652370       USE control_parameters,                                                 &
    23662371           ONLY:  message_string
     2372
     2373       USE indices,                                                            &
     2374           ONLY:  nx, ny
     2375
     2376       USE pmc_interface,                                                      &
     2377           ONLY:  nested_run
    23672378   
    23682379       IMPLICIT NONE
     2380
     2381       LOGICAL      ::  init_soil_dynamically_in_child !< flag controlling initialization of soil in child domains
    23692382
    23702383       INTEGER(iwp) ::  i                       !< running index
    … …  
    23842397
    23852398       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
     2399       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
    23862400
    23872401!
    … …  
    38943908       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    38953909!
     3910!--       In case of nested runs, check if soil is initialized via dynamic
     3911!--       input file in root domain and distribute this information
     3912!--       to all embedded child domains. This case, soil moisture and
     3913!--       temperature will be initialized via root domain. 
     3914          init_soil_dynamically_in_child = .FALSE.
     3915          IF ( nested_run )  THEN
     3916#if defined( __parallel )
     3917             CALL MPI_ALLREDUCE( init_3d%from_file_tsoil  .OR.                 &
     3918                                 init_3d%from_file_msoil,                      &
     3919                                 init_soil_dynamically_in_child,               &
     3920                                 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
     3921#endif
     3922          ENDIF
     3923
     3924!
    38963925!--       First, initialize soil temperature and moisture.
    38973926!--       According to the initialization for surface and soil parameters,
    … …  
    39233952          ENDDO
    39243953!
    3925 !--       Level 2, if soil moisture and/or temperature  are
    3926 !--       provided from file, interpolate / extrapolate the provided data
    3927 !--       onto the respective soil layers. Please note, both, zs as well as
    3928 !--       init_3d%z_soil indicate a depth with positive values, so that no
    3929 !--       distinction between atmosphere is required concerning interpolation.
    3930 !--       Start with soil moisture
     3954!--       Initialization of soil moisture and temperature from file.
     3955!--       In case of nested runs, only root parent reads dynamic input file.
     3956!--       This case, transfer respective soil information provide
     3957!--       by dynamic input file from root parent domain onto all other domains.
     3958          IF ( init_soil_dynamically_in_child )  THEN
     3959!
     3960!--          Child domains will be only initialized with horizontally
     3961!--          averaged soil profiles in parent domain (for sake of simplicity).
     3962!--          If required, average soil data on root parent domain before
     3963!--          distribute onto child domains.
     3964             IF ( init_3d%from_file_msoil  .AND.  init_3d%lod_msoil == 2 )     &
     3965             THEN
     3966                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
     3967
     3968                DO  k = 0, init_3d%nzs-1
     3969                   pr_soil_init(k) = SUM( init_3d%msoil(k,nys:nyn,nxl:nxr)  )
     3970                ENDDO
     3971!
     3972!--             Allocate 1D array for soil-moisture profile (will not be
     3973!--             allocated in lod==2 case).
     3974                ALLOCATE( init_3d%msoil_init(0:init_3d%nzs-1) )
     3975                init_3d%msoil_init = 0.0_wp
     3976#if defined( __parallel )
     3977                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%msoil_init(0),    &
     3978                                    SIZE(pr_soil_init),                        &
     3979                                    MPI_REAL, MPI_SUM, comm2d, ierr )
     3980#endif
     3981                init_3d%msoil_init = init_3d%msoil_init /                      &
     3982                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
     3983                DEALLOCATE( pr_soil_init )
     3984             ENDIF
     3985             IF ( init_3d%from_file_tsoil  .AND.  init_3d%lod_tsoil == 2 )  THEN
     3986                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
     3987
     3988                DO  k = 0, init_3d%nzs-1
     3989                   pr_soil_init(k) = SUM( init_3d%tsoil(k,nys:nyn,nxl:nxr)  )
     3990                ENDDO
     3991!
     3992!--             Allocate 1D array for soil-temperature profile (will not be
     3993!--             allocated in lod==2 case).
     3994                ALLOCATE( init_3d%tsoil_init(0:init_3d%nzs-1) )
     3995                init_3d%tsoil_init = 0.0_wp
     3996#if defined( __parallel )
     3997                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%tsoil_init(0),    &
     3998                                    SIZE(pr_soil_init),                        &
     3999                                    MPI_REAL, MPI_SUM, comm2d, ierr )
     4000#endif
     4001                init_3d%tsoil_init = init_3d%tsoil_init /                      &
     4002                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
     4003                DEALLOCATE( pr_soil_init )
     4004
     4005             ENDIF
     4006
     4007#if defined( __parallel )
     4008!
     4009!--          Distribute soil grid information on file from root to all childs.
     4010!--          Only process with rank 0 sends the information.
     4011             CALL MPI_BCAST( init_3d%nzs,    1,                                &
     4012                             MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
     4013
     4014             IF ( .NOT.  ALLOCATED( init_3d%z_soil ) )                         &
     4015                ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
     4016
     4017             CALL MPI_BCAST( init_3d%z_soil, SIZE(init_3d%z_soil),             &
     4018                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
     4019#endif
     4020!
     4021!--          ALLOCATE arrays on child domains and set control attributes.
     4022             IF ( .NOT. init_3d%from_file_msoil )  THEN
     4023                ALLOCATE( init_3d%msoil_init(0:init_3d%nzs-1) )
     4024                init_3d%lod_msoil = 1
     4025                init_3d%from_file_msoil = .TRUE.
     4026             ENDIF
     4027             IF ( .NOT. init_3d%from_file_tsoil )  THEN
     4028                ALLOCATE( init_3d%tsoil_init(0:init_3d%nzs-1) )
     4029                init_3d%lod_tsoil = 1
     4030                init_3d%from_file_tsoil = .TRUE.
     4031             ENDIF
     4032
     4033
     4034#if defined( __parallel )
     4035!
     4036!--          Distribute soil profiles from root to all childs
     4037             CALL MPI_BCAST( init_3d%msoil_init, SIZE(init_3d%msoil_init),     &
     4038                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
     4039             CALL MPI_BCAST( init_3d%tsoil_init, SIZE(init_3d%tsoil_init),     &
     4040                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
     4041#endif
     4042
     4043          ENDIF
     4044!
     4045!--       Proceed with Level 2 initialization. Information from dynamic input
     4046!--       is now available on all processes.
    39314047          IF ( init_3d%from_file_msoil )  THEN
    39324048
    … …  
    39834099                ENDDO
    39844100             ENDIF
    3985 
    39864101          ENDIF
    39874102!

  • palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90

    r2718 r2938  
    2525! -----------------
    2626! $Id$
     27! Further improvements for nesting in larger-scale model
     28!
     29! 2863 2018-03-08 11:36:25Z suehring
    2730! Corrected "Former revisions" section
    2831!
    … …  
    210213       ENDIF
    211214!
    212 !--       Top boundary
     215!--    Top boundary
    213216       k = nzt
    214217       DO  i = nxl, nxr
    … …  
    235238          ENDDO
    236239       ENDDO
     240
     241write(9,*) "w correction", w_correct
     242flush(9)
    237243
    238244    END SUBROUTINE forcing_bc_mass_conservation
    … …  
    581587!--    identical to the Inifor grid. At the top boundary an extrapolation is
    582588!--    not possible.
    583        IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
    584           DO  i = nxlu, nxr
    585              DO  j = nys, nyn
    586                 u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *          &
    587                            ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *       &
    588                               MERGE( 1.0_wp, 0.0_wp,                           &
    589                                      BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
    590              ENDDO
    591           ENDDO
    592 
    593           DO  i = nxl, nxr
    594              DO  j = nysv, nyn
    595                 v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *          &
    596                            ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *       &
    597                               MERGE( 1.0_wp, 0.0_wp,                           &
    598                                      BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
    599              ENDDO
    600           ENDDO
    601 
     589       DO  i = nxlu, nxr
     590          DO  j = nys, nyn
     591             u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *             &
     592                        ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *          &
     593                           MERGE( 1.0_wp, 0.0_wp,                             &
     594                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
     595          ENDDO
     596       ENDDO
     597
     598       DO  i = nxl, nxr
     599          DO  j = nysv, nyn
     600             v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *             &
     601                        ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *          &
     602                           MERGE( 1.0_wp, 0.0_wp,                              &
     603                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
     604          ENDDO
     605       ENDDO
     606
     607       DO  i = nxl, nxr
     608          DO  j = nys, nyn
     609             w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *         &
     610                        ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *          &
     611                           MERGE( 1.0_wp, 0.0_wp,                              &
     612                                  BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
     613          ENDDO
     614       ENDDO
     615
     616
     617       IF ( .NOT. neutral )  THEN
    602618          DO  i = nxl, nxr
    603619             DO  j = nys, nyn
    604                 w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *      &
    605                            ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *       &
    606                               MERGE( 1.0_wp, 0.0_wp,                           &
    607                                      BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
    608              ENDDO
    609           ENDDO
    610 
    611 
    612           IF ( .NOT. neutral )  THEN
    613              DO  i = nxl, nxr
    614                 DO  j = nys, nyn
    615                    pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *     &
    616                            ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
    617                 ENDDO
    618              ENDDO
    619           ENDIF
    620 
    621           IF ( humidity )  THEN
    622              DO  i = nxl, nxr
    623                 DO  j = nys, nyn
    624                    q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *       &
    625                            ( force%q_top(1,j,i) - force%q_top(0,j,i) )
    626                 ENDDO
    627              ENDDO
    628           ENDIF
     620                pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *        &
     621                        ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
     622             ENDDO
     623          ENDDO
     624       ENDIF
     625
     626       IF ( humidity )  THEN
     627          DO  i = nxl, nxr
     628             DO  j = nys, nyn
     629                q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *          &
     630                        ( force%q_top(1,j,i) - force%q_top(0,j,i) )
     631             ENDDO
     632          ENDDO
     633       ENDIF
     634!
     635!--    At the edges( left-south, left-north, right-south and right-north) set
     636!--    data on ghost points.
     637       IF ( force_bound_l  .AND.  force_bound_s )  THEN
     638          DO  i = 1, nbgp
     639             u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
     640             w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
     641             IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
     642             IF ( humidity )       q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
     643          ENDDO
     644          DO  i = 1, nbgp+1
     645             v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
     646          ENDDO
     647       ENDIF
     648       IF ( force_bound_l  .AND.  force_bound_n )  THEN
     649          DO  i = 1, nbgp
     650             u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
     651             v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
     652             w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
     653             IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
     654             IF ( humidity )       q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
     655          ENDDO
     656       ENDIF
     657       IF ( force_bound_r  .AND.  force_bound_s )  THEN
     658          DO  i = 1, nbgp
     659             u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
     660             w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
     661             IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
     662             IF ( humidity )       q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
     663          ENDDO
     664          DO  i = 1, nbgp+1
     665             v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
     666          ENDDO
     667       ENDIF
     668       IF ( force_bound_r  .AND.  force_bound_n )  THEN
     669          DO  i = 1, nbgp
     670             u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
     671             v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
     672             w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
     673             IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
     674             IF ( humidity )       q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
     675          ENDDO
    629676       ENDIF
    630677!
    … …  
    941988    SUBROUTINE lsf_init
    942989
     990       USE control_parameters,                                                 &
     991           ONLY:  bc_lr_cyc, bc_ns_cyc
     992
    943993       USE netcdf_data_input_mod,                                              &
    944994           ONLY:  netcdf_data_input_lsf 
    … …  
    9751025
    9761026       IF ( forcing )  THEN
     1027!
     1028!--       Allocate arrays for geostrophic wind components. Arrays will
     1029!--       incorporate 2 time levels in order to interpolate in between. Please
     1030!--       note, forcing using geostrophic wind components is only required in
     1031!--       case of cyclic boundary conditions.
     1032          IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
     1033             ALLOCATE( force%ug(0:1,nzb:nzt+1) )
     1034             ALLOCATE( force%vg(0:1,nzb:nzt+1) )
     1035          ENDIF
    9771036!
    9781037!--       Allocate arrays for reading boundary values. Arrays will incorporate 2

  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r2930 r2938  
    2525! -----------------
    2626! $Id$
     27! Initial read of geostrophic wind components from dynamic driver.
     28!
     29! 2773 2018-01-30 14:12:54Z suehring
    2730! Revise checks for surface_fraction.
    2831!
    … …  
    126129       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos         !< vertical levels at w grid in dynamic input file
    127130
     131       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug         !< domain-averaged geostrophic component
     132       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg         !< domain-averaged geostrophic component
     133
    128134       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_left   !< u-component at left boundary
    129135       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_left   !< v-component at left boundary
    … …  
    180186       LOGICAL ::  from_file_tsoil  = .FALSE. !< flag indicating whether soil temperature is already initialized from file
    181187       LOGICAL ::  from_file_u      = .FALSE. !< flag indicating whether u is already initialized from file
     188       LOGICAL ::  from_file_ug     = .FALSE. !< flag indicating whether ug is already initialized from file
    182189       LOGICAL ::  from_file_v      = .FALSE. !< flag indicating whether v is already initialized from file
    183        LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
     190       LOGICAL ::  from_file_vg     = .FALSE. !< flag indicating whether ug is already initialized from file
     191       LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
     192
    184193
    185194       REAL(wp) ::  fill_msoil       !< fill value for soil moisture
    … …  
    422431       MODULE PROCEDURE get_variable_4d_real
    423432    END INTERFACE get_variable
     433
     434    INTERFACE get_variable_pr
     435       MODULE PROCEDURE get_variable_pr
     436    END INTERFACE get_variable_pr
    424437
    425438    INTERFACE get_attribute
    … …  
    20252038             ENDIF
    20262039!
    2027 !--          Read geostrophic wind components
    2028 !              IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
    2029 !
    2030 !              ENDIF
    2031 !              IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
    2032 !
    2033 !              ENDIF
    2034 
     2040!--          Read initial geostrophic wind components at
     2041!--          t = 0 (index 1 in file).
     2042             IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
     2043                ALLOCATE( init_3d%ug_init(nzb:nzt+1) )
     2044                CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', 1,          &
     2045                                      init_3d%ug_init )
     2046                init_3d%from_file_ug = .TRUE.
     2047             ELSE
     2048                init_3d%from_file_ug = .FALSE.
     2049             ENDIF
     2050             IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
     2051                ALLOCATE( init_3d%vg_init(nzb:nzt+1) )
     2052                CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', 1,          &
     2053                                      init_3d%vg_init )
     2054                init_3d%from_file_vg = .TRUE.
     2055             ELSE
     2056                init_3d%from_file_vg = .FALSE.
     2057             ENDIF
    20352058!
    20362059!--          Read inital 3D data of u, v, w, pt and q,
    … …  
    22222245!--          Read soil moisture
    22232246             IF ( land_surface )  THEN
     2247
    22242248                IF ( check_existence( var_names, 'init_soil_m' ) )  THEN
    22252249!
    … …  
    22402264!
    22412265!--                level-of-detail 2 - read 3D initialization data
    2242                    ELSEIF ( init_3d%lod_msoil == 2 )  THEN  ! need to be corrected
     2266                   ELSEIF ( init_3d%lod_msoil == 2 )  THEN
    22432267                      ALLOCATE ( init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
    22442268                      DO  i = nxl, nxr
    … …  
    22722296!
    22732297!--                level-of-detail 2 - read 3D initialization data
    2274                    ELSEIF ( init_3d%lod_tsoil == 2 )  THEN  ! need to be corrected
     2298                   ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
    22752299                      ALLOCATE ( init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
    22762300                      DO  i = nxl, nxr
    … …  
    23462370
    23472371       USE control_parameters,                                                 &
    2348            ONLY:  force_bound_l, force_bound_n, force_bound_r, force_bound_s,  &
     2372           ONLY:  bc_lr_cyc, bc_ns_cyc, force_bound_l, force_bound_n,          &
     2373                  force_bound_r, force_bound_s,                                &
    23492374                  forcing, humidity, message_string, neutral, simulated_time
     2375
    23502376
    23512377       USE indices,                                                            &
    … …  
    24302456             force%tind = MINLOC( ABS( force%time - simulated_time ), DIM = 1 )&
    24312457                          - 1
    2432              force%tind_p = force%tind + 1 
     2458             force%tind_p = force%tind + 1
     2459!
     2460!--          Read geostrophic wind components. In case of forcing, this is only
     2461!--          required if cyclic boundary conditions are applied.
     2462             IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
     2463                DO  t = force%tind, force%tind_p
     2464                   CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', t+1,     &
     2465                                         force%ug(t-force%tind,:) )
     2466                   CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', t+1,     &
     2467                                         force%ug(t-force%tind,:) )
     2468                ENDDO
     2469             ENDIF 
    24332470!
    24342471!--          Read data at lateral and top boundaries. Please note, at left and
    … …  
    38323869    END SUBROUTINE get_variable_1d_real
    38333870
     3871
     3872!------------------------------------------------------------------------------!
     3873! Description:
     3874! ------------
     3875!> Reads a time-dependent 1D float variable from file.
     3876!------------------------------------------------------------------------------!
     3877    SUBROUTINE get_variable_pr( id, variable_name, t, var )
     3878#if defined( __netcdf )
     3879
     3880       USE pegrid
     3881
     3882       IMPLICIT NONE
     3883
     3884       CHARACTER(LEN=*)                      ::  variable_name    !< variable name
     3885
     3886       INTEGER(iwp), INTENT(IN)              ::  id               !< file id
     3887       INTEGER(iwp), DIMENSION(1:2)          ::  id_dim           !< dimension ids
     3888       INTEGER(iwp)                          ::  id_var           !< dimension id
     3889       INTEGER(iwp)                          ::  n_file           !< number of data-points in file along z dimension
     3890       INTEGER(iwp), INTENT(IN)              ::  t                !< timestep number
     3891
     3892       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
     3893
     3894!
     3895!--    First, inquire variable ID
     3896       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
     3897!
     3898!--    Inquire dimension size of vertical dimension
     3899       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
     3900       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = n_file )
     3901!
     3902!--    Read variable.
     3903       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
     3904                               start = (/ 1,      t     /),                    &
     3905                               count = (/ n_file, 1     /) )
     3906       CALL handle_error( 'get_variable_pr', 527 )
     3907
     3908#endif
     3909    END SUBROUTINE get_variable_pr
     3910
     3911
    38343912!------------------------------------------------------------------------------!
    38353913! Description:

  • palm/trunk/SOURCE/parin.f90

    r2932 r2938  
    2525! -----------------
    2626! $Id$
     27! Change initialization in case child domain should be initialized with Inifor.
     28!
     29! 2936 2018-03-27 14:49:27Z suehring
    2730! inipar renamed to initialization_parameters.
    2831! d3par renamed to runtime_parameters.
    … …  
    907910!--       are set properly. An exception is made in case of restart runs and
    908911!--       if user decides to do everything by its own.
    909 !--       MS: is this really necessary?
    910912          IF ( nest_domain  .AND.  .NOT. (                                     &
    911913               TRIM( initializing_actions ) == 'read_restart_data'  .OR.       &
    912                TRIM( initializing_actions ) == 'by_user'            .OR.       &
    913                TRIM( initializing_actions ) == 'inifor' ) )  THEN
     914               TRIM( initializing_actions ) == 'by_user' ) )  THEN
    914915             initializing_actions = 'set_constant_profiles'
    915916          ENDIF

  • palm/trunk/SOURCE/pmc_interface_mod.f90

    r2903 r2938  
    2525! -----------------
    2626! $Id$
     27! - Nesting for RANS mode implemented
     28!    - Interpolation of TKE onto child domain only if both parent and child are
     29!      either in LES mode or in RANS mode
     30!    - Interpolation of dissipation if both parent and child are in RANS mode
     31!      and TKE-epsilon closure is applied
     32!    - Enable anterpolation of TKE and dissipation rate in case parent and
     33!      child operate in RANS mode
     34!
     35! - Some unused variables removed from ONLY list
     36! - Some formatting adjustments for particle nesting
     37!
     38! 2936 2018-03-27 14:49:27Z suehring
    2739! Control logics improved to allow nesting also in cases with
    2840! constant_flux_layer = .F. or constant_diffusion = .T.
    … …  
    254266#if defined( __nopointer )
    255267    USE arrays_3d,                                                             &
    256         ONLY:  dzu, dzw, e, e_p, nc, nr, pt, pt_p, q, q_p, qc, qr, s, u, u_p,  &
     268        ONLY:  diss, dzu, dzw, e, e_p, nc, nr, pt, q, qc, qr, s, u, u_p,       &
    257269               v, v_p, w, w_p, zu, zw
    258270#else
    259271   USE arrays_3d,                                                              &
    260         ONLY:  dzu, dzw, e, e_p, e_1, e_2, nc, nc_2, nc_p, nr, nr_2, nr_p, pt, &
    261                pt_p, pt_1, pt_2, q, q_p, q_1, q_2, qc, qc_2, qr, qr_2, s, s_2, &
    262                u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
     272        ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2, &
     273               pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                      &
     274               u, u_p, u_2, v, v_p, v_2, w, w_p, w_2, zu, zw
    263275#endif
    264276
    … …  
    269281               microphysics_morrison, microphysics_seifert,                    &
    270282               nest_bound_l, nest_bound_r, nest_bound_s, nest_bound_n,         &
    271                nest_domain, neutral, passive_scalar, roughness_length,         &
    272                simulated_time, topography, volume_flow
     283               nest_domain, neutral, passive_scalar, rans_mode, rans_tke_e,    &
     284               roughness_length, simulated_time, topography, volume_flow
    273285
    274286    USE chem_modules,                                                          &
    … …  
    359371
    360372    LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
     373    LOGICAL       ::  rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode)
    361374
    362375    REAL(wp), SAVE ::  anterp_relax_length_l = -1.0_wp   !<
    … …  
    387400    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !<
    388401
     402    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< coarse grid array on child domain - dissipation rate
    389403    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec   !<
    390404    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc  !<
    … …  
    401415    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
    402416
    403 
    404     REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c   !< child coarse data array for chemical species
     417    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c !< coarse grid array on child domain - chemical species
    405418
    406419!
    … …  
    616629           anterp_relax_length_t, child_to_parent, comm_world_nesting,         &
    617630           cpl_id, nested_run, nesting_datatransfer_mode, nesting_mode,        &
    618            parent_to_child
     631           parent_to_child, rans_mode_parent
    619632
    620633    PUBLIC pmci_boundary_conds
    … …  
    804817!   Initialize the pmc parent
    805818    CALL pmc_parentinit
     819
    806820!
    807821!-- Corners of all children of the present parent
    … …  
    821835
    822836       child_id = pmc_parent_for_child(m)
    823        IF ( myid == 0 )  THEN       
     837
     838       IF ( myid == 0 )  THEN
    824839
    825840          CALL pmc_recv_from_child( child_id, val,  size(val),  0, 123, ierr )
    … …  
    932947          DEALLOCATE( cl_coord_x )
    933948          DEALLOCATE( cl_coord_y )
     949
     950!
     951!--       Send information about operating mode (LES or RANS) to child. This will be
     952!--       used to control TKE nesting and setting boundary conditions properly.
     953          CALL pmc_send_to_child( child_id, rans_mode, 1, 0, 19, ierr ) 
    934954!
    935955!--       Send coarse grid information to child
    … …  
    9901010          ENDIF
    9911011       ENDDO
     1012
    9921013       CALL pmc_s_setind_and_allocmem( child_id )
    9931014    ENDDO
    … …  
    11951216       CALL pmc_set_dataarray_name( 'coarse', 'v'  ,'fine', 'v',  ierr )
    11961217       CALL pmc_set_dataarray_name( 'coarse', 'w'  ,'fine', 'w',  ierr )
     1218!
     1219!--    Set data array name for TKE. Please note, nesting of TKE is actually
     1220!--    only done if both parent and child are in LES or in RANS mode. Due to
     1221!--    design of model coupler, however, data array names must be already
     1222!--    available at this point.
    11971223       CALL pmc_set_dataarray_name( 'coarse', 'e'  ,'fine', 'e',  ierr )
     1224!
     1225!--    Nesting of dissipation rate only if both parent and child are in RANS
     1226!--    mode and TKE-epsilo closure is applied. Please so also comment for TKE
     1227!--    above.
     1228       CALL pmc_set_dataarray_name( 'coarse', 'diss'  ,'fine', 'diss',  ierr )
    11981229
    11991230       IF ( .NOT. neutral )  THEN
    … …  
    12601291          CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
    12611292          CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
     1293
     1294          CALL pmc_recv_from_parent( rans_mode_parent, 1, 0, 19, ierr )
     1295!
    12621296!
    12631297!--       Receive Coarse grid information.
    … …  
    13191353       CALL MPI_BCAST( cg%zu, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
    13201354       CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
    1321        
     1355       CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )
     1356 
    13221357!
    13231358!--    Find the index bounds for the nest domain in the coarse-grid index space
    … …  
    13571392       ENDIF
    13581393!
    1359 !--    Define the SGS-TKE scaling factor based on the grid-spacing ratio
    1360        IF ( .NOT. constant_diffusion )  THEN
    1361           CALL pmci_init_tkefactor
    1362        ENDIF
     1394!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio. Only
     1395!--    if both parent and child are in LES mode or in RANS mode.
     1396!--    Please note, in case parent and child are in RANS mode, TKE weighting
     1397!--    factor is simply one.
     1398       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
     1399            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
     1400               .NOT. constant_diffusion ) )  CALL pmci_init_tkefactor
    13631401!
    13641402!--    Two-way coupling for general and vertical nesting.
    … …  
    29953033!--    energy spectrum. Near the surface, the reduction of TKE is made
    29963034!--    smaller than further away from the surface.
     3035!--    Please note, in case parent and child model operate in RANS mode,
     3036!--    TKE is not grid depenedent and weighting factor is one.
    29973037
    29983038       IMPLICIT NONE
    … …  
    30123052       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !<       
    30133053
    3014        IF ( nest_bound_l )  THEN
    3015           ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
    3016           tkefactor_l = 0.0_wp
    3017           i = nxl - 1
    3018           DO  j = nysg, nyng
    3019              k_wall = get_topography_top_index_ji( j, i, 's' )
    3020 
    3021              DO  k = k_wall + 1, nzt
    3022 
    3023                 kc     = kco(k) + 1
    3024                 glsf   = ( dx * dy * dzu(k) )**p13
    3025                 glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
    3026                 height = zu(k) - zu(k_wall)
    3027                 fw     = EXP( -cfw * height / glsf )
    3028                 tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
    3029                                               ( glsf / glsc )**p23 )
     3054!
     3055       IF ( .NOT. rans_mode  .AND.  .NOT. rans_mode_parent )  THEN
     3056          IF ( nest_bound_l )  THEN
     3057             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
     3058             tkefactor_l = 0.0_wp
     3059             i = nxl - 1
     3060             DO  j = nysg, nyng
     3061                k_wall = get_topography_top_index_ji( j, i, 's' )
     3062
     3063                DO  k = k_wall + 1, nzt
     3064                   kc     = kco(k) + 1
     3065                   glsf   = ( dx * dy * dzu(k) )**p13
     3066                   glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
     3067                   height = zu(k) - zu(k_wall)
     3068                   fw     = EXP( -cfw * height / glsf )
     3069                   tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
     3070                                                 ( glsf / glsc )**p23 )
     3071                ENDDO
     3072                tkefactor_l(k_wall,j) = c_tkef * fw0
    30303073             ENDDO
    3031              tkefactor_l(k_wall,j) = c_tkef * fw0
    3032           ENDDO
    3033        ENDIF
    3034 
    3035        IF ( nest_bound_r )  THEN
    3036           ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
    3037           tkefactor_r = 0.0_wp
    3038           i = nxr + 1
    3039           DO  j = nysg, nyng
    3040              k_wall = get_topography_top_index_ji( j, i, 's' )
    3041 
    3042              DO  k = k_wall + 1, nzt
     3074          ENDIF
     3075
     3076          IF ( nest_bound_r )  THEN
     3077             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
     3078             tkefactor_r = 0.0_wp
     3079             i = nxr + 1
     3080             DO  j = nysg, nyng
     3081                k_wall = get_topography_top_index_ji( j, i, 's' )
     3082
     3083                DO  k = k_wall + 1, nzt
     3084                   kc     = kco(k) + 1
     3085                   glsf   = ( dx * dy * dzu(k) )**p13
     3086                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
     3087                   height = zu(k) - zu(k_wall)
     3088                   fw     = EXP( -cfw * height / glsf )
     3089                   tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
     3090                                                 ( glsf / glsc )**p23 )
     3091                ENDDO
     3092                tkefactor_r(k_wall,j) = c_tkef * fw0
     3093             ENDDO
     3094          ENDIF
     3095
     3096          IF ( nest_bound_s )  THEN
     3097             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
     3098             tkefactor_s = 0.0_wp
     3099             j = nys - 1
     3100             DO  i = nxlg, nxrg
     3101                k_wall = get_topography_top_index_ji( j, i, 's' )
     3102               
     3103                DO  k = k_wall + 1, nzt
     3104   
     3105                   kc     = kco(k) + 1
     3106                   glsf   = ( dx * dy * dzu(k) )**p13
     3107                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
     3108                   height = zu(k) - zu(k_wall)
     3109                   fw     = EXP( -cfw*height / glsf )
     3110                   tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
     3111                        ( glsf / glsc )**p23 )
     3112                ENDDO
     3113                tkefactor_s(k_wall,i) = c_tkef * fw0
     3114             ENDDO
     3115          ENDIF
     3116
     3117          IF ( nest_bound_n )  THEN
     3118             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
     3119             tkefactor_n = 0.0_wp
     3120             j = nyn + 1
     3121             DO  i = nxlg, nxrg
     3122                k_wall = get_topography_top_index_ji( j, i, 's' )
     3123
     3124                DO  k = k_wall + 1, nzt
     3125
     3126                   kc     = kco(k) + 1
     3127                   glsf   = ( dx * dy * dzu(k) )**p13
     3128                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
     3129                   height = zu(k) - zu(k_wall)
     3130                   fw     = EXP( -cfw * height / glsf )
     3131                   tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
     3132                                                 ( glsf / glsc )**p23 )
     3133                ENDDO
     3134                tkefactor_n(k_wall,i) = c_tkef * fw0
     3135             ENDDO
     3136          ENDIF
     3137
     3138          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
     3139          k = nzt
     3140
     3141          DO  i = nxlg, nxrg
     3142             DO  j = nysg, nyng
     3143!
     3144!--             Determine vertical index for local topography top
     3145                k_wall = get_topography_top_index_ji( j, i, 's' )
    30433146
    30443147                kc     = kco(k) + 1
    … …  
    30473150                height = zu(k) - zu(k_wall)
    30483151                fw     = EXP( -cfw * height / glsf )
    3049                 tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
     3152                tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
    30503153                                              ( glsf / glsc )**p23 )
    30513154             ENDDO
    3052              tkefactor_r(k_wall,j) = c_tkef * fw0
    30533155          ENDDO
    3054        ENDIF
    3055 
    3056        IF ( nest_bound_s )  THEN
    3057           ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
    3058           tkefactor_s = 0.0_wp
    3059           j = nys - 1
    3060           DO  i = nxlg, nxrg
    3061              k_wall = get_topography_top_index_ji( j, i, 's' )
    3062              
    3063              DO  k = k_wall + 1, nzt
    3064  
    3065                 kc     = kco(k) + 1
    3066                 glsf   = ( dx * dy * dzu(k) )**p13
    3067                 glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
    3068                 height = zu(k) - zu(k_wall)
    3069                 fw     = EXP( -cfw*height / glsf )
    3070                 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
    3071                      ( glsf / glsc )**p23 )
    3072              ENDDO
    3073              tkefactor_s(k_wall,i) = c_tkef * fw0
    3074           ENDDO
    3075        ENDIF
    3076 
    3077        IF ( nest_bound_n )  THEN
    3078           ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
    3079           tkefactor_n = 0.0_wp
    3080           j = nyn + 1
    3081           DO  i = nxlg, nxrg
    3082              k_wall = get_topography_top_index_ji( j, i, 's' )
    3083 
    3084              DO  k = k_wall + 1, nzt
    3085 
    3086                 kc     = kco(k) + 1
    3087                 glsf   = ( dx * dy * dzu(k) )**p13
    3088                 glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
    3089                 height = zu(k) - zu(k_wall)
    3090                 fw     = EXP( -cfw * height / glsf )
    3091                 tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
    3092                                               ( glsf / glsc )**p23 )
    3093              ENDDO
    3094              tkefactor_n(k_wall,i) = c_tkef * fw0
    3095           ENDDO
    3096        ENDIF
    3097 
    3098        ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
    3099        k = nzt
    3100 
    3101        DO  i = nxlg, nxrg
    3102           DO  j = nysg, nyng
    3103 !
    3104 !--          Determine vertical index for local topography top
    3105              k_wall = get_topography_top_index_ji( j, i, 's' )
    3106 
    3107              kc     = kco(k) + 1
    3108              glsf   = ( dx * dy * dzu(k) )**p13
    3109              glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
    3110              height = zu(k) - zu(k_wall)
    3111              fw     = EXP( -cfw * height / glsf )
    3112              tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
    3113                                            ( glsf / glsc )**p23 )
    3114           ENDDO
    3115        ENDDO
     3156!
     3157!--    RANS mode
     3158       ELSE
     3159          IF ( nest_bound_l )  THEN
     3160             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
     3161             tkefactor_l = 1.0_wp
     3162          ENDIF
     3163          IF ( nest_bound_r )  THEN
     3164             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
     3165             tkefactor_r = 1.0_wp
     3166          ENDIF
     3167          IF ( nest_bound_s )  THEN
     3168             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
     3169             tkefactor_s = 1.0_wp
     3170          ENDIF
     3171          IF ( nest_bound_n )  THEN
     3172             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
     3173             tkefactor_n = 1.0_wp
     3174          ENDIF
     3175
     3176          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
     3177          tkefactor_t = 1.0_wp
     3178
     3179       ENDIF
    31163180     
    31173181    END SUBROUTINE pmci_init_tkefactor
    … …  
    31773241!-- List of array names, which can be coupled.
    31783242!-- In case of 3D please change also the second array for the pointer version
    3179     IF ( TRIM(name) == "u"  )  p_3d => u
    3180     IF ( TRIM(name) == "v"  )  p_3d => v
    3181     IF ( TRIM(name) == "w"  )  p_3d => w
    3182     IF ( TRIM(name) == "e"  )  p_3d => e
    3183     IF ( TRIM(name) == "pt" )  p_3d => pt
    3184     IF ( TRIM(name) == "q"  )  p_3d => q
    3185     IF ( TRIM(name) == "qc" )  p_3d => qc
    3186     IF ( TRIM(name) == "qr" )  p_3d => qr
    3187     IF ( TRIM(name) == "nr" )  p_3d => nr
    3188     IF ( TRIM(name) == "nc" )  p_3d => nc
    3189     IF ( TRIM(name) == "s"  )  p_3d => s
    3190     IF ( TRIM(name) == "nr_part"  )   i_2d => nr_part
    3191     IF ( TRIM(name) == "part_adr"  )  i_2d => part_adr
     3243    IF ( TRIM(name) == "u"          )  p_3d => u
     3244    IF ( TRIM(name) == "v"          )  p_3d => v
     3245    IF ( TRIM(name) == "w"          )  p_3d => w
     3246    IF ( TRIM(name) == "e"          )  p_3d => e
     3247    IF ( TRIM(name) == "pt"         )  p_3d => pt
     3248    IF ( TRIM(name) == "q"          )  p_3d => q
     3249    IF ( TRIM(name) == "qc"         )  p_3d => qc
     3250    IF ( TRIM(name) == "qr"         )  p_3d => qr
     3251    IF ( TRIM(name) == "nr"         )  p_3d => nr
     3252    IF ( TRIM(name) == "nc"         )  p_3d => nc
     3253    IF ( TRIM(name) == "s"          )  p_3d => s
     3254    IF ( TRIM(name) == "diss"       )  p_3d => diss   
     3255    IF ( TRIM(name) == "nr_part"    )   i_2d => nr_part
     3256    IF ( TRIM(name) == "part_adr"   )  i_2d => part_adr
    31923257    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d => chem_species(n)%conc
    31933258
    … …  
    32193284    ENDIF
    32203285#else
    3221     IF ( TRIM(name) == "u"  )  p_3d_sec => u_2
    3222     IF ( TRIM(name) == "v"  )  p_3d_sec => v_2
    3223     IF ( TRIM(name) == "w"  )  p_3d_sec => w_2
    3224     IF ( TRIM(name) == "e"  )  p_3d_sec => e_2
    3225     IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
    3226     IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
    3227     IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
    3228     IF ( TRIM(name) == "qr" )  p_3d_sec => qr_2
    3229     IF ( TRIM(name) == "nr" )  p_3d_sec => nr_2
    3230     IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
    3231     IF ( TRIM(name) == "s"  )  p_3d_sec => s_2
     3286    IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
     3287    IF ( TRIM(name) == "v"    )  p_3d_sec => v_2
     3288    IF ( TRIM(name) == "w"    )  p_3d_sec => w_2
     3289    IF ( TRIM(name) == "e"    )  p_3d_sec => e_2
     3290    IF ( TRIM(name) == "pt"   )  p_3d_sec => pt_2
     3291    IF ( TRIM(name) == "q"    )  p_3d_sec => q_2
     3292    IF ( TRIM(name) == "qc"   )  p_3d_sec => qc_2
     3293    IF ( TRIM(name) == "qr"   )  p_3d_sec => qr_2
     3294    IF ( TRIM(name) == "nr"   )  p_3d_sec => nr_2
     3295    IF ( TRIM(name) == "nc"   )  p_3d_sec => nc_2
     3296    IF ( TRIM(name) == "s"    )  p_3d_sec => s_2
     3297    IF ( TRIM(name) == "diss" )  p_3d_sec => diss_2
    32323298    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d_sec => spec_conc_2(:,:,:,n)
    32333299
    … …  
    33583424       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1,js:je,is:ie) )
    33593425       p_3d => ec
     3426    ELSEIF ( TRIM( name ) == "diss" )  THEN
     3427       IF ( .NOT. ALLOCATED( dissc ) )  ALLOCATE( dissc(0:nzc+1,js:je,is:ie) )
     3428       p_3d => dissc
    33603429    ELSEIF ( TRIM( name ) == "pt")  THEN
    33613430       IF ( .NOT. ALLOCATED( ptc ) )  ALLOCATE( ptc(0:nzc+1,js:je,is:ie) )
    … …  
    33793448       IF ( .NOT. ALLOCATED( sc ) )  ALLOCATE( sc(0:nzc+1,js:je,is:ie) )
    33803449       p_3d => sc
    3381     ELSEIF (trim(name) == "nr_part") then
    3382        IF (.not.allocated(nr_partc))  allocate(nr_partc(js:je, is:ie))
     3450    ELSEIF ( TRIM( name ) == "nr_part") THEN
     3451       IF ( .NOT. ALLOCATED( nr_partc ) )  ALLOCATE( nr_partc(js:je,is:ie) )
    33833452       i_2d => nr_partc
    3384     ELSEIF (trim(name) == "part_adr") then
    3385        IF (.not.allocated(part_adrc))  allocate(part_adrc(js:je, is:ie))
     3453    ELSEIF ( TRIM( name ) == "part_adr") THEN
     3454       IF ( .NOT. ALLOCATED( part_adrc ) )  ALLOCATE( part_adrc(js:je,is:ie) )
    33863455       i_2d => part_adrc
    33873456    ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
    … …  
    34843553                                   r2yo, r1zw, r2zw, 'w' )
    34853554
    3486        IF ( .NOT. constant_diffusion )  THEN
    3487           CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
    3488                                       r1yo, r2yo, r1zo, r2zo, 'e' )
    3489        ENDIF
    3490        
     3555       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.          &
     3556            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.           &
     3557               .NOT. constant_diffusion ) )  THEN
     3558          CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
     3559                                      r2yo, r1zo, r2zo, 'e' )
     3560       ENDIF
     3561
     3562       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     3563          CALL pmci_interp_tril_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
     3564                                      r1yo, r2yo, r1zo, r2zo, 's' )
     3565       ENDIF
     3566
    34913567       IF ( .NOT. neutral )  THEN
    34923568          CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
    … …  
    42074283                                       nzt_topo_nestbc_l, 'l', 'w' )
    42084284
    4209              IF ( .NOT. constant_diffusion )  THEN
     4285             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
     4286                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
     4287                     .NOT. constant_diffusion ) )  THEN
    42104288                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
    42114289                                          r1yo, r2yo, r1zo, r2zo,              &
    … …  
    42144292                                          nzt_topo_nestbc_l, 'l', 'e' )
    42154293             ENDIF
    4216              
     4294
     4295             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     4296                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
     4297                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
     4298                                          logc_w_l, logc_ratio_w_l,            &
     4299                                          logc_kbounds_w_l,                    &
     4300                                          nzt_topo_nestbc_l, 'l', 's' )
     4301             ENDIF
     4302
    42174303             IF ( .NOT. neutral )  THEN
    42184304                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
    … …  
    43044390                                       nzt_topo_nestbc_r, 'r', 'w' )
    43054391
    4306              IF ( .NOT. constant_diffusion )  THEN
     4392             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
     4393                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
     4394                     .NOT. constant_diffusion ) )  THEN
    43074395                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
    43084396                                          r1yo,r2yo, r1zo, r2zo,               &
    … …  
    43104398                                          logc_kbounds_w_r,                    &
    43114399                                          nzt_topo_nestbc_r, 'r', 'e' )
     4400
     4401             ENDIF
     4402
     4403             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     4404                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
     4405                                          r2xo, r1yo,r2yo, r1zo, r2zo,         &
     4406                                          logc_w_r, logc_ratio_w_r,            &
     4407                                          logc_kbounds_w_r,                    &
     4408                                          nzt_topo_nestbc_r, 'r', 's' )
     4409
    43124410             ENDIF
    43134411
    … …  
    43824480                ENDDO
    43834481             ENDIF
    4384 
    43854482          ENDIF
    43864483!
    … …  
    44064503                                       nzt_topo_nestbc_s, 's','w' )
    44074504
    4408              IF ( .NOT. constant_diffusion )  THEN
     4505             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
     4506                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
     4507                     .NOT. constant_diffusion ) )  THEN
    44094508                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
    44104509                                          r1yo, r2yo, r1zo, r2zo,              &
    … …  
    44124511                                          logc_kbounds_w_s,                    &
    44134512                                          nzt_topo_nestbc_s, 's', 'e' )
     4513
    44144514             ENDIF
    4415          
     4515
     4516             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     4517                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
     4518                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
     4519                                          logc_w_s, logc_ratio_w_s,            &
     4520                                          logc_kbounds_w_s,                    &
     4521                                          nzt_topo_nestbc_s, 's', 's' )
     4522
     4523             ENDIF
     4524
    44164525             IF ( .NOT. neutral )  THEN
    44174526                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
    … …  
    44824591                ENDDO
    44834592             ENDIF
    4484 
    44854593          ENDIF
    44864594!
    … …  
    45054613                                       logc_kbounds_w_n,                       &
    45064614                                       nzt_topo_nestbc_n, 'n', 'w' )
    4507              IF ( .NOT. constant_diffusion )  THEN
     4615
     4616             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
     4617                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
     4618                     .NOT. constant_diffusion ) )  THEN
    45084619                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
    45094620                                          r1yo, r2yo, r1zo, r2zo,              &
    … …  
    45114622                                          logc_kbounds_w_n,                    &
    45124623                                          nzt_topo_nestbc_n, 'n', 'e' )
     4624
    45134625             ENDIF
    4514          
     4626
     4627             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     4628                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
     4629                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
     4630                                          logc_w_n, logc_ratio_w_n,            &
     4631                                          logc_kbounds_w_n,                    &
     4632                                          nzt_topo_nestbc_n, 'n', 's' )
     4633
     4634             ENDIF
     4635
    45154636             IF ( .NOT. neutral )  THEN
    45164637                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
    … …  
    45814702                ENDDO
    45824703             ENDIF
    4583              
    45844704          ENDIF
    45854705
    … …  
    45934713       CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
    45944714                                r2yo, r1zw, r2zw, 'w' )
    4595        IF ( .NOT. constant_diffusion )  THEN
     4715
     4716       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
     4717            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
     4718               .NOT. constant_diffusion ) )  THEN
    45964719          CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
    45974720                                   r2yo, r1zo, r2zo, 'e' )
    45984721       ENDIF
     4722
     4723       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
     4724          CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
     4725                                   r1yo, r2yo, r1zo, r2zo, 's' )
     4726       ENDIF
     4727
    45994728       IF ( .NOT. neutral )  THEN
    46004729          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
    … …  
    46444773          ENDDO
    46454774       ENDIF
    4646    
     4775
    46474776   END SUBROUTINE pmci_interpolation
    46484777
    … …  
    46664795      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
    46674796                               kfuw, ijfc_s, kfc_w, 'w' )
     4797!
     4798!--   Anterpolation of TKE and dissipation rate if parent and child are in
     4799!--   RANS mode.
     4800      IF ( rans_mode_parent  .AND.  rans_mode )  THEN
     4801         CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
     4802                                  kfuo, ijfc_s, kfc_s, 'e' )
     4803!
     4804!--      Anterpolation of dissipation rate only if TKE-e closure is applied.
     4805         IF ( rans_tke_e )  THEN
     4806            CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
     4807                                     kflo, kfuo, ijfc_s, kfc_s, 'diss' )
     4808         ENDIF
     4809
     4810      ENDIF
    46684811
    46694812      IF ( .NOT. neutral )  THEN
    … …  
    55695712 END SUBROUTINE pmci_boundary_conds
    55705713
     5714
    55715715END MODULE pmc_interface

  • palm/trunk/SOURCE/pmc_mpi_wrapper_mod.f90

    r2841 r2938  
    2626! -----------------
    2727! $Id$
     28! Extent interface by logical buffer
     29!
     30! 2936 2018-03-27 14:49:27Z suehring
    2831! Bugfix: wrong placement of include 'mpif.h' corrected
    2932!
    … …  
    103106
    104107    INTERFACE pmc_recv_from_parent
     108       MODULE PROCEDURE pmc_recv_from_parent_logical
    105109       MODULE PROCEDURE pmc_recv_from_parent_integer
    106110       MODULE PROCEDURE pmc_recv_from_parent_real_r1
    … …  
    110114
    111115    INTERFACE pmc_send_to_child
     116       MODULE PROCEDURE pmc_send_to_child_logical
    112117       MODULE PROCEDURE pmc_send_to_child_integer
    113118       MODULE PROCEDURE pmc_send_to_child_real_r1
    … …  
    149154
    150155
     156 SUBROUTINE pmc_recv_from_parent_logical( buf, n, parent_rank, tag, ierr )
     157
     158    IMPLICIT NONE
     159
     160    INTEGER, INTENT(IN)                ::  n            !<
     161    INTEGER, INTENT(IN)                ::  parent_rank  !<
     162    INTEGER, INTENT(IN)                ::  tag          !<
     163    INTEGER, INTENT(OUT)               ::  ierr         !<
     164
     165    LOGICAL, INTENT(OUT)               ::  buf          !<
     166   
     167    ierr = 0
     168    CALL MPI_RECV( buf, n, MPI_LOGICAL, parent_rank, tag, m_to_parent_comm,    &
     169                   MPI_STATUS_IGNORE, ierr )
     170
     171 END SUBROUTINE pmc_recv_from_parent_logical
     172
     173 SUBROUTINE pmc_send_to_child_logical( child_id, buf, n, child_rank, tag,      &
     174                                       ierr )
     175
     176    IMPLICIT NONE
     177
     178    INTEGER, INTENT(IN)               ::  child_id     !<
     179    INTEGER, INTENT(IN)               ::  n            !<
     180    INTEGER, INTENT(IN)               ::  child_rank   !<
     181    INTEGER, INTENT(IN)               ::  tag          !<
     182    INTEGER, INTENT(OUT)              ::  ierr         !<
     183
     184    LOGICAL, INTENT(IN)               ::  buf          !<
     185   
     186    ierr = 0
     187    CALL MPI_SEND( buf, n, MPI_LOGICAL, child_rank, tag,                       &
     188                   m_to_child_comm(child_id), ierr )
     189
     190 END SUBROUTINE pmc_send_to_child_logical
     191
     192
    151193 SUBROUTINE pmc_send_to_parent_integer( buf, n, parent_rank, tag, ierr )
    152194
    … …  
    167209
    168210
    169 
    170211 SUBROUTINE pmc_recv_from_parent_integer( buf, n, parent_rank, tag, ierr )
    171212
    … …  
    312353
    313354 END SUBROUTINE pmc_recv_from_parent_real_r3
    314 
    315355
    316356

  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r2894 r2938  
    2525! -----------------
    2626! $Id$
     27! Apply turbulence generator at all non-cyclic lateral boundaries in case of
     28! realistic Inifor large-scale forcing or RANS-LES nesting
     29!
     30! 2936 2018-03-27 14:49:27Z suehring
    2731! variable named found has been introduced for checking if restart data was found,
    2832! reading of restart strings has been moved completely to read_restart_data_mod,
    … …  
    101105
    102106    USE control_parameters,                                                    &
    103         ONLY:  initializing_actions, message_string,                           &
    104                syn_turb_gen
     107        ONLY:  initializing_actions, message_string, syn_turb_gen
    105108
    106109    USE cpulog,                                                                &
    … …  
    108111
    109112    USE indices,                                                               &
    110         ONLY:  nbgp, nzb, nzt, nyng, nysg
     113        ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
    111114
    112115    USE kinds
    … …  
    117120
    118121    USE pegrid,                                                                &
    119         ONLY:  comm1dx, comm1dy, comm2d, id_inflow, ierr, myidx, pdims
     122        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
    120123
    121124    USE transpose_indices,                                                     &
    … …  
    132135    LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
    133136
    134     INTEGER(iwp) :: stg_type_yz        !< MPI type for full z range
    135     INTEGER(iwp) :: stg_type_yz_small  !< MPI type for small z range
    136     INTEGER(iwp) :: merg               !< maximum length scale (in gp)
    137     INTEGER(iwp) :: mergp              !< merg + nbgp
     137    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator 
     138    INTEGER(iwp) ::  id_stg_north       !< north lateral boundary core id in case of turbulence generator 
     139    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator 
     140    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator 
     141    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
     142    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
     143    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
     144    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
     145    INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
     146    INTEGER(iwp) ::  mergp              !< merg + nbgp
     147    INTEGER(iwp) ::  nzb_x_stg          !<
     148    INTEGER(iwp) ::  nzt_x_stg          !<
     149    INTEGER(iwp) ::  nzb_y_stg          !<
     150    INTEGER(iwp) ::  nzt_y_stg          !<
    138151
    139152    REAL(wp) :: mc_factor    !< mass flux correction factor
    140153
    141     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs      !< displacement for MPI_GATHERV
    142     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count  !< receive count for MPI_GATHERV
    143     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy         !< length scale of u in y direction (in gp)
    144     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz         !< length scale of u in z direction (in gp)
    145     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy         !< length scale of v in y direction (in gp)
    146     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz         !< length scale of v in z direction (in gp)
    147     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy         !< length scale of w in y direction (in gp)
    148     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz         !< length scale of w in z direction (in gp)
     154    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz      !< displacement for MPI_GATHERV
     155    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz  !< receive count for MPI_GATHERV
     156    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz      !< displacement for MPI_GATHERV
     157    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz  !< receive count for MPI_GATHERV
     158    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux            !< length scale of u in x direction (in gp)
     159    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy            !< length scale of u in y direction (in gp)
     160    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz            !< length scale of u in z direction (in gp)
     161    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx            !< length scale of v in x direction (in gp)
     162    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy            !< length scale of v in y direction (in gp)
     163    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz            !< length scale of v in z direction (in gp)
     164    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx            !< length scale of w in x direction (in gp)
     165    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy            !< length scale of w in y direction (in gp)
     166    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz            !< length scale of w in z direction (in gp)
    149167
    150168    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed        !< seed of random number for rn-generator
    … …  
    160178    REAL(wp), DIMENSION(:), ALLOCATABLE :: tw              !< Lagrangian time scale of w
    161179
     180    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux           !< filter function for u in x direction
    162181    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy           !< filter function for u in y direction
    163182    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz           !< filter function for u in z direction
     183    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx           !< filter function for v in x direction
    164184    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy           !< filter function for v in y direction
    165185    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz           !< filter function for v in z direction
     186    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx           !< filter function for w in y direction
    166187    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy           !< filter function for w in y direction
    167188    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz           !< filter function for w in z direction
    168     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu            !< velocity seed for u
    169     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo           !< velocity seed for u with new random number
    170     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv            !< velocity seed for v
    171     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo           !< velocity seed for v with new random number
    172     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw            !< velocity seed for w
    173     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo           !< velocity seed for w with new random number
     189    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz         !< velocity seed for u at xz plane
     190    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz        !< velocity seed for u at xz plane with new random number
     191    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz         !< velocity seed for u at yz plane
     192    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz        !< velocity seed for u at yz plane with new random number
     193    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz         !< velocity seed for v at xz plane
     194    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz        !< velocity seed for v at xz plane with new random number
     195    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz         !< velocity seed for v at yz plane
     196    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz        !< velocity seed for v at yz plane with new random number
     197    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz         !< velocity seed for w at xz plane
     198    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz        !< velocity seed for w at xz plane with new random number
     199    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz         !< velocity seed for w at yz plane
     200    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz        !< velocity seed for w at yz plane with new random number
    174201
    175202
    … …  
    188215
    189216!
    190 !-- Generate velocity seeds
     217!-- Generate velocity seeds at south and north domain boundary
     218    INTERFACE stg_generate_seed_xz
     219       MODULE PROCEDURE stg_generate_seed_xz
     220    END INTERFACE stg_generate_seed_xz
     221!
     222!-- Generate velocity seeds at left and/or right domain boundary
    191223    INTERFACE stg_generate_seed_yz
    192224       MODULE PROCEDURE stg_generate_seed_yz
    … …  
    240272!
    241273!-- Public variables
    242     PUBLIC  use_syn_turb_gen
     274    PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
     275            use_syn_turb_gen
    243276
    244277
    … …  
    255288
    256289    USE control_parameters,                                                    &
    257         ONLY:  bc_lr, bc_ns, turbulent_inflow
     290        ONLY:  bc_lr, bc_ns, forcing, nest_domain, rans_mode, turbulent_inflow
     291
     292    USE pmc_interface,                                                         &
     293        ONLY : rans_mode_parent
    258294
    259295
    260296    IMPLICIT NONE
    261297
     298    IF ( .NOT. use_syn_turb_gen  .AND.  .NOT. rans_mode  .AND.  forcing )  THEN
     299       message_string = 'Synthetic turbulence generator has to be applied ' // &
     300                        'when forcing is used and model operates in LES mode.'
     301       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     302    ENDIF
     303
     304    IF ( .NOT. use_syn_turb_gen  .AND.  nest_domain                            &
     305         .AND. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
     306       message_string = 'Synthetic turbulence generator has to be applied ' // &
     307                        'when nesting is applied and parent operates in '  //  &
     308                        'RANS-mode but current child in LES mode.'
     309       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
     310    ENDIF
     311
    262312    IF ( use_syn_turb_gen )  THEN
    263313
    264        IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0  .AND. &
    265             INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
    266           message_string = 'Using synthetic turbulence generator ' //          &
    267                            'requires &initializing_actions = '            //   &
    268                            '"set_constant_profiles" or "read_restart_data"'
    269           CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
     314       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN
     315
     316          IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0     &
     317        .AND.  INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
     318             message_string = 'Using synthetic turbulence generator ' //       &
     319                              'requires &initializing_actions = '         //   &
     320                              '"set_constant_profiles" or "read_restart_data"'
     321             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
     322          ENDIF
     323
     324          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
     325             message_string = 'Using synthetic turbulence generator ' //       &
     326                              'requires &bc_lr = "dirichlet/radiation"'
     327             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
     328          ENDIF
     329          IF ( bc_ns /= 'cyclic' )  THEN
     330             message_string = 'Using synthetic turbulence generator ' //       &
     331                              'requires &bc_ns = "cyclic"'
     332             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
     333          ENDIF
     334
    270335       ENDIF
    271336
    272        IF ( bc_lr /= 'dirichlet/radiation' )  THEN
    273           message_string = 'Using synthetic turbulence generator ' //          &
    274                            'requires &bc_lr = "dirichlet/radiation"'
    275           CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
    276        ENDIF
    277        IF ( bc_ns /= 'cyclic' )  THEN
    278           message_string = 'Using synthetic turbulence generator ' //          &
    279                            'requires &bc_ns = "cyclic"'
    280           CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
    281        ENDIF
    282337       IF ( turbulent_inflow )  THEN
    283338          message_string = 'Using synthetic turbulence generator ' //          &
    284                            'in combination &with turbulent_inflow = .T. ' //   &
    285                            'is not allowed'
     339                           'in combination &with turbulent_inflow = .T. '//    &
     340                              'is not allowed'
    286341          CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
    287342       ENDIF
    … …  
    330385
    331386    USE arrays_3d,                                                             &
    332         ONLY:  u_init, v_init
     387        ONLY:  ddzw, u_init, v_init, zu
    333388
    334389    USE control_parameters,                                                    &
    335         ONLY:  coupling_char, dz, e_init
     390        ONLY:  coupling_char, dz, e_init, forcing, nest_domain
    336391
    337392    USE grid_variables,                                                        &
    338         ONLY:  ddy
     393        ONLY:  ddx, ddy
     394
     395    USE indices,                                                               &
     396        ONLY:  nz
    339397
    340398
    341399    IMPLICIT NONE
     400
     401    LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
    342402
    343403#if defined( __parallel )
    … …  
    346406#endif
    347407
     408    INTEGER(iwp) :: i                        !> grid index in x-direction
    348409    INTEGER(iwp) :: j                        !> loop index
    349410    INTEGER(iwp) :: k                        !< index
    350411    INTEGER(iwp) :: newtype                  !< dummy MPI type
    351412    INTEGER(iwp) :: realsize                 !< size of REAL variables
     413    INTEGER(iwp) :: nnz                      !< increment used to determine processor decomposition of z-axis along x and y direction 
    352414    INTEGER(iwp) :: nseed                    !< dimension of random number seed
    353415    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
    … …  
    362424    REAL(wp) :: d21     !< vertical interpolation for a21
    363425    REAL(wp) :: d22     !< vertical interpolation for a22
     426    REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
    364427    REAL(wp) :: luy     !< length scale for u in y direction
    365428    REAL(wp) :: luz     !< length scale for u in z direction
    … …  
    370433    REAL(wp) :: zz      !< height
    371434
     435    REAL(wp) :: length_scale_surface, r_ii_0, time_scale, length_scale_z
     436
    372437    REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
    373438    REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
    … …  
    389454    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
    390455               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
    391                nuy(nzb:nzt+1), nuz(nzb:nzt+1), nvy(nzb:nzt+1),                 &
    392                nvz(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
     456               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
     457               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
     458               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
    393459               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
    394460
    395461#if defined( __parallel )
     462!
     463!-- Determine processor decomposition of z-axis along x- and y-direction
     464    nnz = nz / pdims(1)
     465    nzb_x_stg = 1 + myidx * nnz
     466    nzt_x_stg = ( myidx + 1 ) * nnz
     467
     468    IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             & 
     469       nzt_x_stg = myidx * nnz + MOD( nz , pdims(1) )
     470
     471    IF ( forcing  .OR.  nest_domain )  THEN
     472       nnz = nz / pdims(2)
     473       nzb_y_stg = 1 + myidy * nnz
     474       nzt_y_stg = ( myidy + 1 ) * nnz
     475
     476       IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          & 
     477          nzt_y_stg = myidy * nnz + MOD( nz , pdims(2) )
     478    ENDIF
     479
    396480!
    397481!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
    … …  
    399483    CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
    400484    extent = 1 * realsize
    401 
    402     ! stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
     485!
     486!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
     487!-- layer
     488!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
    403489    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
    404490            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    … …  
    407493    CALL MPI_TYPE_FREE( newtype, ierr )
    408494
    409     ! stg_type_yz_small: yz-slice with vertical bounds nzb_x:nzt_x+1
    410     CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x-nzb_x+2,nyng-nysg+1],             &
     495    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
     496    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],             &
    411497            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    412498    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
    … …  
    415501
    416502    ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
    417     ALLOCATE( recv_count(pdims(1)), displs(pdims(1)) )
    418 
    419     recv_count           = nzt_x-nzb_x + 1
    420     recv_count(pdims(1)) = recv_count(pdims(1)) + 1
     503    ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
     504
     505    recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
     506    recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
    421507
    422508    DO  j = 1, pdims(1)
    423        displs(j) = 0 + (nzt_x-nzb_x+1) * (j-1)
    424     ENDDO
     509       displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
     510    ENDDO
     511!
     512!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
     513!-- layer
     514!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
     515    IF ( forcing  .OR.  nest_domain)  THEN
     516       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],                 &
     517               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     518       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
     519       CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
     520       CALL MPI_TYPE_FREE( newtype, ierr )
     521
     522       ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
     523       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],             &
     524               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     525       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
     526       CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
     527       CALL MPI_TYPE_FREE( newtype, ierr )
     528
     529       ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
     530       ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
     531
     532       recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
     533       recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
     534
     535       DO  j = 1, pdims(2)
     536          displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
     537       ENDDO
     538
     539    ENDIF
     540
    425541#endif
    426 
    427542!
    428543!-- Define seed of random number
    … …  
    435550    CALL RANDOM_SEED(put=seed)
    436551
    437 !
    438552!-- Read inflow profile
    439553!-- Height levels of profiles in input profile is as follows:
    … …  
    441555!-- zw: lwy, lwz, tw, r31, r32, r33, d3
    442556!-- WARNING: zz is not used at the moment
    443     OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',        &
    444                    FORM='FORMATTED')
    445 
    446     ! Skip header
    447     READ( 90, * )
    448 
    449     DO  k = nzb, nzt+1
    450        READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k),    &
    451                      r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),           &
    452                      d1, d2, d3, d5
    453 
    454 !
    455 !--    Convert length scales from meter to number of grid points
    456        nuy(k) = INT( luy * ddy )
    457        nuz(k) = INT( luz / dz  )
    458        nvy(k) = INT( lvy * ddy )
    459        nvz(k) = INT( lvz / dz  )
    460        nwy(k) = INT( lwy * ddy )
    461        nwz(k) = INT( lwz / dz  )
    462 
    463 !
    464 !--    Save Mean inflow profiles
    465        IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
    466           mean_inflow_profiles(k,1) = d1
    467           mean_inflow_profiles(k,2) = d2
    468          !  mean_inflow_profiles(k,4) = d4
    469           mean_inflow_profiles(k,5) = d5
    470        ENDIF
    471     ENDDO
    472 
    473     CLOSE( 90 )
    474 
    475 !
    476 !--    Assign initial profiles
    477     u_init = mean_inflow_profiles(:,1)
    478     v_init = mean_inflow_profiles(:,2)
    479    !  pt_init = mean_inflow_profiles(:,4)
    480     e_init = MAXVAL( mean_inflow_profiles(:,5) )
    481 
     557    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ),                   &
     558             EXIST = file_stg_exist  )
     559
     560    IF ( file_stg_exist )  THEN
     561
     562       OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',     &
     563                      FORM='FORMATTED')
     564!
     565!--    Skip header
     566       READ( 90, * )
     567
     568       DO  k = nzb, nzt+1
     569          READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
     570                        r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),        &
     571                        d1, d2, d3, d5
     572
     573!
     574!--       Convert length scales from meter to number of grid points
     575          nuy(k) = INT( luy * ddy )
     576          nuz(k) = INT( luz / dz  )
     577          nvy(k) = INT( lvy * ddy )
     578          nvz(k) = INT( lvz / dz  )
     579          nwy(k) = INT( lwy * ddy )
     580          nwz(k) = INT( lwz / dz  )
     581!
     582!--       Workaround, assume isotropic turbulence
     583          nwx(k) = nwy(k)
     584          nvx(k) = nvy(k)
     585          nux(k) = nuy(k)
     586!
     587!--       Save Mean inflow profiles
     588          IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
     589             mean_inflow_profiles(k,1) = d1
     590             mean_inflow_profiles(k,2) = d2
     591            !  mean_inflow_profiles(k,4) = d4
     592             mean_inflow_profiles(k,5) = d5
     593          ENDIF
     594       ENDDO
     595
     596       CLOSE( 90 )
     597   
     598    ELSE
     599!
     600!--    Set-up defaul length scales. Assume exponentially decreasing length
     601!--    scales and isotropic turbulence.
     602!--    Typical length (time) scales of 100 m (s) should be a goog compromise
     603!--    between all stratrifications. Near-surface variances are fixed to
     604!--    0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
     605!--    Vertical fluxes 
     606       length_scale_surface = 100.0_wp
     607       r_ii_0               = 0.1_wp
     608       time_scale           = 100.0_wp
     609       DO  k = nzb+1, nzt+1
     610          dum_exp        = MERGE( -zu(k) / ( 0.3* zu(nzt) ),                   &
     611                                  0.0_wp,                                      &
     612                                  zu(k) > 0.3 * zu(nzt)                        &
     613                                )
     614          length_scale_z = length_scale_surface * EXP( dum_exp )
     615
     616          nux(k) = MAX( INT( length_scale_z * ddx     ), 1 )
     617          nuy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
     618          nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
     619          nvx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
     620          nvy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
     621          nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
     622          nwx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
     623          nwy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
     624          nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
     625
     626          r11(k) = r_ii_0 * EXP( dum_exp )
     627          r22(k) = r_ii_0 * EXP( dum_exp )
     628          r33(k) = r_ii_0 * EXP( dum_exp )
     629
     630          r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
     631          r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
     632          r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
     633
     634          tu(k)  = time_scale
     635          tv(k)  = time_scale
     636          tw(k)  = time_scale
     637
     638       ENDDO
     639       nux(nzb) = nux(nzb+1)
     640       nuy(nzb) = nuy(nzb+1)
     641       nuz(nzb) = nuz(nzb+1)
     642       nvx(nzb) = nvx(nzb+1)
     643       nvy(nzb) = nvy(nzb+1)
     644       nvz(nzb) = nvz(nzb+1)
     645       nwx(nzb) = nwx(nzb+1)
     646       nwy(nzb) = nwy(nzb+1)
     647       nwz(nzb) = nwz(nzb+1)
     648
     649       r11(nzb) = r11(nzb+1)
     650       r22(nzb) = r22(nzb+1)
     651       r33(nzb) = r33(nzb+1)
     652
     653       r21(nzb) = r11(nzb+1)
     654       r31(nzb) = r31(nzb+1)
     655       r32(nzb) = r32(nzb+1)
     656
     657       tu(nzb)  = time_scale
     658       tv(nzb)  = time_scale
     659       tw(nzb)  = time_scale
     660
     661    ENDIF
     662
     663!
     664!-- Assign initial profiles
     665    IF ( .NOT. forcing  .AND.  .NOT.  nest_domain )  THEN
     666       u_init = mean_inflow_profiles(:,1)
     667       v_init = mean_inflow_profiles(:,2)
     668      !pt_init = mean_inflow_profiles(:,4)
     669       e_init = MAXVAL( mean_inflow_profiles(:,5) )
     670    ENDIF
    482671!
    483672!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
    … …  
    534723
    535724    ENDDO
    536 
    537725!
    538726!-- Define the size of the filter functions and allocate them.
    … …  
    541729    ! arrays must be large enough to cover the largest length scale
    542730    DO  k = nzb, nzt+1
    543        j = MAX( ABS(nuy(k)), ABS(nuz(k)), &
    544                 ABS(nvy(k)), ABS(nvz(k)), &
    545                 ABS(nwy(k)), ABS(nwz(k))  )
     731       j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
     732                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
     733                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
    546734       IF ( j > merg )  merg = j
    547735    ENDDO
    … …  
    550738    mergp = merg + nbgp
    551739
    552     ALLOCATE ( buy(-merg:merg,nzb:nzt+1), buz(-merg:merg,nzb:nzt+1), &
    553                bvy(-merg:merg,nzb:nzt+1), bvz(-merg:merg,nzb:nzt+1), &
    554                bwy(-merg:merg,nzb:nzt+1), bwz(-merg:merg,nzb:nzt+1)  )
    555 
    556 !
    557 !-- Allocate velocity seeds
    558     ALLOCATE ( fu( nzb:nzt+1,nysg:nyng), fuo(nzb:nzt+1,nysg:nyng), &
    559                fv( nzb:nzt+1,nysg:nyng), fvo(nzb:nzt+1,nysg:nyng), &
    560                fw( nzb:nzt+1,nysg:nyng), fwo(nzb:nzt+1,nysg:nyng)  )
    561 
    562     fu  = 0._wp
    563     fuo = 0._wp
    564     fv  = 0._wp
    565     fvo = 0._wp
    566     fw  = 0._wp
    567     fwo = 0._wp
     740    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
     741               buy(-merg:merg,nzb:nzt+1),                                      &
     742               buz(-merg:merg,nzb:nzt+1),                                      &
     743               bvx(-merg:merg,nzb:nzt+1),                                      &
     744               bvy(-merg:merg,nzb:nzt+1),                                      &
     745               bvz(-merg:merg,nzb:nzt+1),                                      &
     746               bwx(-merg:merg,nzb:nzt+1),                                      &
     747               bwy(-merg:merg,nzb:nzt+1),                                      &
     748               bwz(-merg:merg,nzb:nzt+1)  )
     749
     750!
     751!-- Allocate velocity seeds for turbulence at xz-layer
     752    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
     753               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
     754               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
     755
     756!
     757!-- Allocate velocity seeds for turbulence at yz-layer
     758    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
     759               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
     760               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
     761
     762    fu_xz  = 0.0_wp
     763    fuo_xz = 0.0_wp
     764    fv_xz  = 0.0_wp
     765    fvo_xz = 0.0_wp
     766    fw_xz  = 0.0_wp
     767    fwo_xz = 0.0_wp
     768
     769    fu_yz  = 0.0_wp
     770    fuo_yz = 0.0_wp
     771    fv_yz  = 0.0_wp
     772    fvo_yz = 0.0_wp
     773    fw_yz  = 0.0_wp
     774    fwo_yz = 0.0_wp
    568775
    569776!
    570777!-- Create filter functions
     778    CALL stg_filter_func( nux, bux ) !filter ux
    571779    CALL stg_filter_func( nuy, buy ) !filter uy
    572780    CALL stg_filter_func( nuz, buz ) !filter uz
     781    CALL stg_filter_func( nvx, bvx ) !filter vx
    573782    CALL stg_filter_func( nvy, bvy ) !filter vy
    574783    CALL stg_filter_func( nvz, bvz ) !filter vz
     784    CALL stg_filter_func( nwx, bwx ) !filter wx
    575785    CALL stg_filter_func( nwy, bwy ) !filter wy
    576786    CALL stg_filter_func( nwz, bwz ) !filter wz
    … …  
    581791
    582792!
    583 !--    In case of restart, calculate velocity seeds fu, fv, fw from former
    584 !      time step.
    585 !      Bug: fu, fv, fw are different in those heights where a11, a22, a33
    586 !           are 0 compared to the prerun. This is mostly for k=nzt+1.
     793!-- In case of restart, calculate velocity seeds fu, fv, fw from former
     794!   time step.
     795!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
     796!        are 0 compared to the prerun. This is mostly for k=nzt+1.
    587797    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
    588        IF ( myidx == id_inflow )  THEN
     798       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
     799
     800          IF ( myidx == id_stg_left  )  i = -1
     801          IF ( myidx == id_stg_right )  i = nxr+1
     802
    589803          DO  j = nysg, nyng
    590804             DO  k = nzb, nzt+1
    591805
    592806                IF  ( a11(k) .NE. 0._wp ) THEN
    593                 fu(k,j) = ( u(k,j,-1) / mc_factor -                    &
    594                             mean_inflow_profiles(k,1) ) / a11(k)
     807                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
    595808                ELSE
    596                    fu(k,j) = 0._wp
     809                   fu_yz(k,j) = 0._wp
    597810                ENDIF
    598811
    599812                IF  ( a22(k) .NE. 0._wp ) THEN
    600                 fv(k,j) = ( v(k,j,-1) / mc_factor - a21(k) * fu(k,j) - &
    601                             mean_inflow_profiles(k,2) ) / a22(k)
     813                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
     814                               v_init(k) ) / a22(k)
    602815                ELSE
    603                    fv(k,j) = 0._wp
     816                   fv_yz(k,j) = 0._wp
    604817                ENDIF
    605818
    606819                IF  ( a33(k) .NE. 0._wp ) THEN
    607                 fw(k,j) = ( w(k,j,-1) / mc_factor - a31(k) * fu(k,j) - &
    608                             a32(k) * fv(k,j) ) / a33(k)
     820                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
     821                               a32(k) * fv_yz(k,j) ) / a33(k)
    609822                ELSE
    610                    fw = 0._wp
     823                   fw_yz = 0._wp
    611824                ENDIF
    612825
    … …  
    759972    IMPLICIT NONE
    760973
    761 
    762974    CALL wrd_write_string( 'mc_factor' )
    763975    WRITE ( 14 )  mc_factor
    … …  
    785997
    786998    USE control_parameters,                                                    &
    787         ONLY:  dt_3d, intermediate_timestep_count, simulated_time,             &
    788                volume_flow_initial
     999        ONLY:  dt_3d, forcing, intermediate_timestep_count,  nest_domain,      &
     1000               simulated_time, volume_flow_initial
     1001
     1002    USE grid_variables,                                                        &
     1003        ONLY:  dx, dy
    7891004
    7901005    USE indices,                                                               &
    791         ONLY:  nyn, nys
     1006        ONLY:  wall_flags_0
    7921007
    7931008    USE statistics,                                                            &
    … …  
    7971012    IMPLICIT NONE
    7981013
     1014    INTEGER(iwp) :: i           !< grid index in x-direction
    7991015    INTEGER(iwp) :: j           !< loop index in y-direction
    8001016    INTEGER(iwp) :: k           !< loop index in z-direction
    … …  
    8021018    REAL(wp) :: dt_stg          !< wheighted subtimestep
    8031019    REAL(wp) :: mc_factor_l     !< local mass flux correction factor
    804 
    805     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) :: inflow_dist !< disturbance for inflow profiles
     1020    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
     1021    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
     1022
     1023    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
     1024    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
    8061025
    8071026
    … …  
    8151034!-- Initial value of fu, fv, fw
    8161035    IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
    817        CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu )
    818        CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv )
    819        CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw )
     1036       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
     1037       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
     1038       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
     1039
     1040       IF ( forcing  .OR.  nest_domain )  THEN
     1041!
     1042!--       Generate turbulence at right boundary
     1043          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
     1044          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
     1045          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
     1046!
     1047!--       Generate turbulence at north boundary
     1048          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
     1049          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
     1050          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
     1051!
     1052!--       Generate turbulence at south boundary
     1053          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
     1054          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
     1055          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
     1056       ENDIF
    8201057       velocity_seed_initialized = .TRUE.
    8211058    ENDIF
    8221059!
    8231060!-- New set of fu, fv, fw
    824     CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo )
    825     CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo )
    826     CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo )
    827 
    828     IF ( myidx == id_inflow )  THEN
     1061    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
     1062    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
     1063    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
     1064
     1065    IF ( forcing  .OR.  nest_domain )  THEN
     1066!
     1067!--       Generate turbulence at right boundary
     1068          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
     1069          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
     1070          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
     1071!
     1072!--       Generate turbulence at north boundary
     1073          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
     1074          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
     1075          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
     1076!
     1077!--       Generate turbulence at south boundary
     1078          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
     1079          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
     1080          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1081    ENDIF
     1082!
     1083!-- Turbulence generation at left and or right boundary
     1084    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
    8291085
    8301086       DO  j = nysg, nyng
    … …  
    8331089!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
    8341090             IF ( tu(k) == 0.0_wp )  THEN
    835                 fu(k,j) = fuo(k,j)
     1091                fu_yz(k,j) = fuo_yz(k,j)
    8361092             ELSE
    837                 fu(k,j) = fu(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
    838                          fuo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
     1093                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
     1094                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
    8391095             ENDIF
    8401096
    8411097             IF ( tv(k) == 0.0_wp )  THEN
    842                 fv(k,j) = fvo(k,j)
     1098                fv_yz(k,j) = fvo_yz(k,j)
    8431099             ELSE
    844                 fv(k,j) = fv(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
    845                          fvo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1100                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
     1101                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    8461102             ENDIF
    8471103
    8481104             IF ( tw(k) == 0.0_wp )  THEN
    849                 fw(k,j) = fwo(k,j)
     1105                fw_yz(k,j) = fwo_yz(k,j)
    8501106             ELSE
    851                 fw(k,j) = fw(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
    852                          fwo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1107                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
     1108                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    8531109             ENDIF
    8541110!
    … …  
    8561112!--          Additional factors are added to improve the variance of v and w
    8571113             IF( k == 0 )  THEN
    858                 inflow_dist(k,j,1,1) = 0.0_wp
    859                 inflow_dist(k,j,2,1) = 0.0_wp
    860                 inflow_dist(k,j,3,1) = 0.0_wp
    861                 inflow_dist(k,j,4,1) = 0.0_wp
    862                 inflow_dist(k,j,5,1) = 0.0_wp
     1114                dist_yz(k,j,1) = 0.0_wp
     1115                dist_yz(k,j,2) = 0.0_wp
     1116                dist_yz(k,j,3) = 0.0_wp
     1117!                 dist_yz(k,j,4) = 0.0_wp
     1118!                 dist_yz(k,j,5) = 0.0_wp
    8631119             ELSE
    864                 inflow_dist(k,j,1,1) = a11(k) * fu(k,j)
     1120                dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
    8651121                !experimental test of 1.2
    866                 inflow_dist(k,j,2,1) = ( SQRT( a22(k) / MAXVAL(a22) )          &
     1122                dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                   &
    8671123                                         * 1.2_wp )                            &
    868                                        * (   a21(k) * fu(k,j)                  &
    869                                            + a22(k) * fv(k,j) )
    870                 inflow_dist(k,j,3,1) = ( SQRT(a33(k) / MAXVAL(a33) )           &
     1124                                       * (   a21(k) * fu_yz(k,j)                  &
     1125                                           + a22(k) * fv_yz(k,j) )
     1126                dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                    &
    8711127                                         * 1.3_wp )                            &
    872                                        * (   a31(k) * fu(k,j)                  &
    873                                            + a32(k) * fv(k,j)                  &
    874                                            + a33(k) * fw(k,j) )
     1128                                       * (   a31(k) * fu_yz(k,j)                  &
     1129                                           + a32(k) * fv_yz(k,j)                  &
     1130                                           + a33(k) * fw_yz(k,j) )
    8751131                ! Calculation for pt and e not yet implemented
    876                 inflow_dist(k,j,4,1) = 0.0_wp
    877                 inflow_dist(k,j,5,1) = 0.0_wp
     1132!                 dist_yz(k,j,4) = 0.0_wp
     1133!                 dist_yz(k,j,5) = 0.0_wp
    8781134             ENDIF
    8791135
    … …  
    8851141!--    This correction factor insures that the mass flux is preserved at the
    8861142!--    inflow boundary
    887        mc_factor_l = 0.0_wp
    888        mc_factor   = 0.0_wp
    889        DO  j = nys, nyn
     1143       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN
     1144          mc_factor_l = 0.0_wp
     1145          mc_factor   = 0.0_wp
     1146          DO  j = nys, nyn
     1147             DO  k = nzb+1, nzt
     1148                mc_factor_l = mc_factor_l + dzw(k)  *                          &
     1149                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
     1150             ENDDO
     1151          ENDDO
     1152
     1153#if defined( __parallel )
     1154          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,  &
     1155                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1156#else
     1157          mc_factor = mc_factor_l
     1158#endif
     1159
     1160          mc_factor = volume_flow_initial(1) / mc_factor
     1161
     1162!
     1163!--       Add disturbance at the inflow
     1164          DO  j = nysg, nyng
     1165             DO  k = nzb, nzt+1
     1166                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
     1167                                      dist_yz(k,j,1)             ) * mc_factor
     1168                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
     1169                                      dist_yz(k,j,2)             ) * mc_factor
     1170                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
     1171             ENDDO
     1172          ENDDO
     1173
     1174       ELSE
     1175!
     1176!--       First, calculate volume flow at yz boundary
     1177          IF ( myidx == id_stg_left  )  i = nxl
     1178          IF ( myidx == id_stg_right )  i = nxr+1
     1179
     1180          volume_flow_l = 0.0_wp
     1181          volume_flow   = 0.0_wp
     1182          mc_factor_l   = 0.0_wp
     1183          mc_factor     = 0.0_wp
     1184          DO  j = nys, nyn
     1185             DO  k = nzb+1, nzt
     1186                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
     1187                                     * MERGE( 1.0_wp, 0.0_wp,                  &   
     1188                                              BTEST( wall_flags_0(k,j,i), 1 ) )
     1189
     1190                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
     1191                                                         * dzw(k) * dy         &
     1192                                     * MERGE( 1.0_wp, 0.0_wp,                  &   
     1193                                              BTEST( wall_flags_0(k,j,i), 1 ) )
     1194             ENDDO
     1195          ENDDO
     1196#if defined( __parallel )
     1197          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
     1198                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1199          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
     1200                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1201#else
     1202          volume_flow = volume_flow_l
     1203          mc_factor   = mc_factor_l
     1204#endif
     1205
     1206          mc_factor = volume_flow / mc_factor
     1207
     1208!
     1209!--       Add disturbances
     1210          IF ( myidx == id_stg_left  )  THEN
     1211
     1212             DO  j = nysg, nyng
     1213                DO  k = nzb, nzt+1
     1214                   u(k,j,-nbgp+1:0) = ( u(k,j,-nbgp+1:0) + dist_yz(k,j,1) )    &
     1215                                        * mc_factor
     1216                   v(k,j,-nbgp:-1)  = ( v(k,j,-nbgp:-1)  + dist_yz(k,j,2)  )   &
     1217                                        * mc_factor
     1218                   w(k,j,-nbgp:-1)  = ( w(k,j,-nbgp:-1)  + dist_yz(k,j,3)  )   &
     1219                                        * mc_factor
     1220                ENDDO
     1221             ENDDO
     1222          ENDIF
     1223          IF ( myidx == id_stg_right  )  THEN
     1224
     1225             DO  j = nysg, nyng
     1226                DO  k = nzb, nzt+1
     1227                   u(k,j,nxr+1:nxr+nbgp) = ( u(k,j,nxr+1:nxr+nbgp) +           &
     1228                                             dist_yz(k,j,1) ) * mc_factor
     1229                   v(k,j,nxr+1:nxr+nbgp) = ( v(k,j,nxr+1:nxr+nbgp) +           &
     1230                                             dist_yz(k,j,2) ) * mc_factor
     1231                   w(k,j,nxr+1:nxr+nbgp) = ( w(k,j,nxr+1:nxr+nbgp) +           &
     1232                                             dist_yz(k,j,3) ) * mc_factor
     1233                ENDDO
     1234             ENDDO
     1235          ENDIF
     1236       ENDIF
     1237
     1238    ENDIF
     1239!
     1240!-- Turbulence generation at north and south boundary
     1241    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
     1242
     1243       DO  i = nxlg, nxrg
     1244          DO  k = nzb, nzt + 1
     1245!
     1246!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
     1247             IF ( tu(k) == 0.0_wp )  THEN
     1248                fu_xz(k,i) = fuo_xz(k,i)
     1249             ELSE
     1250                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
     1251                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
     1252             ENDIF
     1253
     1254             IF ( tv(k) == 0.0_wp )  THEN
     1255                fv_xz(k,i) = fvo_xz(k,i)
     1256             ELSE
     1257                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
     1258                         fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1259             ENDIF
     1260
     1261             IF ( tw(k) == 0.0_wp )  THEN
     1262                fw_xz(k,i) = fwo_xz(k,i)
     1263             ELSE
     1264                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
     1265                         fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1266             ENDIF
     1267!
     1268!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
     1269!--          Additional factors are added to improve the variance of v and w
     1270             IF( k == 0 )  THEN
     1271                dist_xz(k,i,1) = 0.0_wp
     1272                dist_xz(k,i,2) = 0.0_wp
     1273                dist_xz(k,i,3) = 0.0_wp
     1274
     1275             ELSE
     1276                dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
     1277                !experimental test of 1.2
     1278                dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
     1279                                         * 1.2_wp )                            &
     1280                                       * (   a21(k) * fu_xz(k,i)               &
     1281                                           + a22(k) * fv_xz(k,i) )
     1282                dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
     1283                                         * 1.3_wp )                            &
     1284                                       * (   a31(k) * fu_xz(k,i)               &
     1285                                           + a32(k) * fv_xz(k,i)               &
     1286                                           + a33(k) * fw_xz(k,i) )
     1287             ENDIF
     1288
     1289          ENDDO
     1290       ENDDO
     1291!
     1292!--    Mass flux correction following Kim et al. (2013)
     1293!--    This correction factor insures that the mass flux is preserved at the
     1294!--    inflow boundary.
     1295!--    First, calculate volume flow at xz boundary
     1296       IF ( myidy == id_stg_south  ) j = nys
     1297       IF ( myidy == id_stg_north )  j = nyn+1
     1298
     1299       volume_flow_l = 0.0_wp
     1300       volume_flow   = 0.0_wp
     1301       mc_factor_l   = 0.0_wp
     1302       mc_factor     = 0.0_wp
     1303       DO  i = nxl, nxr
    8901304          DO  k = nzb+1, nzt
    891              mc_factor_l = mc_factor_l + dzw(k)  *                             &
    892                            ( mean_inflow_profiles(k,1) + inflow_dist(k,j,1,1) )
     1305             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx         &
     1306                                  * MERGE( 1.0_wp, 0.0_wp,                  &   
     1307                                           BTEST( wall_flags_0(k,j,i), 2 ) )
     1308
     1309             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )  &
     1310                                                      * dzw(k) * dx         &
     1311                                  * MERGE( 1.0_wp, 0.0_wp,                  &   
     1312                                           BTEST( wall_flags_0(k,j,i), 2 ) )
    8931313          ENDDO
    8941314       ENDDO
    895 
    8961315#if defined( __parallel )
    897        CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,  &
    898                            1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1316       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
     1317                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
     1318       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
     1319                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
    8991320#else
    900        mc_factor = mc_factor_l
     1321       volume_flow = volume_flow_l
     1322       mc_factor   = mc_factor_l
    9011323#endif
    9021324
    903        mc_factor = volume_flow_initial(1) / mc_factor
    904 
    905 !
    906 !--    Add disturbance at the inflow
    907        DO  j = nysg, nyng
    908           DO  k = nzb, nzt+1
    909               u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +                 &
    910                                    inflow_dist(k,j,1,1)      ) * mc_factor
    911               v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +                 &
    912                                    inflow_dist(k,j,2,1)      ) * mc_factor
    913               w(k,j,-nbgp:-1)  =   inflow_dist(k,j,3,1)        * mc_factor
     1325       mc_factor = volume_flow / mc_factor
     1326
     1327!
     1328!--    Add disturbances
     1329       IF ( myidy == id_stg_south  )  THEN
     1330
     1331          DO  i = nxlg, nxrg
     1332             DO  k = nzb, nzt+1
     1333                u(k,-nbgp:-1,i) = ( u(k,-nbgp:-1,i) + dist_xz(k,i,1) )         &
     1334                                     * mc_factor
     1335                v(k,-nbgp:0,i)  = ( v(k,-nbgp:0,i)  + dist_xz(k,i,2)  )        &
     1336                                     * mc_factor
     1337                w(k,-nbgp:-1,i) = ( w(k,-nbgp:-1,i) + dist_xz(k,i,3)  )        &
     1338                                     * mc_factor
     1339             ENDDO
    9141340          ENDDO
    915        ENDDO
    916 
     1341       ENDIF
     1342       IF ( myidy == id_stg_north  )  THEN
     1343
     1344          DO  i = nxlg, nxrg
     1345             DO  k = nzb, nzt+1
     1346                u(k,nyn+1:nyn+nbgp,i) = ( u(k,nyn+1:nyn+nbgp,i) +              &
     1347                                          dist_xz(k,i,1) ) * mc_factor
     1348                v(k,nyn+1:nyn+nbgp,i) = ( v(k,nyn+1:nyn+nbgp,i) +              &
     1349                                          dist_xz(k,i,2) ) * mc_factor
     1350                w(k,nyn+1:nyn+nbgp,i) = ( w(k,nyn+1:nyn+nbgp,i) +              &
     1351                                          dist_xz(k,i,3) ) * mc_factor
     1352             ENDDO
     1353          ENDDO
     1354       ENDIF
    9171355    ENDIF
    9181356
    … …  
    9201358
    9211359 END SUBROUTINE stg_main
    922 
    9231360
    9241361!------------------------------------------------------------------------------!
    … …  
    9311368!> parts are collected to form the full array.
    9321369!------------------------------------------------------------------------------!
    933  SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n )
     1370 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
    9341371
    9351372
    … …  
    9401377    IMPLICIT NONE
    9411378
     1379    INTEGER(iwp) :: id          !< core ids at respective boundaries
    9421380    INTEGER(iwp) :: j           !< loop index in y-direction
    9431381    INTEGER(iwp) :: jj          !< loop index in y-direction
    … …  
    9571395    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
    9581396    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
    959     REAL(wp), DIMENSION(nzb_x:nzt_x+1,nysg:nyng) :: f_n_l   !<  local velocity seed
     1397    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
    9601398    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
    9611399    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
    … …  
    10181456
    10191457    DO  j = nysg, nyng
    1020        DO  k = nzb_x, nzt_x+1
     1458       DO  k = nzb_x_stg, nzt_x_stg+1
    10211459          DO  jj = -n_y2(k), n_y2(k)
    10221460             DO  kk = -n_z2(k), n_z2(k)
    … …  
    10321470!
    10331471!-- Gather velocity seeds of full subdomain
    1034     send_count = nzt_x - nzb_x + 1
    1035     IF ( nzt_x == nzt )  send_count = send_count + 1
     1472    send_count = nzt_x_stg - nzb_x_stg + 1
     1473    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
    10361474
    10371475#if defined( __parallel )
    1038     CALL MPI_GATHERV( f_n_l(nzb_x,nysg), send_count, stg_type_yz_small,        &
    1039                       f_n(nzb+1,nysg), recv_count, displs, stg_type_yz,        &
    1040                       id_inflow, comm1dx, ierr )
     1476    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,        &
     1477                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,        &
     1478                      id, comm1dx, ierr )
    10411479#else
    1042     f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x:nzt_x+1,nysg:nyng)
     1480    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
    10431481#endif
    10441482
    … …  
    10461484 END SUBROUTINE stg_generate_seed_yz
    10471485
     1486
     1487
     1488
     1489!------------------------------------------------------------------------------!
     1490! Description:
     1491! ------------
     1492!> Generate a set of random number rand_it wich is equal on each PE
     1493!> and calculate the velocity seed f_n.
     1494!> f_n is splitted in vertical direction by the number of PEs in y-direction and
     1495!> and each PE calculates a vertical subsection of f_n. At the the end, all
     1496!> parts are collected to form the full array.
     1497!------------------------------------------------------------------------------!
     1498 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
     1499
     1500
     1501    USE indices,                                                               &
     1502        ONLY: nx
     1503
     1504
     1505    IMPLICIT NONE
     1506
     1507    INTEGER(iwp) :: id          !< core ids at respective boundaries
     1508    INTEGER(iwp) :: i           !< loop index in x-direction
     1509    INTEGER(iwp) :: ii          !< loop index in x-direction
     1510    INTEGER(iwp) :: k           !< loop index in z-direction
     1511    INTEGER(iwp) :: kk          !< loop index in z-direction
     1512    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
     1513
     1514    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
     1515    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
     1516    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
     1517    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
     1518
     1519    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
     1520    REAL(wp) :: rand_av         !< average of random number
     1521    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
     1522
     1523    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
     1524    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
     1525    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
     1526    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
     1527    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
     1528
     1529
     1530!
     1531!-- Generate random numbers using a seed generated in stg_init.
     1532!-- The set of random numbers are modified to have an average of 0 and
     1533!-- unit variance.
     1534    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
     1535
     1536    rand_av        = 0.0_wp
     1537    rand_sigma_inv = 0.0_wp
     1538    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
     1539
     1540    DO  i = 0, nx
     1541       DO  k = nzb, nzt+1
     1542          CALL RANDOM_NUMBER( rand_it(k,i) )
     1543          rand_av = rand_av + rand_it(k,i)
     1544       ENDDO
     1545    ENDDO
     1546
     1547    rand_av = rand_av * nxz_inv
     1548
     1549    DO  i = 0, nx
     1550       DO  k = nzb, nzt+1
     1551          rand_it(k,i)   = rand_it(k,i) - rand_av
     1552          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
     1553       ENDDO
     1554    ENDDO
     1555
     1556    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
     1557
     1558    DO  i = 0, nx
     1559       DO  k = nzb, nzt+1
     1560          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
     1561       ENDDO
     1562    ENDDO
     1563
     1564!
     1565!-- Periodic fill of random number in space
     1566    DO  i = 0, nx
     1567       DO  k = 1, mergp
     1568          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
     1569          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
     1570       ENDDO
     1571    ENDDO
     1572    DO  i = 1, mergp
     1573       DO  k = nzb-mergp, nzt+1+mergp
     1574          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
     1575          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
     1576       ENDDO
     1577    ENDDO
     1578
     1579!
     1580!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
     1581    n_x2 = n_x * 2
     1582    n_z2 = n_z * 2
     1583    f_n_l  = 0.0_wp
     1584
     1585    DO  i = nxlg, nxrg
     1586       DO  k = nzb_y_stg, nzt_y_stg+1
     1587          DO  ii = -n_x2(k), n_x2(k)
     1588             DO  kk = -n_z2(k), n_z2(k)
     1589                f_n_l(k,i) = f_n_l(k,i)                                        &
     1590                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
     1591             ENDDO
     1592          ENDDO
     1593       ENDDO
     1594    ENDDO
     1595
     1596    DEALLOCATE( rand_it )
     1597
     1598!
     1599!-- Gather velocity seeds of full subdomain
     1600    send_count = nzt_y_stg - nzb_y_stg + 1
     1601    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
     1602
     1603
     1604#if defined( __parallel )
     1605    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,        &
     1606                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
     1607                      id, comm1dy, ierr )
     1608#else
     1609    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
     1610#endif
     1611
     1612
     1613 END SUBROUTINE stg_generate_seed_xz
     1614
    10481615 END MODULE synthetic_turbulence_generator_mod

  • palm/trunk/SOURCE/time_integration.f90

    r2934 r2938  
    2525! -----------------
    2626! $Id$
     27! Nesting of dissipation rate in case of RANS mode and TKE-e closure is applied
     28!
     29! 2936 2018-03-27 14:49:27Z suehring
    2730! Little formatting adjustment.
    2831!
    … …  
    358361               microphysics_morrison, microphysics_seifert, mid, nest_domain,  &
    359362               neutral, nr_timesteps_this_run, nudging,                        &
    360                ocean, passive_scalar,                                          &
    361                prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
     363               ocean, passive_scalar, prho_reference, pt_reference,            &
     364               pt_slope_offset, random_heatflux, rans_mode,                    &
    362365               rans_tke_e, run_coupled, simulated_time, simulated_time_chr,    &
    363366               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
    … …  
    782785                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
    783786
     787                IF ( .NOT. constant_diffusion  .AND.  rans_mode  .AND.         &
     788                                                      rans_tke_e )             &
     789                   CALL exchange_horiz( diss, nbgp )
     790
    784791                IF ( air_chemistry )  THEN
    785792                   DO  n = 1, nspec     
    … …  
    820827!
    821828!--       Impose a turbulent inflow using synthetic generated turbulence
    822           IF ( use_syn_turb_gen ) THEN
    823              CALL  stg_main
    824           ENDIF
     829          IF ( use_syn_turb_gen )  CALL  stg_main
    825830
    826831!
    … …  
    862867                               intermediate_timestep_count_max  )  THEN
    863868             CALL forcing_bc
    864           ENDIF
    865 !
    866 !--       Moreover, ensure mass conservation at each RK substep.
    867           IF ( forcing )  CALL forcing_bc_mass_conservation
     869!
     870!--          Moreover, ensure mass conservation
     871             CALL forcing_bc_mass_conservation
     872          ENDIF
    868873
    869874!

  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r2918 r2938  
    2525! -----------------
    2626! $Id$
     27! Further todo's
     28!
     29! 2936 2018-03-27 14:49:27Z suehring
    2730! - defined l_grid only within this module
    2831! - Moved l_wall definition from modules.f90
    … …  
    6972!>       add OpenMP directives whereever possible
    7073!>       remove debug output variables (dummy1, dummy2, dummy3)
     74!> @todo Move initialization of wall-mixing length from init_grid
     75!> @todo Check for random disturbances
    7176!> @note <Enter notes on the module>
    7277!> @bug  TKE-e closure still crashes due to too small dt
    … …  
    284289
    285290    USE control_parameters,                                                    &
    286         ONLY:  message_string, neutral, turbulent_inflow, turbulent_outflow
     291        ONLY:  message_string, nest_domain, neutral, turbulent_inflow,         &
     292               turbulent_outflow
    287293
    288294    IMPLICIT NONE
    … …  
    302308             rans_tke_e = .TRUE.
    303309
    304              IF ( INDEX( initializing_actions, 'set_1d-model_profiles' )       &
    305                   == 0 )  THEN
     310             IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0  &
     311                  .AND.  .NOT.  nest_domain )  THEN
    306312                message_string = 'Initializing without 1D model while ' //     &
    307313                                 'using TKE-e closure&' //                     &
    … …  
    826832        ONLY:  use_sgs_for_particles, wang_kernel
    827833
     834    USE pmc_interface,                                                         &
     835        ONLY:  nested_run
     836
    828837    IMPLICIT NONE
    829838
    … …  
    847856    ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    848857#endif
    849 
     858!
     859!-- Allocate arrays required for dissipation.
     860!-- Please note, if it is a nested run, arrays need to be allocated even if
     861!-- they do not necessarily need to be transferred, which is attributed to
     862!-- the design of the model coupler which allocates memory for each variable.
    850863    IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
    851          collision_turbulence )  THEN
     864         collision_turbulence  .OR.  nested_run )  THEN
    852865#if defined( __nopointer )
    853866       ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    … …  
    858871#else
    859872       ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    860        IF ( rans_tke_e )  THEN
     873       IF ( rans_tke_e  .OR.  nested_run )  THEN
    861874          ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    862875          ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    … …  
    871884
    872885    IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.     &
    873          wang_kernel  .OR.  collision_turbulence )  THEN
     886         wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
    874887       diss => diss_1
    875        IF ( rans_tke_e )  THEN
     888       IF ( rans_tke_e  .OR.  nested_run )  THEN
    876889       diss_p => diss_2; tdiss_m => diss_3
    877890       ENDIF
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |