source: palm/trunk/SOURCE/pmc_interface_mod.f90 @ 4329

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

  • Property svn:keywords set to Id
File size: 221.6 KB
Line 
1!> @file pmc_interface_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: pmc_interface_mod.f90 4329 2019-12-10 15:46:36Z motisi $
27! Renamed wall_flags_0 to wall_flags_static_0
28!
29! 4273 2019-10-24 13:40:54Z monakurppa
30! Add a logical switch nesting_chem and rename nest_salsa to nesting_salsa
31!
32! 4260 2019-10-09 14:04:03Z hellstea
33! Rest of the possibly round-off-error sensitive grid-line matching tests
34! changed to round-off-error tolerant forms throughout the module.
35!
36! 4249 2019-10-01 12:27:47Z hellstea
37! Several grid-line matching tests changed to a round-off-error tolerant form
38! in pmci_setup_parent, pmci_define_index_mapping and pmci_check_grid_matching.
39!
40! 4182 2019-08-22 15:20:23Z scharf
41! Corrected "Former revisions" section
42!
43! 4168 2019-08-16 13:50:17Z suehring
44! Replace function get_topography_top_index by topo_top_ind
45!
46! 4029 2019-06-14 14:04:35Z raasch
47! nest_chemistry switch removed
48!
49! 4026 2019-06-12 16:50:15Z suehring
50! Masked topography at boundary grid points in mass conservation, in order to
51! avoid that mean velocities within topography are imposed
52!
53! 4011 2019-05-31 14:34:03Z hellstea
54! Mass (volume) flux correction included to ensure global mass conservation for child domains.
55!
56! 3987 2019-05-22 09:52:13Z kanani
57! Introduce alternative switch for debug output during timestepping
58!
59! 3984 2019-05-16 15:17:03Z hellstea
60! Commenting improved, pmci_map_fine_to_coarse_grid renamed as pmci_map_child_grid_to_parent_grid,
61! set_child_edge_coords renamed as pmci_set_child_edge_coords, some variables renamed, etc.
62!
63! 3979 2019-05-15 13:54:29Z hellstea
64! Bugfix in pmc_interp_1sto_sn. This bug had effect only in case of 1-d domain
65! decomposition with npex = 1.
66!
67! 3976 2019-05-15 11:02:34Z hellstea
68! Child initialization also for the redundant ghost points behind the nested
69! boundaries added (2nd and 3rd ghost-point layers and corners).
70!
71! 3948 2019-05-03 14:49:57Z hellstea
72! Some variables renamed, a little cleaning up and some commenting improvements
73!
74! 3947 2019-05-03 07:56:44Z hellstea
75! The checks included in 3946 are extended for the z-direction and moved into its
76! own subroutine called from pmci_define_index_mapping.
77!
78! 3946 2019-05-02 14:18:59Z hellstea
79! Check added for child domains too small in terms of number of parent-grid cells so
80! that anterpolation is not possible. Checks added for too wide anterpolation buffer
81! for the same reason. Some minor code reformatting done.
82!
83! 3945 2019-05-02 11:29:27Z raasch
84!
85! 3932 2019-04-24 17:31:34Z suehring
86! Add missing if statements for call of pmc_set_dataarray_name for TKE and
87! dissipation.
88!
89! 3888 2019-04-12 09:18:10Z hellstea
90! Variables renamed, commenting improved etc.
91!
92! 3885 2019-04-11 11:29:34Z kanani
93! Changes related to global restructuring of location messages and introduction
94! of additional debug messages
95!
96! 3883 2019-04-10 12:51:50Z hellstea
97! Checks and error messages improved and extended. All the child index bounds in the
98! parent-grid index space are made module variables. Function get_number_of_childs
99! renamed get_number_of_children. A number of variables renamed
100! and qite a lot of other code reshaping made all around the module.
101!
102! 3876 2019-04-08 18:41:49Z knoop
103! Implemented nesting for salsa variables.
104!
105! 3833 2019-03-28 15:04:04Z forkel
106! replaced USE chem_modules by USE chem_gasphase_mod
107!
108! 3822 2019-03-27 13:10:23Z hellstea
109! Temporary increase of the vertical dimension of the parent-grid arrays and
110! workarrc_t is cancelled as unnecessary.
111!
112! 3819 2019-03-27 11:01:36Z hellstea
113! Adjustable anterpolation buffer introduced on all nest boundaries, it is controlled
114! by the new nesting_parameters parameter anterpolation_buffer_width.
115!
116! 3804 2019-03-19 13:46:20Z hellstea
117! Anterpolation domain is lowered from kct-1 to kct-3 to avoid exessive       
118! kinetic energy from building up in CBL flows.
119!
120! 3803 2019-03-19 13:44:40Z hellstea
121! A bug fixed in lateral boundary interpolations. Dimension of val changed from 
122! 5 to 3 in pmci_setup_parent and pmci_setup_child.
123!
124! 3794 2019-03-15 09:36:33Z raasch
125! two remaining unused variables removed
126!
127! 3792 2019-03-14 16:50:07Z hellstea
128! Interpolations improved. Large number of obsolete subroutines removed.
129! All unused variables removed.
130!
131! 3741 2019-02-13 16:24:49Z hellstea
132! Interpolations and child initialization adjusted to handle set ups with child
133! pe-subdomain dimension not integer divisible by the grid-spacing ratio in the
134! respective direction. Set ups with pe-subdomain dimension smaller than the
135! grid-spacing ratio in the respective direction are now forbidden.
136!
137! 3708 2019-01-30 12:58:13Z hellstea
138! Checks for parent / child grid line matching introduced.
139! Interpolation of nest-boundary-tangential velocity components revised.
140!
141! 3697 2019-01-24 17:16:13Z hellstea
142! Bugfix: upper k-bound in the child initialization interpolation
143! pmci_interp_1sto_all corrected.
144! Copying of the nest boundary values into the redundant 2nd and 3rd ghost-node
145! layers is added to the pmci_interp_1sto_*-routines.
146!
147! 3681 2019-01-18 15:06:05Z hellstea
148! Linear interpolations are replaced by first order interpolations. The linear
149! interpolation routines are still included but not called. In the child
150! inititialization the interpolation is also changed to 1st order and the linear
151! interpolation is not kept.
152! Subroutine pmci_map_fine_to_coarse_grid is rewritten.
153! Several changes in pmci_init_anterp_tophat.
154! Child's parent-grid arrays (uc, vc,...) are made non-overlapping on the PE-
155! subdomain boundaries in order to allow grid-spacing ratios higher than nbgp.
156! Subroutine pmci_init_tkefactor is removed as unnecessary.
157!
158! 3655 2019-01-07 16:51:22Z knoop
159! Remove unused variable simulated_time
160!
161! 1762 2016-02-25 12:31:13Z hellstea
162! Initial revision by A. Hellsten
163!
164! Description:
165! ------------
166! Domain nesting interface routines. The low-level inter-domain communication   
167! is conducted by the PMC-library routines.
168!
169! @todo Remove array_3d variables from USE statements thate not used in the
170!       routine
171! @todo Data transfer of qc and nc is prepared but not activated
172!------------------------------------------------------------------------------!
173 MODULE pmc_interface
174
175    USE ISO_C_BINDING
176
177
178    USE arrays_3d,                                                             &
179        ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2,  &
180               pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                   &
181               u, u_p, u_2, v, v_p, v_2, w, w_p, w_2, zu, zw
182
183    USE control_parameters,                                                    &
184        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
185               bc_dirichlet_s, child_domain,                                   &
186               constant_diffusion, constant_flux_layer,                        &
187               coupling_char,                                                  &
188               debug_output_timestep,                                          &
189               dt_3d, dz, humidity, message_string,                            &
190               neutral, passive_scalar, rans_mode, rans_tke_e,                 &
191               roughness_length, salsa, topography, volume_flow
192
193    USE chem_gasphase_mod,                                                     &
194        ONLY:  nspec
195
196    USE chem_modules,                                                          &
197        ONLY:  chem_species, nesting_chem
198
199    USE chemistry_model_mod,                                                   &
200        ONLY:  spec_conc_2
201
202    USE cpulog,                                                                &
203        ONLY:  cpu_log, log_point_s
204
205    USE grid_variables,                                                        &
206        ONLY:  dx, dy
207
208    USE indices,                                                               &
209        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
210               nysv, nz, nzb, nzt, topo_top_ind, wall_flags_static_0
211
212    USE bulk_cloud_model_mod,                                                  &
213        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
214
215    USE particle_attributes,                                                   &
216        ONLY:  particle_advection
217
218    USE kinds
219
220#if defined( __parallel )
221#if !defined( __mpifh )
222    USE MPI
223#endif
224
225    USE pegrid,                                                                &
226        ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
227               numprocs, pdims, pleft, pnorth, pright, psouth, status
228
229    USE pmc_child,                                                             &
230        ONLY:  pmc_childinit, pmc_c_clear_next_array_list,                     &
231               pmc_c_getnextarray, pmc_c_get_2d_index_list, pmc_c_getbuffer,   &
232               pmc_c_putbuffer, pmc_c_setind_and_allocmem,                     &
233               pmc_c_set_dataarray, pmc_set_dataarray_name
234
235    USE pmc_general,                                                           &
236        ONLY:  da_namelen
237
238    USE pmc_handle_communicator,                                               &
239        ONLY:  pmc_get_model_info, pmc_init_model, pmc_is_rootmodel,           &
240               pmc_no_namelist_found, pmc_parent_for_child, m_couplers
241
242    USE pmc_mpi_wrapper,                                                       &
243        ONLY:  pmc_bcast, pmc_recv_from_child, pmc_recv_from_parent,           &
244               pmc_send_to_child, pmc_send_to_parent
245
246    USE pmc_parent,                                                            &
247        ONLY:  pmc_parentinit, pmc_s_clear_next_array_list, pmc_s_fillbuffer,  &
248               pmc_s_getdata_from_buffer, pmc_s_getnextarray,                  &
249               pmc_s_setind_and_allocmem, pmc_s_set_active_data_array,         &
250               pmc_s_set_dataarray, pmc_s_set_2d_index_list
251
252#endif
253   
254    USE salsa_mod,                                                             &
255        ONLY:  aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol,  &
256               ncomponents_mass, nconc_2, nesting_salsa, ngases_salsa,         &
257               salsa_gas, salsa_gases_from_chem
258
259    USE surface_mod,                                                           &
260        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
261
262    IMPLICIT NONE
263
264#if defined( __parallel )
265#if defined( __mpifh )
266    INCLUDE "mpif.h"
267#endif
268#endif
269
270    PRIVATE
271!
272!-- Constants
273    INTEGER(iwp), PARAMETER ::  child_to_parent = 2   !< Parameter for pmci_parent_datatrans indicating the direction of transfer
274    INTEGER(iwp), PARAMETER ::  parent_to_child = 1   !< Parameter for pmci_parent_datatrans indicating the direction of transfer
275    INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries
276    INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary
277
278    REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp                 !< Relative tolerence for grid-line matching tests and comparisons
279!
280!-- Coupler setup
281    INTEGER(iwp), SAVE      ::  comm_world_nesting    !< Global nesting communicator
282    INTEGER(iwp), SAVE      ::  cpl_id  = 1           !<
283    INTEGER(iwp), SAVE      ::  cpl_npe_total         !<
284    INTEGER(iwp), SAVE      ::  cpl_parent_id         !<
285   
286    CHARACTER(LEN=32), SAVE ::  cpl_name              !<
287
288!
289!-- Control parameters
290    INTEGER(iwp),     SAVE ::  anterpolation_buffer_width = 2       !< Boundary buffer width for anterpolation
291    CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !< steering parameter for data-transfer mode
292    CHARACTER(LEN=8), SAVE ::  nesting_mode = 'two-way'             !< steering parameter for 1- or 2-way nesting
293   
294    LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
295    LOGICAL       ::  rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode)
296!
297!-- Geometry
298    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_x            !< Array for the absolute x-coordinates
299    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_y            !< Array for the absolute y-coordinates
300    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_x !< x-coordinate of the lower left corner of the domain
301    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_y !< y-coordinate of the lower left corner of the domain
302!
303!-- Children's parent-grid arrays
304    INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC    ::  parent_bound        !< subdomain index bounds for children's parent-grid arrays
305
306    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< Parent-grid array on child domain - dissipation rate
307    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec    !< Parent-grid array on child domain - SGS TKE
308    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc   !< Parent-grid array on child domain - potential temperature
309    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uc    !< Parent-grid array on child domain - velocity component u
310    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vc    !< Parent-grid array on child domain - velocity component v
311    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc    !< Parent-grid array on child domain - velocity component w
312    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_c   !< Parent-grid array on child domain -
313    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc   !< Parent-grid array on child domain -
314    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qrc   !< Parent-grid array on child domain -
315    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nrc   !< Parent-grid array on child domain -
316    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc   !< Parent-grid array on child domain -
317    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc    !< Parent-grid array on child domain -
318    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  nr_partc    !<
319    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
320
321    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c  !< Parent-grid array on child domain - chemical species
322
323    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_mass_c    !< Aerosol mass
324    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_number_c  !< Aerosol number
325    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gas_c       !< SALSA gases
326!
327!-- Grid-spacing ratios.
328    INTEGER(iwp), SAVE ::  igsr    !< Integer grid-spacing ratio in i-direction
329    INTEGER(iwp), SAVE ::  jgsr    !< Integer grid-spacing ratio in j-direction
330    INTEGER(iwp), SAVE ::  kgsr    !< Integer grid-spacing ratio in k-direction
331!
332!-- Global parent-grid index bounds
333    INTEGER(iwp), SAVE ::  iplg    !< Leftmost parent-grid array ip index of the whole child domain
334    INTEGER(iwp), SAVE ::  iprg    !< Rightmost parent-grid array ip index of the whole child domain
335    INTEGER(iwp), SAVE ::  jpsg    !< Southmost parent-grid array jp index of the whole child domain
336    INTEGER(iwp), SAVE ::  jpng    !< Northmost parent-grid array jp index of the whole child domain
337!
338!-- Local parent-grid index bounds. Different sets of index bounds are needed for parent-grid arrays (uc, etc),
339!-- for index mapping arrays (iflu, etc) and for work arrays (workarr_lr, etc). This is because these arrays
340!-- have different dimensions depending on the location of the subdomain relative to boundaries and corners.
341    INTEGER(iwp), SAVE ::  ipl     !< Left index limit for children's parent-grid arrays
342    INTEGER(iwp), SAVE ::  ipla    !< Left index limit for allocation of index-mapping and other auxiliary arrays
343    INTEGER(iwp), SAVE ::  iplw    !< Left index limit for children's parent-grid work arrays
344    INTEGER(iwp), SAVE ::  ipr     !< Right index limit for children's parent-grid arrays
345    INTEGER(iwp), SAVE ::  ipra    !< Right index limit for allocation of index-mapping and other auxiliary arrays
346    INTEGER(iwp), SAVE ::  iprw    !< Right index limit for children's parent-grid work arrays
347    INTEGER(iwp), SAVE ::  jpn     !< North index limit for children's parent-grid arrays
348    INTEGER(iwp), SAVE ::  jpna    !< North index limit for allocation of index-mapping and other auxiliary arrays
349    INTEGER(iwp), SAVE ::  jpnw    !< North index limit for children's parent-grid work arrays
350    INTEGER(iwp), SAVE ::  jps     !< South index limit for children's parent-grid arrays
351    INTEGER(iwp), SAVE ::  jpsa    !< South index limit for allocation of index-mapping and other auxiliary arrays
352    INTEGER(iwp), SAVE ::  jpsw    !< South index limit for children's parent-grid work arrays
353!
354!-- Highest prognostic parent-grid k-indices.
355    INTEGER(iwp), SAVE ::  kcto     !< Upper bound for k in anterpolation of variables other than w.
356    INTEGER(iwp), SAVE ::  kctw     !< Upper bound for k in anterpolation of w.
357!
358!-- Child-array indices to be precomputed and stored for anterpolation.
359    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !< child index indicating left bound of parent grid box on u-grid
360    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !< child index indicating right bound of parent grid box on u-grid
361    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !< child index indicating left bound of parent grid box on scalar-grid
362    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !< child index indicating right bound of parent grid box on scalar-grid
363    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !< child index indicating south bound of parent grid box on v-grid
364    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !< child index indicating north bound of parent grid box on v-grid
365    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !< child index indicating south bound of parent grid box on scalar-grid
366    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !< child index indicating north bound of parent grid box on scalar-grid
367    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !< child index indicating lower bound of parent grid box on w-grid
368    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !< child index indicating upper bound of parent grid box on w-grid
369    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !< child index indicating lower bound of parent grid box on scalar-grid
370    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !< child index indicating upper bound of parent grid box on scalar-grid
371!
372!-- Number of child-grid nodes within anterpolation cells to be precomputed for anterpolation.
373    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid points contributing to a parent grid
374                                                                   !< node in anterpolation, u-grid
375    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid points contributing to a parent grid
376                                                                   !< node in anterpolation, v-grid
377    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid points contributing to a parent grid
378                                                                   !< node in anterpolation, w-grid
379    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid points contributing to a parent grid
380                                                                   !< node in anterpolation, scalar-grid
381!   
382!-- Work arrays for interpolation and user-defined type definitions for horizontal work-array exchange   
383    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_lr
384    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_sn
385    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_t
386
387    INTEGER(iwp) :: workarr_lr_exchange_type
388    INTEGER(iwp) :: workarr_sn_exchange_type
389    INTEGER(iwp) :: workarr_t_exchange_type_x
390    INTEGER(iwp) :: workarr_t_exchange_type_y
391 
392    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !< Array for communicating the parent-grid dimensions
393                                                                    !< to its children.
394
395    REAL(wp), DIMENSION(6)              ::  face_area               !< Surface area of each boundary face
396    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !< Array for communicating the real-type parent-grid
397                                                                    !< parameters to its children.
398
399    TYPE parentgrid_def
400       INTEGER(iwp)                        ::  nx                 !<
401       INTEGER(iwp)                        ::  ny                 !<
402       INTEGER(iwp)                        ::  nz                 !<
403       REAL(wp)                            ::  dx                 !<
404       REAL(wp)                            ::  dy                 !<
405       REAL(wp)                            ::  dz                 !<
406       REAL(wp)                            ::  lower_left_coord_x !<
407       REAL(wp)                            ::  lower_left_coord_y !<
408       REAL(wp)                            ::  xend               !<
409       REAL(wp)                            ::  yend               !<
410       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_x            !<
411       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_y            !<
412       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu                !<
413       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw                !<
414       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu                 !<
415       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw                 !<
416    END TYPE parentgrid_def
417
418    TYPE(parentgrid_def), SAVE, PUBLIC     ::  pg                 !< Parent-grid information package of type parentgrid_def
419!
420!-- Variables for particle coupling
421    TYPE, PUBLIC :: childgrid_def
422       INTEGER(iwp)                        ::  nx                   !<
423       INTEGER(iwp)                        ::  ny                   !<
424       INTEGER(iwp)                        ::  nz                   !<
425       REAL(wp)                            ::  dx                   !<
426       REAL(wp)                            ::  dy                   !<
427       REAL(wp)                            ::  dz                   !<
428       REAL(wp)                            ::  lx_coord, lx_coord_b !<   ! split onto separate lines
429       REAL(wp)                            ::  rx_coord, rx_coord_b !<
430       REAL(wp)                            ::  sy_coord, sy_coord_b !<
431       REAL(wp)                            ::  ny_coord, ny_coord_b !<
432       REAL(wp)                            ::  uz_coord, uz_coord_b !<
433    END TYPE childgrid_def
434
435    TYPE(childgrid_def), SAVE, ALLOCATABLE, DIMENSION(:), PUBLIC ::  childgrid  !<
436
437    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  nr_part  !<
438    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  part_adr !<
439
440   
441    INTERFACE pmci_boundary_conds
442       MODULE PROCEDURE pmci_boundary_conds
443    END INTERFACE pmci_boundary_conds
444   
445    INTERFACE pmci_check_setting_mismatches
446       MODULE PROCEDURE pmci_check_setting_mismatches
447    END INTERFACE
448
449    INTERFACE pmci_child_initialize
450       MODULE PROCEDURE pmci_child_initialize
451    END INTERFACE
452
453    INTERFACE pmci_synchronize
454       MODULE PROCEDURE pmci_synchronize
455    END INTERFACE
456
457    INTERFACE pmci_datatrans
458       MODULE PROCEDURE pmci_datatrans
459    END INTERFACE pmci_datatrans
460
461    INTERFACE pmci_ensure_nest_mass_conservation
462       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
463    END INTERFACE pmci_ensure_nest_mass_conservation
464
465    INTERFACE pmci_ensure_nest_mass_conservation_vertical
466       MODULE PROCEDURE pmci_ensure_nest_mass_conservation_vertical
467    END INTERFACE pmci_ensure_nest_mass_conservation_vertical
468
469    INTERFACE pmci_init
470       MODULE PROCEDURE pmci_init
471    END INTERFACE
472
473    INTERFACE pmci_modelconfiguration
474       MODULE PROCEDURE pmci_modelconfiguration
475    END INTERFACE
476
477    INTERFACE pmci_parent_initialize
478       MODULE PROCEDURE pmci_parent_initialize
479    END INTERFACE
480
481    INTERFACE get_number_of_children
482       MODULE PROCEDURE get_number_of_children
483    END  INTERFACE get_number_of_children
484
485    INTERFACE get_childid
486       MODULE PROCEDURE get_childid
487    END  INTERFACE get_childid
488
489    INTERFACE get_child_edges
490       MODULE PROCEDURE get_child_edges
491    END  INTERFACE get_child_edges
492
493    INTERFACE get_child_gridspacing
494       MODULE PROCEDURE get_child_gridspacing
495    END  INTERFACE get_child_gridspacing
496
497    INTERFACE pmci_set_swaplevel
498       MODULE PROCEDURE pmci_set_swaplevel
499    END INTERFACE pmci_set_swaplevel
500
501    PUBLIC child_to_parent, comm_world_nesting, cpl_id, nested_run,                                 &
502           nesting_datatransfer_mode, nesting_mode, parent_to_child, rans_mode_parent
503
504    PUBLIC pmci_boundary_conds
505    PUBLIC pmci_child_initialize
506    PUBLIC pmci_datatrans
507    PUBLIC pmci_init
508    PUBLIC pmci_modelconfiguration
509    PUBLIC pmci_parent_initialize
510    PUBLIC pmci_synchronize
511    PUBLIC pmci_set_swaplevel
512    PUBLIC get_number_of_children, get_childid, get_child_edges, get_child_gridspacing
513    PUBLIC pmci_ensure_nest_mass_conservation
514    PUBLIC pmci_ensure_nest_mass_conservation_vertical
515   
516 CONTAINS
517
518
519 SUBROUTINE pmci_init( world_comm )
520
521    USE control_parameters,                                                                         &
522        ONLY:  message_string
523
524    IMPLICIT NONE
525
526    INTEGER(iwp), INTENT(OUT) ::  world_comm   !<
527
528#if defined( __parallel )
529
530    INTEGER(iwp) ::  pmc_status   !<
531
532
533    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,                       &
534                         anterpolation_buffer_width, pmc_status )
535
536    IF ( pmc_status == pmc_no_namelist_found )  THEN
537!
538!--    This is not a nested run
539       world_comm = MPI_COMM_WORLD
540       cpl_id     = 1
541       cpl_name   = ""
542
543       RETURN
544
545    ENDIF
546!
547!-- Check steering parameter values
548    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                                                   &
549         TRIM( nesting_mode ) /= 'two-way'  .AND.                                                   &
550         TRIM( nesting_mode ) /= 'vertical' )                                                       &
551    THEN
552       message_string = 'illegal nesting mode: ' // TRIM( nesting_mode )
553       CALL message( 'pmci_init', 'PA0417', 3, 2, 0, 6, 0 )
554    ENDIF
555
556    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                                      &
557         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                                      &
558         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                                           &
559    THEN
560       message_string = 'illegal nesting datatransfer mode: ' // TRIM( nesting_datatransfer_mode )
561       CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
562    ENDIF
563!
564!-- Set the general steering switch which tells PALM that its a nested run
565    nested_run = .TRUE.
566!
567!-- Get some variables required by the pmc-interface (and in some cases in the
568!-- PALM code out of the pmci) out of the pmc-core
569    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,                               &
570                             cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,                        &
571                             cpl_name = cpl_name, npe_total = cpl_npe_total,                        &
572                             lower_left_x = lower_left_coord_x,                                     &
573                             lower_left_y = lower_left_coord_y )
574!
575!-- Set the steering switch which tells the models that they are nested (of
576!-- course the root domain is not nested)
577    IF ( .NOT.  pmc_is_rootmodel() )  THEN
578       child_domain = .TRUE.
579       WRITE( coupling_char, '(A2,I2.2)') '_N', cpl_id
580    ENDIF
581
582!
583!-- Message that communicators for nesting are initialized.
584!-- Attention: myid has been set at the end of pmc_init_model in order to
585!-- guarantee that only PE0 of the root domain does the output.
586    CALL location_message( 'initialize model nesting', 'finished' )
587!
588!-- Reset myid to its default value
589    myid = 0
590#else
591!
592!-- Nesting cannot be used in serial mode. cpl_id is set to root domain (1)
593!-- because no location messages would be generated otherwise.
594!-- world_comm is given a dummy value to avoid compiler warnings (INTENT(OUT)
595!-- must get an explicit value).
596!-- Note that this branch is only to avoid compiler warnings. The actual
597!-- execution never reaches here because the call of this subroutine is
598!-- already enclosed by  #if defined( __parallel ).
599    cpl_id     = 1
600    nested_run = .FALSE.
601    world_comm = 1
602#endif
603
604 END SUBROUTINE pmci_init
605
606
607
608 SUBROUTINE pmci_modelconfiguration
609
610    IMPLICIT NONE
611
612    INTEGER(iwp) ::  ncpl   !<  number of nest domains
613
614   
615#if defined( __parallel )
616    CALL location_message( 'setup the nested model configuration', 'start' )
617    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'start' )
618!
619!-- Compute absolute coordinates for all models
620    CALL pmci_setup_coordinates         ! CONTAIN THIS
621!
622!-- Determine the number of coupled arrays
623    CALL pmci_num_arrays                ! CONTAIN THIS
624!
625!-- Initialize the child (must be called before pmc_setup_parent)
626!-- Klaus, extend this comment to explain why it must be called before   
627    CALL pmci_setup_child               ! CONTAIN THIS
628!
629!-- Initialize PMC parent
630    CALL pmci_setup_parent              ! CONTAIN THIS
631!
632!-- Check for mismatches between settings of master and child variables
633!-- (e.g., all children have to follow the end_time settings of the root master)
634    CALL pmci_check_setting_mismatches  ! CONTAIN THIS
635!
636!-- Set flag file for combine_plot_fields for precessing the nest output data
637    OPEN( 90, FILE='3DNESTING', FORM='FORMATTED' )
638    CALL pmc_get_model_info( ncpl = ncpl )
639    WRITE( 90, '(I2)' )  ncpl
640    CLOSE( 90 )
641
642    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'stop' )
643    CALL location_message( 'setup the nested model configuration', 'finished' )
644#endif
645
646 END SUBROUTINE pmci_modelconfiguration
647
648
649
650 SUBROUTINE pmci_setup_parent
651
652#if defined( __parallel )
653    IMPLICIT NONE
654
655    INTEGER(iwp) ::  child_id           !< Child id-number for the child m
656    INTEGER(iwp) ::  ierr               !< MPI-error code
657    INTEGER(iwp) ::  kp                 !< Parent-grid index n the z-direction
658    INTEGER(iwp) ::  lb = 1             !< Running index for aerosol size bins
659    INTEGER(iwp) ::  lc = 1             !< Running index for aerosol mass bins
660    INTEGER(iwp) ::  lg = 1             !< Running index for SALSA gases
661    INTEGER(iwp) ::  m                  !< Loop index over all children of the current parent
662    INTEGER(iwp) ::  msib               !< Loop index over all other children than m in case of siblings (parallel children)
663    INTEGER(iwp) ::  n = 1              !< Running index for chemical species
664    INTEGER(iwp) ::  nest_overlap = 0   !< Tag for parallel child-domains' overlap situation (>0 if overlap found)
665    INTEGER(iwp) ::  nomatch = 0        !< Tag for child-domain mismatch situation (>0 if mismatch found)
666    INTEGER(iwp) ::  nx_child           !< Number of child-grid points in the x-direction
667    INTEGER(iwp) ::  ny_child           !< Number of child-grid points in the y-direction
668    INTEGER(iwp) ::  nz_child           !< Number of child-grid points in the z-direction
669   
670    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for receiving the child-grid dimensions from the children
671   
672    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_left   !< Minimum x-coordinate of the child domain including the ghost
673                                                           !< point layers
674    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_right  !< Maximum x-coordinate of the child domain including the ghost
675                                                           !< point layers   
676    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_south  !< Minimum y-coordinate of the child domain including the ghost
677                                                           !< point layers
678    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_north  !< Maximum y-coordinate of the child domain including the ghost
679                                                           !< point layers
680    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_x  !< Child domain x-coordinate array
681    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_y  !< Child domain y-coordinate array
682   
683    REAL(wp), DIMENSION(5) ::  child_grid_info  !< Array for receiving the child-grid spacings etc from the children
684   
685    REAL(wp) ::  child_height         !< Height of the child domain defined on the child side as zw(nzt+1)
686    REAL(wp) ::  dx_child             !< Child-grid spacing in the x-direction
687    REAL(wp) ::  dy_child             !< Child-grid spacing in the y-direction
688    REAL(wp) ::  dz_child             !< Child-grid spacing in the z-direction
689    REAL(wp) ::  left_limit           !< Left limit for the absolute x-coordinate of the child left boundary
690    REAL(wp) ::  north_limit          !< North limit for the absolute y-coordinate of the child north boundary
691    REAL(wp) ::  right_limit          !< Right limit for the absolute x-coordinate of the child right boundary
692    REAL(wp) ::  south_limit          !< South limit for the absolute y-coordinate of the child south boundary
693    REAL(wp) ::  upper_right_coord_x  !< Absolute x-coordinate of the upper right corner of the child domain
694    REAL(wp) ::  upper_right_coord_y  !< Absolute y-coordinate of the upper right corner of the child domain 
695    REAL(wp) ::  xez                  !< Minimum separation in the x-direction required between the child and
696                                      !< parent boundaries (left or right)
697    REAL(wp) ::  yez                  !< Minimum separation in the y-direction required between the child and
698                                      !< parent boundaries (south or north)
699    REAL(wp)     ::  tolex            !< Tolerance for grid-line matching in x-direction
700    REAL(wp)     ::  toley            !< Tolerance for grid-line matching in y-direction
701    REAL(wp)     ::  tolez            !< Tolerance for grid-line matching in z-direction   
702
703    CHARACTER(LEN=32) ::  myname      !< String for variable name such as 'u'
704
705    LOGICAL :: m_left_in_msib         !< Logical auxiliary parameter for the overlap test: true if the left border
706                                      !< of the child m is within the x-range of the child msib
707    LOGICAL :: m_right_in_msib        !< Logical auxiliary parameter for the overlap test: true if the right border
708                                      !< of the child m is within the x-range of the child msib
709    LOGICAL :: msib_left_in_m         !< Logical auxiliary parameter for the overlap test: true if the left border
710                                      !< of the child msib is within the x-range of the child m
711    LOGICAL :: msib_right_in_m        !< Logical auxiliary parameter for the overlap test: true if the right border
712                                      !< of the child msib is within the x-range of the child m
713    LOGICAL :: m_south_in_msib        !< Logical auxiliary parameter for the overlap test: true if the south border
714                                      !< of the child m is within the y-range of the child msib
715    LOGICAL :: m_north_in_msib        !< Logical auxiliary parameter for the overlap test: true if the north border
716                                      !< of the child m is within the y-range of the child msib
717    LOGICAL :: msib_south_in_m        !< Logical auxiliary parameter for the overlap test: true if the south border
718                                      !< of the child msib is within the y-range of the child m
719    LOGICAL :: msib_north_in_m        !< Logical auxiliary parameter for the overlap test: true if the north border
720                                      !< of the child msib is within the y-range of the child m
721
722
723!
724!-- Grid-line tolerances.
725    tolex = tolefac * dx
726    toley = tolefac * dy
727    tolez = tolefac * dz(1)   
728!
729!-- Initialize the current pmc parent.
730    CALL pmc_parentinit
731!
732!-- Corners of all children of the present parent. Note that
733!-- SIZE( pmc_parent_for_child ) = 1 if we have no children.
734    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 )  .AND.  myid == 0 )  THEN
735       ALLOCATE( child_x_left(1:SIZE( pmc_parent_for_child ) - 1) )
736       ALLOCATE( child_x_right(1:SIZE( pmc_parent_for_child ) - 1) )
737       ALLOCATE( child_y_south(1:SIZE( pmc_parent_for_child ) - 1) )
738       ALLOCATE( child_y_north(1:SIZE( pmc_parent_for_child ) - 1) )
739    ENDIF
740    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) )  THEN
741       ALLOCATE( childgrid(1:SIZE( pmc_parent_for_child ) - 1) )
742    ENDIF
743!
744!-- Get coordinates from all children and check that the children match the parent
745!-- domain and each others. Note that SIZE( pmc_parent_for_child ) = 1
746!-- if we have no children, thence the loop is not executed at all.
747    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
748
749       child_id = pmc_parent_for_child(m)
750
751       IF ( myid == 0 )  THEN
752
753          CALL pmc_recv_from_child( child_id, child_grid_dim,  SIZE(child_grid_dim), 0, 123, ierr )
754          CALL pmc_recv_from_child( child_id, child_grid_info, SIZE(child_grid_info), 0, 124, ierr )
755         
756          nx_child     = child_grid_dim(1)
757          ny_child     = child_grid_dim(2)
758          dx_child     = child_grid_info(3)
759          dy_child     = child_grid_info(4)
760          dz_child     = child_grid_info(5)
761          child_height = child_grid_info(1)
762!
763!--       Find the highest child-domain level in the parent grid for the reduced z transfer
764          DO  kp = 1, nzt                 
765             IF ( zw(kp) - child_height > tolez )  THEN                   
766                nz_child = kp
767                EXIT
768             ENDIF
769          ENDDO
770!   
771!--       Get absolute coordinates from the child
772          ALLOCATE( child_coord_x(-nbgp:nx_child+nbgp) )
773          ALLOCATE( child_coord_y(-nbgp:ny_child+nbgp) )
774         
775          CALL pmc_recv_from_child( child_id, child_coord_x, SIZE( child_coord_x ), 0, 11, ierr )
776          CALL pmc_recv_from_child( child_id, child_coord_y, SIZE( child_coord_y ), 0, 12, ierr )
777         
778          parent_grid_info_real(1) = lower_left_coord_x
779          parent_grid_info_real(2) = lower_left_coord_y
780          parent_grid_info_real(3) = dx
781          parent_grid_info_real(4) = dy
782
783          upper_right_coord_x      = lower_left_coord_x + ( nx + 1 ) * dx
784          upper_right_coord_y      = lower_left_coord_y + ( ny + 1 ) * dy
785          parent_grid_info_real(5) = upper_right_coord_x
786          parent_grid_info_real(6) = upper_right_coord_y
787          parent_grid_info_real(7) = dz(1)
788
789          parent_grid_info_int(1)  = nx
790          parent_grid_info_int(2)  = ny
791          parent_grid_info_int(3)  = nz_child
792!
793!--       Check that the child domain matches its parent domain.
794          IF ( nesting_mode == 'vertical' )  THEN
795!
796!--          In case of vertical nesting, the lateral boundaries must match exactly.
797             right_limit = upper_right_coord_x
798             north_limit = upper_right_coord_y
799             IF ( ( ABS( child_coord_x(nx_child+1) - right_limit ) > tolex )  .OR.                  &
800                  ( ABS( child_coord_y(ny_child+1) - north_limit ) > toley ) )  THEN               
801                nomatch = 1
802             ENDIF
803          ELSE       
804!
805!--          In case of 3-D nesting, check that the child domain is completely
806!--          inside its parent domain.
807             xez = ( nbgp + 1 ) * dx
808             yez = ( nbgp + 1 ) * dy
809             left_limit  = lower_left_coord_x + xez
810             right_limit = upper_right_coord_x - xez
811             south_limit = lower_left_coord_y + yez
812             north_limit = upper_right_coord_y - yez
813             IF ( ( left_limit - child_coord_x(0) > tolex )  .OR.                                   &
814                  ( child_coord_x(nx_child+1) - right_limit > tolex )  .OR.                         &
815                  ( south_limit - child_coord_y(0) > toley )  .OR.                                  &
816                  ( child_coord_y(ny_child+1) - north_limit > toley ) )  THEN
817                nomatch = 1
818             ENDIF
819          ENDIF
820!             
821!--       Child domain must be lower than the parent domain such that the top ghost
822!--       layer of the child grid does not exceed the parent domain top boundary.
823          IF ( child_height - zw(nzt) > tolez ) THEN
824             nomatch = 1
825          ENDIF
826!
827!--       If parallel child domains (siblings) do exist ( m > 1 ),
828!--       check that they do not overlap.
829          child_x_left(m)  = child_coord_x(-nbgp)
830          child_x_right(m) = child_coord_x(nx_child+nbgp)
831          child_y_south(m) = child_coord_y(-nbgp)
832          child_y_north(m) = child_coord_y(ny_child+nbgp)
833
834          IF ( nesting_mode /= 'vertical' )  THEN
835!
836!--          Note that the msib-loop is executed only if ( m > 1 ). 
837!--          Also note that the tests have to be made both ways (m vs msib and msib vs m)
838!--          in order to detect all the possible overlap situations.
839             DO  msib = 1, m - 1
840!
841!--             Set some logical auxiliary parameters to simplify the IF-condition.                 
842                m_left_in_msib  = ( child_x_left(m)  >= child_x_left(msib)  - tolex )  .AND.        &
843                                  ( child_x_left(m)  <= child_x_right(msib) + tolex )
844                m_right_in_msib = ( child_x_right(m) >= child_x_left(msib)  - tolex )  .AND.        &
845                                  ( child_x_right(m) <= child_x_right(msib) + tolex )
846                msib_left_in_m  = ( child_x_left(msib)  >= child_x_left(m)  - tolex )  .AND.        &
847                                  ( child_x_left(msib)  <= child_x_right(m) + tolex )
848                msib_right_in_m = ( child_x_right(msib) >= child_x_left(m)  - tolex )  .AND.        &
849                                  ( child_x_right(msib) <= child_x_right(m) + tolex )
850                m_south_in_msib = ( child_y_south(m) >= child_y_south(msib) - toley )  .AND.        &
851                                  ( child_y_south(m) <= child_y_north(msib) + toley )
852                m_north_in_msib = ( child_y_north(m) >= child_y_south(msib) - toley )  .AND.        &
853                                  ( child_y_north(m) <= child_y_north(msib) + toley )
854                msib_south_in_m = ( child_y_south(msib) >= child_y_south(m) - toley )  .AND.        &
855                                  ( child_y_south(msib) <= child_y_north(m) + toley )
856                msib_north_in_m = ( child_y_north(msib) >= child_y_south(m) - toley )  .AND.        &
857                                  ( child_y_north(msib) <= child_y_north(m) + toley )
858               
859                IF ( ( m_left_in_msib  .OR.  m_right_in_msib  .OR.                                  &
860                       msib_left_in_m  .OR.  msib_right_in_m )                                      &
861                     .AND.                                                                          &
862                     ( m_south_in_msib  .OR.  m_north_in_msib  .OR.                                 &
863                       msib_south_in_m  .OR.  msib_north_in_m ) )  THEN
864                     nest_overlap = 1
865                ENDIF
866
867             ENDDO
868          ENDIF         
869
870          CALL pmci_set_child_edge_coords
871
872          DEALLOCATE( child_coord_x )
873          DEALLOCATE( child_coord_y )
874!
875!--       Send information about operating mode (LES or RANS) to child. This will be
876!--       used to control TKE nesting and setting boundary conditions properly.
877          CALL pmc_send_to_child( child_id, rans_mode, 1, 0, 19, ierr ) 
878!
879!--       Send parent grid information to child
880          CALL pmc_send_to_child( child_id, parent_grid_info_real,                                  &
881                                  SIZE( parent_grid_info_real ), 0, 21,                             &
882                                  ierr )
883          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,                            &
884                                  22, ierr )
885!
886!--       Send local grid to child
887          CALL pmc_send_to_child( child_id, coord_x, nx+1+2*nbgp, 0, 24,                            &
888                                  ierr )
889          CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,                            &
890                                  ierr )
891!
892!--       Also send the dzu-, dzw-, zu- and zw-arrays here
893          CALL pmc_send_to_child( child_id, dzu, nz_child + 1, 0, 26, ierr )
894          CALL pmc_send_to_child( child_id, dzw, nz_child + 1, 0, 27, ierr )
895          CALL pmc_send_to_child( child_id, zu,  nz_child + 2, 0, 28, ierr )
896          CALL pmc_send_to_child( child_id, zw,  nz_child + 2, 0, 29, ierr )
897
898          IF ( nomatch /= 0 )  THEN
899             WRITE ( message_string, * ) 'nested child domain does not fit into its parent domain'
900             CALL message( 'pmci_setup_parent', 'PA0425', 3, 2, 0, 6, 0 )
901          ENDIF
902 
903          IF ( nest_overlap /= 0  .AND.  nesting_mode /= 'vertical' )  THEN
904             WRITE ( message_string, * ) 'nested parallel child domains overlap'
905             CALL message( 'pmci_setup_parent', 'PA0426', 3, 2, 0, 6, 0 )
906          ENDIF
907         
908       ENDIF  ! ( myid == 0 )
909
910       CALL MPI_BCAST( nz_child, 1, MPI_INTEGER, 0, comm2d, ierr )
911
912       CALL MPI_BCAST( childgrid(m), STORAGE_SIZE(childgrid(1))/8, MPI_BYTE, 0, comm2d, ierr )
913!
914!--    Set up the index-list which is an integer array that maps the child index space on
915!--    the parent index- and subdomain spaces.
916       CALL pmci_create_index_list
917!
918!--    Include couple arrays into parent content.
919!--    The adresses of the PALM 2D or 3D array (here parent grid) which are candidates
920!--    for coupling are stored once into the pmc context. While data transfer, the array do not
921!--    have to be specified again
922       CALL pmc_s_clear_next_array_list
923       DO WHILE ( pmc_s_getnextarray( child_id, myname ) )
924          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
925             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = n )
926             n = n + 1 
927          ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
928             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lb )
929             lb = lb + 1 
930          ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
931             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lc )
932             lc = lc + 1 
933          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.  .NOT. salsa_gases_from_chem )  THEN
934             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lg )
935             lg = lg + 1
936          ELSE
937             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child )
938          ENDIF
939       ENDDO
940
941       CALL pmc_s_setind_and_allocmem( child_id )
942       
943    ENDDO  ! m
944
945    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN
946       DEALLOCATE( child_x_left )
947       DEALLOCATE( child_x_right )
948       DEALLOCATE( child_y_south )
949       DEALLOCATE( child_y_north )
950    ENDIF
951
952   
953 CONTAINS
954
955
956    SUBROUTINE pmci_create_index_list
957
958       IMPLICIT NONE
959
960       INTEGER(iwp) ::  ilist            !< Index-list index running over the child's parent-grid jc,ic-space
961       INTEGER(iwp) ::  index_list_size  !< Dimension 2 of the array index_list
962       INTEGER(iwp) ::  ierr             !< MPI error code
963       INTEGER(iwp) ::  ip               !< Running parent-grid index on the child domain in the x-direction
964       INTEGER(iwp) ::  jp               !< Running parent-grid index on the child domain in the y-direction
965       INTEGER(iwp) ::  n                !< Running index over child subdomains
966       INTEGER(iwp) ::  nrx              !< Parent subdomain dimension in the x-direction
967       INTEGER(iwp) ::  nry              !< Parent subdomain dimension in the y-direction
968       INTEGER(iwp) ::  pex              !< Two-dimensional subdomain (pe) index in the x-direction
969       INTEGER(iwp) ::  pey              !< Two-dimensional subdomain (pe) index in the y-direction
970       INTEGER(iwp) ::  parent_pe        !< Parent subdomain index (one-dimensional)
971
972       INTEGER(iwp), DIMENSION(2) ::  pe_indices_2d                                  !< Array for two-dimensional subdomain (pe)
973                                                                                     !< indices needed for MPI_CART_RANK
974       INTEGER(iwp), DIMENSION(2) ::  size_of_childs_parent_grid_bounds_all          !< Dimensions of childs_parent_grid_bounds_all
975       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  childs_parent_grid_bounds_all  !< Array that contains the child's
976                                                                                     !< parent-grid index bounds for all its
977                                                                                     !< subdomains (pes)
978       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list                     !< Array that maps the child index space on
979                                                                                     !< the parent index- and subdomain spaces
980       
981       IF ( myid == 0 )  THEN
982         
983          CALL pmc_recv_from_child( child_id, size_of_childs_parent_grid_bounds_all,                &
984                                    2, 0, 40, ierr )
985          ALLOCATE( childs_parent_grid_bounds_all(size_of_childs_parent_grid_bounds_all(1),         &
986                                                  size_of_childs_parent_grid_bounds_all(2)) )
987          CALL pmc_recv_from_child( child_id, childs_parent_grid_bounds_all,                        &
988                                    SIZE( childs_parent_grid_bounds_all ), 0, 41, ierr )
989!
990!--       Compute size (dimension) of the index_list.
991          index_list_size = 0         
992          DO  n = 1, size_of_childs_parent_grid_bounds_all(2)
993             index_list_size = index_list_size +                                                    &
994                  ( childs_parent_grid_bounds_all(4,n) - childs_parent_grid_bounds_all(3,n) + 1 ) * &
995                  ( childs_parent_grid_bounds_all(2,n) - childs_parent_grid_bounds_all(1,n) + 1 )
996          ENDDO
997
998          ALLOCATE( index_list(6,index_list_size) )
999
1000          nrx = nxr - nxl + 1
1001          nry = nyn - nys + 1
1002          ilist = 0
1003!
1004!--       Loop over all children PEs
1005          DO  n = 1, size_of_childs_parent_grid_bounds_all(2)           !
1006!
1007!--          Subspace along y required by actual child PE
1008             DO  jp = childs_parent_grid_bounds_all(3,n), childs_parent_grid_bounds_all(4,n)  ! jp = jps, jpn of child PE# n
1009!
1010!--             Subspace along x required by actual child PE
1011                DO  ip = childs_parent_grid_bounds_all(1,n), childs_parent_grid_bounds_all(2,n)  ! ip = ipl, ipr of child PE# n
1012
1013                   pex = ip / nrx
1014                   pey = jp / nry
1015                   pe_indices_2d(1) = pex
1016                   pe_indices_2d(2) = pey
1017                   CALL MPI_CART_RANK( comm2d, pe_indices_2d, parent_pe, ierr )
1018                 
1019                   ilist = ilist + 1
1020!
1021!--                First index in parent array  ! TO_DO: Klaus, please explain better
1022                   index_list(1,ilist) = ip - ( pex * nrx ) + 1 + nbgp
1023!
1024!--                Second index in parent array  ! TO_DO: Klaus, please explain better
1025                   index_list(2,ilist) = jp - ( pey * nry ) + 1 + nbgp
1026!
1027!--                x index of child's parent grid
1028                   index_list(3,ilist) = ip - childs_parent_grid_bounds_all(1,n) + 1
1029!
1030!--                y index of child's parent grid
1031                   index_list(4,ilist) = jp - childs_parent_grid_bounds_all(3,n) + 1
1032!
1033!--                PE number of child
1034                   index_list(5,ilist) = n - 1
1035!
1036!--                PE number of parent
1037                   index_list(6,ilist) = parent_pe
1038
1039                ENDDO
1040             ENDDO
1041          ENDDO
1042!
1043!--       TO_DO: Klaus: comment what is done here
1044          CALL pmc_s_set_2d_index_list( child_id, index_list(:,1:ilist) )
1045
1046       ELSE
1047!
1048!--       TO_DO: Klaus: comment why this dummy allocation is required
1049          ALLOCATE( index_list(6,1) )
1050          CALL pmc_s_set_2d_index_list( child_id, index_list )
1051       ENDIF
1052
1053       DEALLOCATE(index_list)
1054
1055     END SUBROUTINE pmci_create_index_list
1056
1057
1058
1059     SUBROUTINE pmci_set_child_edge_coords
1060        IMPLICIT  NONE
1061
1062        INTEGER(iwp) ::  nbgp_lpm = 1  !< Number of ghost-point layers used for lpm (Klaus, is this correct?)
1063
1064       
1065        nbgp_lpm = MIN( nbgp_lpm, nbgp )
1066
1067        childgrid(m)%nx = nx_child
1068        childgrid(m)%ny = ny_child
1069        childgrid(m)%nz = nz_child
1070        childgrid(m)%dx = dx_child
1071        childgrid(m)%dy = dy_child
1072        childgrid(m)%dz = dz_child
1073
1074        childgrid(m)%lx_coord   = child_coord_x(0)
1075        childgrid(m)%lx_coord_b = child_coord_x(-nbgp_lpm)
1076        childgrid(m)%rx_coord   = child_coord_x(nx_child) + dx_child
1077        childgrid(m)%rx_coord_b = child_coord_x(nx_child+nbgp_lpm) + dx_child
1078        childgrid(m)%sy_coord   = child_coord_y(0)
1079        childgrid(m)%sy_coord_b = child_coord_y(-nbgp_lpm)
1080        childgrid(m)%ny_coord   = child_coord_y(ny_child) + dy_child
1081        childgrid(m)%ny_coord_b = child_coord_y(ny_child+nbgp_lpm) + dy_child
1082        childgrid(m)%uz_coord   = child_grid_info(2)
1083        childgrid(m)%uz_coord_b = child_grid_info(1)
1084
1085     END SUBROUTINE pmci_set_child_edge_coords
1086
1087#endif
1088 END SUBROUTINE pmci_setup_parent
1089
1090
1091
1092 SUBROUTINE pmci_setup_child
1093
1094#if defined( __parallel )
1095    IMPLICIT NONE
1096
1097    INTEGER(iwp) ::  ierr                          !< MPI error code
1098    INTEGER(iwp) ::  lb                            !< Running index for aerosol size bins
1099    INTEGER(iwp) ::  lc                            !< Running index for aerosol mass bins
1100    INTEGER(iwp) ::  lg                            !< Running index for SALSA gases
1101    INTEGER(iwp) ::  n                             !< Running index for number of chemical species
1102    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for sending the child-grid dimensions to parent
1103
1104    REAL(wp), DIMENSION(5) ::  child_grid_info     !< Array for sending the child-grid spacings etc to parent
1105         
1106    CHARACTER( LEN=da_namelen ) ::  myname         !< Name of the variable to be coupled
1107    CHARACTER(LEN=5) ::  salsa_char                !< Name extension for the variable name in case of SALSA variable
1108   
1109!
1110!-- Child setup
1111!-- Root model does not have a parent and is not a child, therefore no child setup on root model
1112    IF ( .NOT. pmc_is_rootmodel() )  THEN
1113!
1114!--    KLaus, add a description here what pmc_childinit does       
1115       CALL pmc_childinit
1116!
1117!--    The arrays, which actually will be exchanged between child and parent
1118!--    are defined Here AND ONLY HERE.
1119!--    If a variable is removed, it only has to be removed from here.
1120!--    Please check, if the arrays are in the list of POSSIBLE exchange arrays
1121!--    in subroutines:
1122!--    pmci_set_array_pointer (for parent arrays)
1123!--    pmci_create_childs_parent_grid_arrays (for child's parent-grid arrays)
1124       CALL pmc_set_dataarray_name( 'parent', 'u', 'child', 'u', ierr )
1125       CALL pmc_set_dataarray_name( 'parent', 'v', 'child', 'v', ierr )
1126       CALL pmc_set_dataarray_name( 'parent', 'w', 'child', 'w', ierr )
1127!
1128!--    Set data array name for TKE. Please note, nesting of TKE is actually
1129!--    only done if both parent and child are in LES or in RANS mode. Due to
1130!--    design of model coupler, however, data array names must be already
1131!--    available at this point.
1132       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                              &
1133            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                               &
1134               .NOT. constant_diffusion ) )  THEN
1135          CALL pmc_set_dataarray_name( 'parent', 'e', 'child', 'e', ierr )
1136       ENDIF
1137!
1138!--    Nesting of dissipation rate only if both parent and child are in RANS
1139!--    mode and TKE-epsilon closure is applied. Please see also comment for TKE
1140!--    above.
1141       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
1142          CALL pmc_set_dataarray_name( 'parent', 'diss', 'child', 'diss', ierr )
1143       ENDIF
1144
1145       IF ( .NOT. neutral )  THEN
1146          CALL pmc_set_dataarray_name( 'parent', 'pt' ,'child', 'pt', ierr )
1147       ENDIF
1148
1149       IF ( humidity )  THEN
1150
1151          CALL pmc_set_dataarray_name( 'parent', 'q', 'child', 'q', ierr )
1152
1153          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1154            CALL pmc_set_dataarray_name( 'parent', 'qc', 'child', 'qc', ierr ) 
1155            CALL pmc_set_dataarray_name( 'parent', 'nc', 'child', 'nc', ierr ) 
1156          ENDIF
1157
1158          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1159             CALL pmc_set_dataarray_name( 'parent', 'qr', 'child', 'qr', ierr )
1160             CALL pmc_set_dataarray_name( 'parent', 'nr', 'child', 'nr', ierr ) 
1161          ENDIF
1162     
1163       ENDIF
1164
1165       IF ( passive_scalar )  THEN
1166          CALL pmc_set_dataarray_name( 'parent', 's', 'child', 's', ierr )
1167       ENDIF
1168
1169       IF ( particle_advection )  THEN
1170          CALL pmc_set_dataarray_name( 'parent', 'nr_part', 'child', 'nr_part', ierr )
1171          CALL pmc_set_dataarray_name( 'parent', 'part_adr', 'child', 'part_adr', ierr )
1172       ENDIF
1173       
1174       IF ( air_chemistry  .AND.  nesting_chem )  THEN
1175          DO n = 1, nspec
1176             CALL pmc_set_dataarray_name( 'parent', 'chem_' // TRIM( chem_species(n)%name ),        &
1177                                          'child',  'chem_' // TRIM( chem_species(n)%name ), ierr )
1178          ENDDO
1179       ENDIF
1180
1181       IF ( salsa  .AND.  nesting_salsa )  THEN
1182          DO  lb = 1, nbins_aerosol
1183             WRITE(salsa_char,'(i0)') lb
1184             CALL pmc_set_dataarray_name( 'parent', 'an_' // TRIM( salsa_char ),                    &
1185                                          'child',  'an_' // TRIM( salsa_char ), ierr )
1186          ENDDO
1187          DO  lc = 1, nbins_aerosol * ncomponents_mass
1188             WRITE(salsa_char,'(i0)') lc
1189             CALL pmc_set_dataarray_name( 'parent', 'am_' // TRIM( salsa_char ),                    &
1190                                          'child',  'am_' // TRIM( salsa_char ), ierr )
1191          ENDDO
1192          IF ( .NOT. salsa_gases_from_chem )  THEN
1193             DO  lg = 1, ngases_salsa
1194                WRITE(salsa_char,'(i0)') lg
1195                CALL pmc_set_dataarray_name( 'parent', 'sg_' // TRIM( salsa_char ),                 &
1196                                             'child',  'sg_' // TRIM( salsa_char ), ierr )
1197             ENDDO
1198          ENDIF
1199       ENDIF
1200
1201       CALL pmc_set_dataarray_name( lastentry = .TRUE. )
1202!
1203!--    Send grid to parent
1204       child_grid_dim(1)  = nx
1205       child_grid_dim(2)  = ny
1206       child_grid_dim(3)  = nz
1207       child_grid_info(1) = zw(nzt+1)
1208       child_grid_info(2) = zw(nzt)
1209       child_grid_info(3) = dx
1210       child_grid_info(4) = dy
1211       child_grid_info(5) = dz(1)
1212
1213       IF ( myid == 0 )  THEN
1214
1215          CALL pmc_send_to_parent( child_grid_dim, SIZE( child_grid_dim ), 0, 123, ierr )
1216          CALL pmc_send_to_parent( child_grid_info, SIZE( child_grid_info ), 0, 124, ierr )
1217          CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
1218          CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
1219
1220          CALL pmc_recv_from_parent( rans_mode_parent, 1, 0, 19, ierr )
1221!
1222!--       Receive parent-grid information.
1223          CALL pmc_recv_from_parent( parent_grid_info_real,                    &
1224                                     SIZE(parent_grid_info_real), 0, 21, ierr )
1225          CALL pmc_recv_from_parent( parent_grid_info_int,  3, 0, 22, ierr )
1226
1227       ENDIF
1228
1229       CALL MPI_BCAST( parent_grid_info_real, SIZE(parent_grid_info_real),     &
1230                       MPI_REAL, 0, comm2d, ierr )
1231       CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr )
1232
1233       pg%dx = parent_grid_info_real(3)
1234       pg%dy = parent_grid_info_real(4)
1235       pg%dz = parent_grid_info_real(7)
1236       pg%nx = parent_grid_info_int(1)
1237       pg%ny = parent_grid_info_int(2)
1238       pg%nz = parent_grid_info_int(3)
1239!
1240!--    Allocate 1-D arrays for parent-grid coordinates and grid-spacings in the z-direction
1241       ALLOCATE( pg%coord_x(-nbgp:pg%nx+nbgp) )
1242       ALLOCATE( pg%coord_y(-nbgp:pg%ny+nbgp) )
1243       ALLOCATE( pg%dzu(1:pg%nz+1) )
1244       ALLOCATE( pg%dzw(1:pg%nz+1) )
1245       ALLOCATE( pg%zu(0:pg%nz+1) )
1246       ALLOCATE( pg%zw(0:pg%nz+1) )
1247!
1248!--    Get parent-grid coordinates and grid-spacings in the z-direction from the parent
1249       IF ( myid == 0)  THEN
1250          CALL pmc_recv_from_parent( pg%coord_x, pg%nx+1+2*nbgp, 0, 24, ierr )
1251          CALL pmc_recv_from_parent( pg%coord_y, pg%ny+1+2*nbgp, 0, 25, ierr )
1252          CALL pmc_recv_from_parent( pg%dzu, pg%nz+1, 0, 26, ierr )
1253          CALL pmc_recv_from_parent( pg%dzw, pg%nz+1, 0, 27, ierr )
1254          CALL pmc_recv_from_parent( pg%zu, pg%nz+2, 0, 28, ierr )
1255          CALL pmc_recv_from_parent( pg%zw, pg%nz+2, 0, 29, ierr )
1256       ENDIF
1257!
1258!--    Broadcast this information
1259       CALL MPI_BCAST( pg%coord_x, pg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1260       CALL MPI_BCAST( pg%coord_y, pg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1261       CALL MPI_BCAST( pg%dzu, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
1262       CALL MPI_BCAST( pg%dzw, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
1263       CALL MPI_BCAST( pg%zu, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
1264       CALL MPI_BCAST( pg%zw, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
1265       CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )       
1266!
1267!--    Find the index bounds for the nest domain in the parent-grid index space
1268       CALL pmci_map_child_grid_to_parent_grid
1269!
1270!--    TO_DO: Klaus give a comment what is happening here
1271       CALL pmc_c_get_2d_index_list
1272!
1273!--    Include couple arrays into child content
1274!--    TO_DO: Klaus: better explain the above comment (what is child content?)
1275       CALL  pmc_c_clear_next_array_list
1276
1277       n  = 1
1278       lb = 1
1279       lc = 1
1280       lg = 1
1281
1282       DO  WHILE ( pmc_c_getnextarray( myname ) )
1283!
1284!--       Note that pg%nz is not the original nz of parent, but the highest
1285!--       parent-grid level needed for nesting.
1286!--       Note that in case of chemical species or SALSA variables an additional
1287!--       parameter needs to be passed. The parameter is required to set the pointer
1288!--       correctlyto the chemical-species or SALSA data structure. Hence, first check if
1289!--       the current variable is a chemical species or a SALSA variable. If so, pass
1290!--       index id of respective sub-variable (species or bin) and increment this subsequently.
1291          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
1292             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, n )
1293             n = n + 1   
1294          ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
1295             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lb )
1296             lb = lb + 1
1297          ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
1298             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lc )
1299             lc = lc + 1
1300          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.  .NOT.  salsa_gases_from_chem )  THEN
1301             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lg )
1302             lg = lg + 1
1303          ELSE
1304             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz )
1305          ENDIF
1306       ENDDO
1307       CALL pmc_c_setind_and_allocmem
1308!
1309!--    Precompute the index-mapping arrays
1310       CALL pmci_define_index_mapping
1311!
1312!--    Check that the child and parent grid lines do match
1313       CALL pmci_check_grid_matching
1314!       
1315!--    Compute surface areas of the nest-boundary faces
1316       CALL pmci_compute_face_areas
1317       
1318    ENDIF
1319
1320 CONTAINS
1321
1322
1323    SUBROUTINE pmci_map_child_grid_to_parent_grid
1324!
1325!--    Determine index bounds of interpolation/anterpolation area in the parent-grid index space
1326       IMPLICIT NONE
1327
1328       INTEGER(iwp), DIMENSION(5,numprocs) ::  parent_bound_all     !< Transfer array for parent-grid index bounds
1329
1330       INTEGER(iwp), DIMENSION(4)          ::  parent_bound_global  !< Transfer array for global parent-grid index bounds
1331       INTEGER(iwp), DIMENSION(2)          ::  size_of_array        !< For sending the dimensions of parent_bound_all to parent
1332
1333       INTEGER(iwp) ::  ip      !< Running parent-grid index in the x-direction
1334       INTEGER(iwp) ::  iauxl   !< Offset between the index bound ipl and the auxiliary index bound ipla
1335       INTEGER(iwp) ::  iauxr   !< Offset between the index bound ipr and the auxiliary index bound ipra
1336       INTEGER(iwp) ::  ijaux   !< Temporary variable for receiving the index bound from the neighbouring subdomain
1337       INTEGER(iwp) ::  jp      !< Running parent-grid index in the y-direction
1338       INTEGER(iwp) ::  jauxs   !< Offset between the index bound jps and the auxiliary index bound jpsa
1339       INTEGER(iwp) ::  jauxn   !< Offset between the index bound jpn and the auxiliary index bound jpna
1340
1341       REAL(wp) ::  tolex       !< Tolerance for grid-line matching in x-direction   
1342       REAL(wp) ::  toley       !< Tolerance for grid-line matching in y-direction   
1343       REAL(wp) ::  xexl        !< Parent-grid array exceedance behind the left edge of the child PE subdomain
1344       REAL(wp) ::  xexr        !< Parent-grid array exceedance behind the right edge of the child PE subdomain
1345       REAL(wp) ::  yexs        !< Parent-grid array exceedance behind the south edge of the child PE subdomain
1346       REAL(wp) ::  yexn        !< Parent-grid array exceedance behind the north edge of the child PE subdomain
1347       REAL(wp) ::  xpl         !< Requested left-edge x-coordinate of the parent-grid array domain (at the internal boundaries
1348                                !< the real edge may differ from this in some cases as explained in the comment block below) 
1349       REAL(wp) ::  xpr         !< Requested right-edge x-coordinate of the parent-grid array domain (at the internal boundaries
1350                                !< the real edge may differ from this in some cases as explained in the comment block below)
1351       REAL(wp) ::  yps         !< Requested south-edge y-coordinate of the parent-grid array domain (at the internal boundaries
1352                                !< the real edge may differ from this in some cases as explained in the comment block below)
1353       REAL(wp) ::  ypn         !< Requested south-edge y-coordinate of the parent-grid array domain (at the internal boundaries
1354                                !< the real edge may differ from this in some cases as explained in the comment block below)
1355
1356!
1357!--    Determine the index limits for the child's parent-grid arrays (such as uc for example).
1358!--    Note that at the outer edges of the child domain (nest boundaries) these arrays exceed
1359!--    the boundary by two parent-grid cells. At the internal boundaries, there are no
1360!--    exceedances and thus no overlaps with the neighbouring subdomain. If at least half
1361!--    of the parent-grid cell is within the current child sub-domain, then it is included
1362!--    in the current sub-domain's parent-grid array. Else the parent-grid cell is
1363!--    included in the neighbouring subdomain's parent-grid array, or not included at all if
1364!--    we are at the outer edge of the child domain. This may occur especially when a large
1365!--    grid-spacing ratio is used.
1366!
1367!--    Tolerances for grid-line matching.
1368       tolex = tolefac * dx
1369       toley = tolefac * dy
1370!
1371!--    Left boundary.
1372!--    Extension by two parent-grid cells behind the boundary, see the comment block above.   
1373       IF ( bc_dirichlet_l )  THEN
1374          xexl  = 2.0_wp * pg%dx
1375          iauxl = 0
1376       ELSE
1377          xexl  = 0.0_wp
1378          iauxl = 1
1379       ENDIF
1380       xpl     = coord_x(nxl) - xexl
1381       DO  ip = 0, pg%nx
1382          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl - tolex )  THEN
1383             ipl = MAX( 0, ip )
1384             EXIT
1385          ENDIF
1386       ENDDO
1387!
1388!--    Right boundary.
1389!--    Extension by two parent-grid cells behind the boundary, see the comment block above.       
1390       IF ( bc_dirichlet_r )  THEN
1391          xexr  = 2.0_wp * pg%dx
1392          iauxr = 0 
1393       ELSE
1394          xexr  = 0.0_wp
1395          iauxr = 1 
1396       ENDIF
1397       xpr  = coord_x(nxr+1) + xexr
1398       DO  ip = pg%nx, 0 , -1
1399          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr + tolex )  THEN
1400             ipr = MIN( pg%nx, MAX( ipl, ip ) )
1401             EXIT
1402          ENDIF
1403       ENDDO
1404!
1405!--    South boundary.
1406!--    Extension by two parent-grid cells behind the boundary, see the comment block above.   
1407       IF ( bc_dirichlet_s )  THEN
1408          yexs  = 2.0_wp * pg%dy
1409          jauxs = 0 
1410       ELSE
1411          yexs  = 0.0_wp
1412          jauxs = 1 
1413       ENDIF
1414       yps  = coord_y(nys) - yexs
1415       DO  jp = 0, pg%ny
1416          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps - toley )  THEN             
1417             jps = MAX( 0, jp )
1418             EXIT
1419          ENDIF
1420       ENDDO
1421!
1422!--    North boundary.
1423!--    Extension by two parent-grid cells behind the boundary, see the comment block above. 
1424       IF  ( bc_dirichlet_n )  THEN
1425          yexn  = 2.0_wp * pg%dy
1426          jauxn = 0
1427       ELSE
1428          yexn  = 0.0_wp
1429          jauxn = 1
1430       ENDIF
1431       ypn  = coord_y(nyn+1) + yexn
1432       DO  jp = pg%ny, 0 , -1
1433          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn + toley )  THEN
1434             jpn = MIN( pg%ny, MAX( jps, jp ) )
1435             EXIT
1436          ENDIF
1437       ENDDO
1438!
1439!--    Make sure that the indexing is contiguous (no gaps, no overlaps).
1440!--    This is a safety measure mainly for cases with high grid-spacing
1441!--    ratio and narrow child subdomains.
1442       IF ( pdims(1) > 1 )  THEN
1443          IF ( nxl == 0 )  THEN
1444             CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
1445          ELSE IF ( nxr == nx )  THEN
1446             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr )
1447             ipl = ijaux + 1
1448          ELSE
1449             CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
1450             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 
1451             ipl = ijaux + 1
1452          ENDIF
1453       ENDIF
1454
1455       IF ( pdims(2) > 1 )  THEN
1456          IF ( nys == 0 )  THEN
1457             CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
1458          ELSE IF ( nyn == ny )  THEN
1459             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr )
1460             jps = ijaux + 1
1461          ELSE
1462             CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
1463             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 
1464             jps = ijaux + 1
1465          ENDIF
1466       ENDIF
1467         
1468       WRITE(9,"('pmci_map_child_grid_to_parent_grid. Parent-grid array bounds: ',4(i4,2x))")             &
1469            ipl, ipr, jps, jpn
1470       FLUSH(9)
1471
1472       parent_bound(1) = ipl
1473       parent_bound(2) = ipr
1474       parent_bound(3) = jps
1475       parent_bound(4) = jpn
1476       parent_bound(5) = myid
1477!
1478!--    The following auxiliary index bounds are used for allocating index mapping and
1479!--    some other auxiliary arrays.
1480       ipla = ipl - iauxl
1481       ipra = ipr + iauxr
1482       jpsa = jps - jauxs
1483       jpna = jpn + jauxn
1484!
1485!--    The index-bounds parent_bound of all subdomains of the current child domain
1486!--    must be sent to the parent in order for the parent to create the index list.
1487!--    For this reason, the parent_bound arrays are packed together in single
1488!--    array parent_bound_all using MPI_GATHER.       
1489!--    Note that MPI_Gather receives data from all processes in the rank order
1490!--    This fact is exploited in creating the index list in pmci_create_index_list.
1491       CALL MPI_GATHER( parent_bound, 5, MPI_INTEGER, parent_bound_all, 5,                          &
1492                        MPI_INTEGER, 0, comm2d, ierr )
1493
1494       IF ( myid == 0 )  THEN
1495          size_of_array(1) = SIZE( parent_bound_all, 1 )
1496          size_of_array(2) = SIZE( parent_bound_all, 2 )
1497          CALL pmc_send_to_parent( size_of_array, 2, 0, 40, ierr )
1498          CALL pmc_send_to_parent( parent_bound_all, SIZE( parent_bound_all ), 0, 41, ierr )
1499!
1500!--       Determine the global parent-grid index bounds       
1501          parent_bound_global(1) = MINVAL( parent_bound_all(1,:) )
1502          parent_bound_global(2) = MAXVAL( parent_bound_all(2,:) )
1503          parent_bound_global(3) = MINVAL( parent_bound_all(3,:) )
1504          parent_bound_global(4) = MAXVAL( parent_bound_all(4,:) )
1505       ENDIF
1506!
1507!--    Broadcast the global parent-grid index bounds to all current child processes
1508       CALL MPI_BCAST( parent_bound_global, 4, MPI_INTEGER, 0, comm2d, ierr )
1509       iplg = parent_bound_global(1)
1510       iprg = parent_bound_global(2)
1511       jpsg = parent_bound_global(3)
1512       jpng = parent_bound_global(4)
1513       WRITE( 9, "('pmci_map_child_grid_to_parent_grid. Global parent-grid index bounds iplg, iprg, jpsg, jpng: ',4(i4,2x))" ) &
1514            iplg, iprg, jpsg, jpng
1515       FLUSH( 9 )
1516       
1517    END SUBROUTINE pmci_map_child_grid_to_parent_grid
1518
1519     
1520     
1521    SUBROUTINE pmci_define_index_mapping
1522!
1523!--    Precomputation of the mapping of the child- and parent-grid indices.
1524
1525       IMPLICIT NONE
1526
1527       INTEGER(iwp) ::  i         !< Child-grid index in the x-direction
1528       INTEGER(iwp) ::  ii        !< Parent-grid index in the x-direction
1529       INTEGER(iwp) ::  istart    !<
1530       INTEGER(iwp) ::  ir        !<
1531       INTEGER(iwp) ::  iw        !< Child-grid index limited to -1 <= iw <= nx+1 for wall_flags_static_0
1532       INTEGER(iwp) ::  j         !< Child-grid index in the y-direction
1533       INTEGER(iwp) ::  jj        !< Parent-grid index in the y-direction
1534       INTEGER(iwp) ::  jstart    !<
1535       INTEGER(iwp) ::  jr        !<
1536       INTEGER(iwp) ::  jw        !< Child-grid index limited to -1 <= jw <= ny+1 for wall_flags_static_0
1537       INTEGER(iwp) ::  k         !< Child-grid index in the z-direction
1538       INTEGER(iwp) ::  kk        !< Parent-grid index in the z-direction
1539       INTEGER(iwp) ::  kstart    !<
1540       INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1 for wall_flags_static_0
1541
1542       REAL(wp)     ::  tolex     !< Tolerance for grid-line matching in x-direction   
1543       REAL(wp)     ::  toley     !< Tolerance for grid-line matching in y-direction   
1544       REAL(wp)     ::  tolez     !< Tolerance for grid-line matching in z-direction   
1545
1546!
1547!--    Grid-line tolerances.
1548       tolex = tolefac * dx
1549       toley = tolefac * dy
1550       tolez = tolefac * dz(1)
1551!
1552!--    Allocate child-grid work arrays for interpolation.
1553       igsr = NINT( pg%dx / dx, iwp )
1554       jgsr = NINT( pg%dy / dy, iwp )
1555       kgsr = NINT( pg%dzw(1) / dzw(1), iwp )
1556       WRITE(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
1557       FLUSH(9)
1558!       
1559!--    Determine index bounds for the parent-grid work arrays for
1560!--    interpolation and allocate them.
1561       CALL pmci_allocate_workarrays
1562!       
1563!--    Define the MPI-datatypes for parent-grid work array
1564!--    exchange between the PE-subdomains.
1565       CALL pmci_create_workarray_exchange_datatypes
1566!
1567!--    First determine kcto and kctw which refer to the uppermost
1568!--    parent-grid levels below the child top-boundary level.
1569!--    Note that these comparison tests are not round-off-error
1570!--    sensitive and therefore tolerance buffering is not needed here.
1571       kk = 0
1572       DO WHILE ( pg%zu(kk) <= zu(nzt) )
1573          kk = kk + 1
1574       ENDDO
1575       kcto = kk - 1
1576
1577       kk = 0
1578       DO WHILE ( pg%zw(kk) <= zw(nzt-1) )
1579          kk = kk + 1
1580       ENDDO
1581       kctw = kk - 1
1582
1583       WRITE( 9, "('kcto, kctw = ', 2(i3,2x))" ) kcto, kctw
1584       FLUSH( 9 )
1585!       
1586!--    In case of two-way coupling, check that the child domain is sufficiently
1587!--    large in terms of the number of parent-grid cells covered. Otherwise
1588!--    anterpolation is not possible.
1589       IF ( nesting_mode == 'two-way')  THEN
1590          CALL pmci_check_child_domain_size
1591       ENDIF
1592       
1593       ALLOCATE( iflu(ipla:ipra) )
1594       ALLOCATE( iflo(ipla:ipra) )
1595       ALLOCATE( ifuu(ipla:ipra) )
1596       ALLOCATE( ifuo(ipla:ipra) )
1597       ALLOCATE( jflv(jpsa:jpna) )
1598       ALLOCATE( jflo(jpsa:jpna) )
1599       ALLOCATE( jfuv(jpsa:jpna) )
1600       ALLOCATE( jfuo(jpsa:jpna) )       
1601       ALLOCATE( kflw(0:pg%nz+1) )
1602       ALLOCATE( kflo(0:pg%nz+1) )
1603       ALLOCATE( kfuw(0:pg%nz+1) )
1604       ALLOCATE( kfuo(0:pg%nz+1) )
1605       ALLOCATE( ijkfc_u(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1606       ALLOCATE( ijkfc_v(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1607       ALLOCATE( ijkfc_w(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1608       ALLOCATE( ijkfc_s(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1609
1610       ijkfc_u = 0
1611       ijkfc_v = 0
1612       ijkfc_w = 0
1613       ijkfc_s = 0
1614!
1615!--    i-indices of u for each ii-index value
1616       istart = nxlg
1617       DO  ii = ipla, ipra
1618!
1619!--       The parent and child grid lines do always match in x, hence we
1620!--       use only the local k,j-child-grid plane for the anterpolation.
1621!--       However, icru still has to be stored separately as these index bounds
1622!--       are passed as arguments to the interpolation and anterpolation
1623!--       subroutines.
1624!--       Note that this comparison test is round-off-error sensitive
1625!--       and therefore tolerance buffering is needed here.
1626          i = istart
1627          DO WHILE ( pg%coord_x(ii) - coord_x(i) > tolex  .AND. i < nxrg )
1628             i = i + 1
1629          ENDDO
1630          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
1631          ifuu(ii) = iflu(ii)
1632          istart   = iflu(ii)
1633!
1634!--       Print out the index bounds for checking and debugging purposes
1635          WRITE( 9, "('pmci_define_index_mapping, ii, iflu, ifuu: ', 3(i4,2x))" )                   &
1636               ii, iflu(ii), ifuu(ii)
1637          FLUSH( 9 )
1638       ENDDO
1639       WRITE( 9, * )
1640!
1641!--    i-indices of others for each ii-index value.
1642!--    Note that these comparison tests are not round-off-error
1643!--    sensitive and therefore tolerance buffering is not needed here.
1644       istart = nxlg
1645       DO  ii = ipla, ipra
1646          i = istart
1647          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.  ( i < nxrg ) )
1648             i  = i + 1
1649          ENDDO
1650          iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
1651          ir = i
1652          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx < pg%coord_x(ii) + pg%dx )  .AND.  ( i < nxrg+1 ) )
1653             i  = i + 1
1654             ir = MIN( i, nxrg )
1655          ENDDO
1656          ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg )
1657          istart = iflo(ii)
1658!
1659!--       Print out the index bounds for checking and debugging purposes
1660          WRITE( 9, "('pmci_define_index_mapping, ii, iflo, ifuo: ', 3(i4,2x))" )                   &
1661               ii, iflo(ii), ifuo(ii)
1662          FLUSH( 9 )
1663       ENDDO
1664       WRITE( 9, * )
1665!
1666!--    j-indices of v for each jj-index value
1667       jstart = nysg
1668       DO  jj = jpsa, jpna
1669!
1670!--       The parent and child grid lines do always match in y, hence we
1671!--       use only the local k,i-child-grid plane for the anterpolation.
1672!--       However, jcnv still has to be stored separately as these index bounds
1673!--       are passed as arguments to the interpolation and anterpolation
1674!--       subroutines.
1675!--       Note that this comparison test is round-off-error sensitive
1676!--       and therefore tolerance buffering is needed here.
1677          j = jstart
1678          DO WHILE ( pg%coord_y(jj) - coord_y(j) > toley  .AND. j < nyng )
1679             j = j + 1
1680          ENDDO
1681          jflv(jj) = MIN( MAX( j, nysg ), nyng )
1682          jfuv(jj) = jflv(jj)
1683          jstart   = jflv(jj)
1684!
1685!--       Print out the index bounds for checking and debugging purposes
1686          WRITE( 9, "('pmci_define_index_mapping, jj, jflv, jfuv: ', 3(i4,2x))" )                   &
1687               jj, jflv(jj), jfuv(jj)
1688          FLUSH(9)
1689       ENDDO
1690       WRITE( 9, * )
1691!
1692!--    j-indices of others for each jj-index value
1693!--    Note that these comparison tests are not round-off-error
1694!--    sensitive and therefore tolerance buffering is not needed here.
1695       jstart = nysg
1696       DO  jj = jpsa, jpna
1697          j = jstart
1698          DO WHILE ( ( coord_y(j) + 0.5_wp * dy < pg%coord_y(jj) ) .AND. ( j < nyng ) )
1699             j  = j + 1
1700          ENDDO
1701          jflo(jj) = MIN( MAX( j, nysg ), nyng )
1702          jr = j
1703          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy < pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
1704             j  = j + 1
1705             jr = MIN( j, nyng )
1706          ENDDO
1707          jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng )
1708          jstart = jflo(jj)
1709!
1710!--       Print out the index bounds for checking and debugging purposes
1711          WRITE( 9, "('pmci_define_index_mapping, jj, jflo, jfuo: ', 3(i4,2x))" )                   &
1712               jj, jflo(jj), jfuo(jj)
1713          FLUSH( 9 )
1714       ENDDO
1715       WRITE( 9, * )
1716!
1717!--    k-indices of w for each kk-index value
1718!--    Note that anterpolation index limits are needed also for the top boundary
1719!--    ghost cell level because they are used also in the interpolation.
1720       kstart  = 0
1721       kflw(0) = 0
1722       kfuw(0) = 0
1723       DO  kk = 1, pg%nz+1
1724!
1725!--       The parent and child grid lines do always match in z, hence we
1726!--       use only the local j,i-child-grid plane for the anterpolation.
1727!--       However, kctw still has to be stored separately as these index bounds
1728!--       are passed as arguments to the interpolation and anterpolation
1729!--       subroutines.
1730!--       Note that this comparison test is round-off-error sensitive
1731!--       and therefore tolerance buffering is needed here.
1732          k = kstart
1733          DO WHILE ( ( pg%zw(kk) - zw(k) > tolez )  .AND.  ( k < nzt+1 ) )
1734             k = k + 1
1735          ENDDO
1736          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
1737          kfuw(kk) = kflw(kk)
1738          kstart   = kflw(kk)
1739!
1740!--       Print out the index bounds for checking and debugging purposes
1741          WRITE( 9, "('pmci_define_index_mapping, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))" )      &
1742               kk, kflw(kk), kfuw(kk), nzt,  pg%zu(kk), pg%zw(kk)
1743          FLUSH( 9 )
1744       ENDDO
1745       WRITE( 9, * )
1746!
1747!--    k-indices of others for each kk-index value
1748       kstart  = 0
1749       kflo(0) = 0
1750       kfuo(0) = 0
1751!
1752!--    Note that anterpolation index limits are needed also for the top boundary
1753!--    ghost cell level because they are used also in the interpolation.
1754!--    Note that these comparison tests are not round-off-error
1755!--    sensitive and therefore tolerance buffering is not needed here.
1756       DO  kk = 1, pg%nz+1
1757          k = kstart
1758          DO WHILE ( ( zu(k) < pg%zw(kk-1) )  .AND.  ( k <= nzt ) )
1759             k = k + 1
1760          ENDDO
1761          kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
1762          DO WHILE ( ( zu(k) < pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
1763             k = k + 1
1764             IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
1765          ENDDO
1766          kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 )
1767          kstart = kflo(kk)
1768       ENDDO
1769!
1770!--    Print out the index bounds for checking and debugging purposes
1771       DO  kk = 1, pg%nz+1
1772          WRITE( 9, "('pmci_define_index_mapping, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))" )      &
1773               kk, kflo(kk), kfuo(kk), nzt,  pg%zu(kk), pg%zw(kk)
1774          FLUSH( 9 )
1775       ENDDO
1776       WRITE( 9, * )
1777!
1778!--    Precomputation of number of child-grid nodes inside parent-grid cells.
1779!--    Note that ii, jj, and kk are parent-grid indices.
1780!--    This information is needed in the anterpolation.
1781!--    The indices for wall_flags_static_0 (kw,jw,iw) must be limited to the range
1782!--    [-1,...,nx/ny/nzt+1] in order to avoid zero values on the outer ghost nodes.
1783       DO  ii = ipla, ipra
1784          DO  jj = jpsa, jpna
1785             DO  kk = 0, pg%nz+1
1786!
1787!--             u-component
1788                DO  i = iflu(ii), ifuu(ii)
1789                   iw = MAX( MIN( i, nx+1 ), -1 )
1790                   DO  j = jflo(jj), jfuo(jj)
1791                      jw = MAX( MIN( j, ny+1 ), -1 )
1792                      DO  k = kflo(kk), kfuo(kk)
1793                         kw = MIN( k, nzt+1 )               
1794                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii)                                      &
1795                              + MERGE( 1, 0, BTEST( wall_flags_static_0(kw,jw,iw), 1 ) )
1796                      ENDDO
1797                   ENDDO
1798                ENDDO
1799!
1800!--             v-component
1801                DO  i = iflo(ii), ifuo(ii)
1802                   iw = MAX( MIN( i, nx+1 ), -1 )
1803                   DO  j = jflv(jj), jfuv(jj)
1804                      jw = MAX( MIN( j, ny+1 ), -1 )
1805                      DO  k = kflo(kk), kfuo(kk)
1806                         kw = MIN( k, nzt+1 )                                       
1807                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii)                                      &
1808                              + MERGE( 1, 0, BTEST( wall_flags_static_0(kw,jw,iw), 2 ) )
1809                      ENDDO
1810                   ENDDO
1811                ENDDO
1812!
1813!--             scalars
1814                DO  i = iflo(ii), ifuo(ii)
1815                   iw = MAX( MIN( i, nx+1 ), -1 )
1816                   DO  j = jflo(jj), jfuo(jj)
1817                      jw = MAX( MIN( j, ny+1 ), -1 )
1818                      DO  k = kflo(kk), kfuo(kk)
1819                         kw = MIN( k, nzt+1 )
1820                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii)                                      &
1821                              + MERGE( 1, 0, BTEST( wall_flags_static_0(kw,jw,iw), 0 ) )
1822                      ENDDO
1823                   ENDDO
1824                ENDDO
1825!
1826!--             w-component
1827                DO  i = iflo(ii), ifuo(ii)
1828                   iw = MAX( MIN( i, nx+1 ), -1 )
1829                   DO  j = jflo(jj), jfuo(jj)
1830                      jw = MAX( MIN( j, ny+1 ), -1 )
1831                      DO  k = kflw(kk), kfuw(kk)
1832                         kw = MIN( k, nzt+1 )
1833                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii)                                      &
1834                              + MERGE( 1, 0, BTEST( wall_flags_static_0(kw,jw,iw), 3 ) )
1835                      ENDDO
1836                   ENDDO
1837                ENDDO
1838
1839             ENDDO  ! kk       
1840          ENDDO  ! jj
1841       ENDDO  ! ii
1842
1843    END SUBROUTINE pmci_define_index_mapping
1844
1845
1846
1847    SUBROUTINE pmci_check_child_domain_size
1848!       
1849!--    Check if the child domain is too small in terms of number of parent-grid cells
1850!--    covered so that anterpolation buffers fill the whole domain so that anterpolation
1851!--    not possible. Also, check that anterpolation_buffer_width is not too large to 
1852!--    prevent anterpolation.
1853       IMPLICIT NONE
1854     
1855!
1856!--    First x-direction
1857       IF ( iplg + 3 + anterpolation_buffer_width > iprg - 3 - anterpolation_buffer_width )  THEN
1858          IF ( iprg - iplg + 1 < 7 )  THEN
1859!
1860!--          Error
1861             WRITE( message_string, * ) 'child domain too narrow for anterpolation in x-direction'
1862             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
1863          ELSE IF ( iprg - iplg + 1 < 11 )  THEN
1864!               
1865!--          Warning
1866             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
1867             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
1868             anterpolation_buffer_width = 0
1869          ELSE
1870!               
1871!--          Informative message
1872             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
1873             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
1874             anterpolation_buffer_width = 2
1875          ENDIF
1876       ENDIF
1877!
1878!--    Then y-direction         
1879       IF ( jpsg + 3 + anterpolation_buffer_width > jpng - 3 - anterpolation_buffer_width )  THEN
1880          IF ( jpng - jpsg + 1 < 7 )  THEN
1881!
1882!--          Error
1883             WRITE( message_string, * ) 'child domain too narrow for anterpolation in y-direction'
1884             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
1885          ELSE IF ( jpng - jpsg + 1 < 11 )  THEN
1886!               
1887!--          Warning
1888             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
1889             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
1890             anterpolation_buffer_width = 0
1891          ELSE
1892!               
1893!--          Informative message
1894             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
1895             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
1896             anterpolation_buffer_width = 2
1897          ENDIF
1898       ENDIF
1899!
1900!--    Finally z-direction               
1901       IF ( kctw - 1 - anterpolation_buffer_width < 1 )  THEN
1902          IF ( kctw - 1 < 1 )  THEN
1903!
1904!--          Error
1905             WRITE( message_string, * ) 'child domain too shallow for anterpolation in z-direction'
1906             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
1907          ELSE IF ( kctw - 3 < 1 )  THEN
1908!               
1909!--          Warning
1910             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
1911             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
1912             anterpolation_buffer_width = 0
1913          ELSE
1914!               
1915!--          Informative message
1916             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
1917             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
1918             anterpolation_buffer_width = 2 
1919          ENDIF
1920       ENDIF
1921
1922    END SUBROUTINE pmci_check_child_domain_size
1923
1924
1925   
1926    SUBROUTINE pmci_allocate_workarrays
1927!
1928!--    Allocate parent-grid work-arrays for interpolation
1929       IMPLICIT NONE
1930
1931!
1932!--    Determine and store the PE-subdomain dependent index bounds
1933       IF ( bc_dirichlet_l )  THEN
1934          iplw = ipl + 1
1935       ELSE
1936          iplw = ipl - 1
1937       ENDIF
1938
1939       IF ( bc_dirichlet_r )  THEN
1940          iprw = ipr - 1
1941       ELSE
1942          iprw = ipr + 1
1943       ENDIF
1944
1945       IF ( bc_dirichlet_s )  THEN
1946          jpsw = jps + 1
1947       ELSE
1948          jpsw = jps - 1
1949       ENDIF
1950
1951       IF ( bc_dirichlet_n )  THEN
1952          jpnw = jpn - 1
1953       ELSE
1954          jpnw = jpn + 1
1955       ENDIF
1956!
1957!--    Left and right boundaries.
1958       ALLOCATE( workarr_lr(0:pg%nz+1,jpsw:jpnw,0:2) )
1959!
1960!--    South and north boundaries.
1961       ALLOCATE( workarr_sn(0:pg%nz+1,0:2,iplw:iprw) )
1962!
1963!--    Top boundary.
1964       ALLOCATE( workarr_t(0:2,jpsw:jpnw,iplw:iprw) )
1965
1966    END SUBROUTINE pmci_allocate_workarrays
1967
1968
1969
1970    SUBROUTINE pmci_create_workarray_exchange_datatypes
1971!
1972!--    Define specific MPI types for workarr-exchange.
1973       IMPLICIT NONE
1974
1975!
1976!--    For the left and right boundaries
1977       CALL MPI_TYPE_VECTOR( 3, pg%nz+2, (jpnw-jpsw+1)*(pg%nz+2), MPI_REAL,                         &
1978            workarr_lr_exchange_type, ierr )
1979       CALL MPI_TYPE_COMMIT( workarr_lr_exchange_type, ierr )
1980!
1981!--    For the south and north boundaries
1982       CALL MPI_TYPE_VECTOR( 1, 3*(pg%nz+2), 3*(pg%nz+2), MPI_REAL,                                 &
1983            workarr_sn_exchange_type, ierr )
1984       CALL MPI_TYPE_COMMIT( workarr_sn_exchange_type, ierr )
1985!
1986!--    For the top-boundary x-slices
1987       CALL MPI_TYPE_VECTOR( iprw-iplw+1, 3, 3*(jpnw-jpsw+1), MPI_REAL,                             &
1988            workarr_t_exchange_type_x, ierr )
1989       CALL MPI_TYPE_COMMIT( workarr_t_exchange_type_x, ierr )
1990!
1991!--    For the top-boundary y-slices
1992       CALL MPI_TYPE_VECTOR( 1, 3*(jpnw-jpsw+1), 3*(jpnw-jpsw+1), MPI_REAL,                         &
1993            workarr_t_exchange_type_y, ierr )
1994       CALL MPI_TYPE_COMMIT( workarr_t_exchange_type_y, ierr )
1995       
1996    END SUBROUTINE pmci_create_workarray_exchange_datatypes
1997
1998
1999
2000    SUBROUTINE pmci_check_grid_matching
2001!
2002!--    Check that the grid lines of child and parent do match.
2003!--    Also check that the child subdomain width is not smaller than
2004!--    the parent grid spacing in the respective direction.
2005       IMPLICIT NONE
2006
2007       INTEGER(iwp) ::  non_matching_height = 0              !< Flag for non-matching child-domain height
2008       INTEGER(iwp) ::  non_matching_lower_left_corner = 0   !< Flag for non-matching lower left corner
2009       INTEGER(iwp) ::  non_matching_upper_right_corner = 0  !< Flag for non-matching upper right corner
2010       INTEGER(iwp) ::  non_int_gsr_x = 0                    !< Flag for non-integer grid-spacing ration in x-direction
2011       INTEGER(iwp) ::  non_int_gsr_y = 0                    !< Flag for non-integer grid-spacing ration in y-direction
2012       INTEGER(iwp) ::  non_int_gsr_z = 0                    !< Flag for non-integer grid-spacing ration in z-direction
2013       INTEGER(iwp) ::  too_narrow_pesd_x = 0                !< Flag for too narrow pe-subdomain in x-direction
2014       INTEGER(iwp) ::  too_narrow_pesd_y = 0                !< Flag for too narrow pe-subdomain in y-direction
2015                                                                                                                 
2016       REAL(wp) ::  child_ngp_x_l                            !< Number of gridpoints in child subdomain in x-direction
2017                                                             !< converted to REAL(wp)
2018       REAL(wp) ::  child_ngp_y_l                            !< Number of gridpoints in child subdomain in y-direction
2019                                                             !< converted to REAL(wp)
2020       REAL(wp) ::  gridline_mismatch_x                      !< Mismatch between the parent and child gridlines in the x-direction
2021       REAL(wp) ::  gridline_mismatch_y                      !< Mismatch between the parent and child gridlines in the y-direction
2022       REAL(wp) ::  gsr_mismatch_x                           !< Deviation of the grid-spacing ratio from the nearest integer value, the x-direction
2023       REAL(wp) ::  gsr_mismatch_y                           !< Deviation of the grid-spacing ratio from the nearest integer value, the y-direction
2024       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
2025       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
2026       REAL(wp) ::  tolex                                    !< Tolerance for grid-line matching in x-direction
2027       REAL(wp) ::  toley                                    !< Tolerance for grid-line matching in y-direction
2028       REAL(wp) ::  tolez                                    !< Tolerance for grid-line matching in z-direction
2029
2030       
2031       IF ( myid == 0 )  THEN
2032
2033          tolex = tolefac * dx
2034          toley = tolefac * dy
2035          tolez = tolefac * dz(1)
2036!
2037!--       First check that the child domain lower left corner matches the parent grid lines.
2038          gridline_mismatch_x = ABS( NINT( lower_left_coord_x / pg%dx ) * pg%dx - lower_left_coord_x )
2039          gridline_mismatch_y = ABS( NINT( lower_left_coord_y / pg%dy ) * pg%dy - lower_left_coord_y )
2040          IF ( gridline_mismatch_x > tolex ) non_matching_lower_left_corner = 1
2041          IF ( gridline_mismatch_y > toley ) non_matching_lower_left_corner = 1
2042!
2043!--       Then check that the child doman upper right corner matches the parent grid lines.
2044          upper_right_coord_x = lower_left_coord_x + ( nx + 1 ) * dx
2045          upper_right_coord_y = lower_left_coord_y + ( ny + 1 ) * dy
2046          gridline_mismatch_x = ABS( NINT( upper_right_coord_x / pg%dx ) * pg%dx - upper_right_coord_x )
2047          gridline_mismatch_y = ABS( NINT( upper_right_coord_y / pg%dy ) * pg%dy - upper_right_coord_y )
2048          IF ( gridline_mismatch_x > tolex ) non_matching_upper_right_corner = 1
2049          IF ( gridline_mismatch_y > toley ) non_matching_upper_right_corner = 1
2050!
2051!--       Also check that the cild domain height matches the parent grid lines.
2052          IF ( MOD( zw(nzt), pg%dz ) > tolez ) non_matching_height = 1
2053!
2054!--       Check that the grid-spacing ratios in each direction are integer valued.   
2055          gsr_mismatch_x = ABS( NINT( pg%dx / dx ) * dx - pg%dx )
2056          gsr_mismatch_y = ABS( NINT( pg%dy / dy ) * dy - pg%dy )
2057          IF ( gsr_mismatch_x > tolex )  non_int_gsr_x = 1
2058          IF ( gsr_mismatch_y > toley )  non_int_gsr_y = 1
2059!
2060!--       In the z-direction, all levels need to be checked separately against grid stretching 
2061!--       which is not allowed.
2062          DO  n = 0, kctw+1
2063             IF ( ABS( pg%zw(n) - zw(kflw(n)) ) > tolez )  non_int_gsr_z = 1
2064          ENDDO
2065
2066          child_ngp_x_l = REAL( nxr - nxl + 1, KIND=wp )
2067          IF ( child_ngp_x_l / REAL( igsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_x = 1
2068          child_ngp_y_l = REAL( nyn - nys + 1, KIND=wp )
2069          IF ( child_ngp_y_l / REAL( jgsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_y = 1
2070         
2071          IF ( non_matching_height > 0 )  THEN
2072             WRITE( message_string, * ) 'nested child domain height must match ',                   &
2073                                        'its parent grid lines'
2074             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2075          ENDIF
2076
2077          IF ( non_matching_lower_left_corner > 0 )  THEN
2078             WRITE( message_string, * ) 'nested child domain lower left ',                          &
2079                                        'corner must match its parent grid lines'
2080             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2081          ENDIF
2082
2083          IF ( non_matching_upper_right_corner > 0 )  THEN
2084             WRITE( message_string, * ) 'nested child domain upper right ',                         &
2085                                        'corner must match its parent grid lines'
2086             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2087          ENDIF
2088
2089          IF ( non_int_gsr_x > 0 )  THEN
2090             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dx / child dx ) ',     &
2091                                        'must have an integer value'
2092             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2093          ENDIF
2094
2095          IF ( non_int_gsr_y > 0 )  THEN
2096             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dy / child dy ) ',     &
2097                                        'must have an integer value'
2098             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2099          ENDIF
2100
2101          IF ( non_int_gsr_z > 0 )  THEN
2102             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dz / child dz ) ',     &
2103                                        'must have an integer value for each z-level'
2104             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2105          ENDIF
2106
2107          IF ( too_narrow_pesd_x > 0 )  THEN
2108            WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
2109                                        'smaller than its parent grid dx. Change the PE-grid ',     &
2110                                        'setting (npex, npey) to satisfy this requirement.' 
2111             CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
2112          ENDIF
2113 
2114          IF ( too_narrow_pesd_y > 0 )  THEN
2115             WRITE( message_string, * ) 'child subdomain width in y-direction must not be ',        &
2116                                        'smaller than its parent grid dy. Change the PE-grid ',     &
2117                                        'setting (npex, npey) to satisfy this requirement.' 
2118             CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
2119          ENDIF
2120                 
2121       ENDIF  !  ( myid == 0 )
2122       
2123    END SUBROUTINE pmci_check_grid_matching
2124
2125
2126
2127    SUBROUTINE pmci_compute_face_areas
2128
2129       IMPLICIT NONE
2130       REAL(wp)  :: face_area_local   !< Local (for the current pe) integral face area of the left boundary
2131       REAL(wp)  :: sub_sum           !< Intermediate sum in order to improve the accuracy of the summation
2132
2133       INTEGER(iwp) :: i              !< Running index in the x-direction
2134       INTEGER(iwp) :: j              !< Running index in the y-direction
2135       INTEGER(iwp) :: k              !< Running index in the z-direction
2136       INTEGER(iwp) :: k_wall         !< Local topography top k-index
2137       INTEGER(iwp) :: n              !< Running index over boundary faces
2138
2139       
2140!
2141!--    Sum up the volume flow through the left boundary
2142       face_area(1) = 0.0_wp
2143       face_area_local = 0.0_wp
2144       IF ( bc_dirichlet_l )  THEN
2145          i = 0
2146          DO  j = nys, nyn
2147             sub_sum = 0.0_wp
2148             k_wall = topo_top_ind(j,i,1)
2149             DO   k = k_wall + 1, nzt
2150                sub_sum = sub_sum + dzw(k)
2151             ENDDO
2152             face_area_local =  face_area_local + dy * sub_sum
2153          ENDDO
2154       ENDIF
2155       
2156#if defined( __parallel )
2157       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2158       CALL MPI_ALLREDUCE( face_area_local, face_area(1), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2159#else
2160       face_area(1) = face_area_local
2161#endif
2162!
2163!--    Sum up the volume flow through the right boundary
2164       face_area(2) = 0.0_wp
2165       face_area_local = 0.0_wp
2166       IF ( bc_dirichlet_r )  THEN
2167          i = nx
2168          DO  j = nys, nyn
2169             sub_sum = 0.0_wp
2170             k_wall = topo_top_ind(j,i,1)
2171             DO   k = k_wall + 1, nzt
2172                sub_sum = sub_sum + dzw(k)
2173             ENDDO
2174             face_area_local =  face_area_local + dy * sub_sum
2175          ENDDO
2176       ENDIF
2177       
2178#if defined( __parallel )
2179       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2180       CALL MPI_ALLREDUCE( face_area_local, face_area(2), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2181#else
2182       face_area(2) = face_area_local
2183#endif
2184!
2185!--    Sum up the volume flow through the south boundary
2186       face_area(3) = 0.0_wp
2187       face_area_local = 0.0_wp
2188       IF ( bc_dirichlet_s )  THEN
2189          j = 0
2190          DO  i = nxl, nxr
2191             sub_sum = 0.0_wp
2192             k_wall = topo_top_ind(j,i,2)
2193             DO  k = k_wall + 1, nzt
2194                sub_sum = sub_sum + dzw(k)
2195             ENDDO
2196             face_area_local = face_area_local + dx * sub_sum
2197          ENDDO
2198       ENDIF
2199       
2200#if defined( __parallel )
2201       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2202       CALL MPI_ALLREDUCE( face_area_local, face_area(3), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2203#else
2204       face_area(3) = face_area_local
2205#endif
2206!
2207!--    Sum up the volume flow through the north boundary
2208       face_area(4) = 0.0_wp
2209       face_area_local = 0.0_wp
2210       IF ( bc_dirichlet_n )  THEN
2211          j = ny
2212          DO  i = nxl, nxr
2213             sub_sum = 0.0_wp
2214             k_wall = topo_top_ind(j,i,2)
2215             DO  k = k_wall + 1, nzt
2216                sub_sum = sub_sum + dzw(k)
2217             ENDDO
2218             face_area_local = face_area_local + dx * sub_sum
2219          ENDDO
2220       ENDIF
2221       
2222#if defined( __parallel )
2223       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2224       CALL MPI_ALLREDUCE( face_area_local, face_area(4), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2225#else
2226       face_area(4) = face_area_local
2227#endif
2228!
2229!--    The top face area does not depend on the topography at all.       
2230       face_area(5) = ( nx + 1 ) * ( ny + 1 ) * dx * dy
2231!
2232!--    The 6th element is used for the total area
2233       face_area(6) = 0.0_wp
2234       DO  n = 1, 5
2235          face_area(6) = face_area(6) + face_area(n)
2236       ENDDO
2237
2238!       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
2239!       flush( 9 )
2240       
2241    END SUBROUTINE pmci_compute_face_areas
2242#endif
2243   
2244 END SUBROUTINE pmci_setup_child
2245
2246
2247
2248 SUBROUTINE pmci_setup_coordinates
2249
2250#if defined( __parallel )
2251    IMPLICIT NONE
2252
2253    INTEGER(iwp) ::  i   !<
2254    INTEGER(iwp) ::  j   !<
2255
2256!
2257!-- Create coordinate arrays.
2258    ALLOCATE( coord_x(-nbgp:nx+nbgp) )
2259    ALLOCATE( coord_y(-nbgp:ny+nbgp) )
2260     
2261    DO  i = -nbgp, nx + nbgp
2262       coord_x(i) = lower_left_coord_x + i * dx
2263    ENDDO
2264
2265    DO  j = -nbgp, ny + nbgp
2266       coord_y(j) = lower_left_coord_y + j * dy
2267    ENDDO
2268
2269#endif
2270 END SUBROUTINE pmci_setup_coordinates
2271
2272!------------------------------------------------------------------------------!
2273! Description:
2274! ------------
2275!> In this subroutine the number of coupled arrays is determined.
2276!------------------------------------------------------------------------------!
2277  SUBROUTINE pmci_num_arrays 
2278               
2279#if defined( __parallel )
2280    USE pmc_general,                                                           &
2281        ONLY:  pmc_max_array
2282
2283    IMPLICIT NONE
2284!
2285!-- The number of coupled arrays depends on the model settings. At least
2286!-- 5 arrays need to be coupled (u, v, w, e, diss).  Please note, actually
2287!-- e and diss (TKE and dissipation rate) are only required if RANS-RANS
2288!-- nesting is applied, but memory is allocated nevertheless. This is because
2289!-- the information whether they are needed or not is retrieved at a later
2290!-- point in time. In case e and diss are not needed, they are also not
2291!-- exchanged between parent and child.
2292    pmc_max_array = 5
2293!
2294!-- pt
2295    IF ( .NOT. neutral )  pmc_max_array = pmc_max_array + 1
2296   
2297    IF ( humidity )  THEN
2298!
2299!--    q
2300       pmc_max_array = pmc_max_array + 1
2301!
2302!--    qc, nc
2303       IF ( bulk_cloud_model  .AND.  microphysics_morrison )                   &
2304          pmc_max_array = pmc_max_array + 2
2305!
2306!--    qr, nr
2307       IF ( bulk_cloud_model  .AND.  microphysics_seifert )                    &
2308          pmc_max_array = pmc_max_array + 2
2309    ENDIF
2310!
2311!-- s
2312    IF ( passive_scalar )  pmc_max_array = pmc_max_array + 1
2313!
2314!-- nr_part, part_adr
2315    IF ( particle_advection )  pmc_max_array = pmc_max_array + 2
2316!
2317!-- Chemistry, depends on number of species
2318    IF ( air_chemistry  .AND.  nesting_chem )  pmc_max_array = pmc_max_array + nspec
2319!
2320!-- SALSA, depens on the number aerosol size bins and chemical components +
2321!-- the number of default gases
2322    IF ( salsa  .AND.  nesting_salsa )  pmc_max_array = pmc_max_array + nbins_aerosol +            &
2323                                                        nbins_aerosol * ncomponents_mass
2324    IF ( .NOT. salsa_gases_from_chem )  pmc_max_array = pmc_max_array + ngases_salsa
2325
2326#endif
2327   
2328 END SUBROUTINE pmci_num_arrays
2329
2330
2331 SUBROUTINE pmci_set_array_pointer( name, child_id, nz_child, n )
2332   
2333    IMPLICIT NONE
2334   
2335    INTEGER(iwp), INTENT(IN) ::  child_id  !<
2336    INTEGER(iwp), INTENT(IN) ::  nz_child  !<
2337   
2338    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n          !< index of chemical species
2339   
2340    CHARACTER(LEN=*), INTENT(IN) ::  name             !<
2341
2342#if defined( __parallel )     
2343!
2344!-- Local variables:       
2345    INTEGER(iwp) ::  ierr                             !< MPI error code
2346
2347    INTEGER(idp), POINTER, DIMENSION(:,:) ::  i_2d    !<
2348       
2349    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d      !<
2350    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d      !<
2351    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d_sec  !<
2352   
2353
2354    NULLIFY( p_3d )
2355    NULLIFY( p_2d )
2356    NULLIFY( i_2d )
2357!
2358!-- List of array names, which can be coupled.
2359!-- In case of 3D please change also the second array for the pointer version
2360    IF ( TRIM(name) == "u"          )  p_3d => u
2361    IF ( TRIM(name) == "v"          )  p_3d => v
2362    IF ( TRIM(name) == "w"          )  p_3d => w
2363    IF ( TRIM(name) == "e"          )  p_3d => e
2364    IF ( TRIM(name) == "pt"         )  p_3d => pt
2365    IF ( TRIM(name) == "q"          )  p_3d => q
2366    IF ( TRIM(name) == "qc"         )  p_3d => qc
2367    IF ( TRIM(name) == "qr"         )  p_3d => qr
2368    IF ( TRIM(name) == "nr"         )  p_3d => nr
2369    IF ( TRIM(name) == "nc"         )  p_3d => nc
2370    IF ( TRIM(name) == "s"          )  p_3d => s
2371    IF ( TRIM(name) == "diss"       )  p_3d => diss   
2372    IF ( TRIM(name) == "nr_part"    )  i_2d => nr_part
2373    IF ( TRIM(name) == "part_adr"   )  i_2d => part_adr
2374    IF ( INDEX( TRIM(name), "chem_" ) /= 0      )  p_3d => chem_species(n)%conc
2375    IF ( INDEX( TRIM(name), "an_" ) /= 0  )  p_3d => aerosol_number(n)%conc
2376    IF ( INDEX( TRIM(name), "am_" ) /= 0 )  p_3d => aerosol_mass(n)%conc
2377    IF ( INDEX( TRIM(name), "sg_" ) /= 0  .AND.  .NOT. salsa_gases_from_chem ) &
2378       p_3d => salsa_gas(n)%conc
2379!
2380!-- Next line is just an example for a 2D array (not active for coupling!)
2381!-- Please note, that z0 has to be declared as TARGET array in modules.f90.
2382!    IF ( TRIM(name) == "z0" )    p_2d => z0
2383    IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
2384    IF ( TRIM(name) == "v"    )  p_3d_sec => v_2
2385    IF ( TRIM(name) == "w"    )  p_3d_sec => w_2
2386    IF ( TRIM(name) == "e"    )  p_3d_sec => e_2
2387    IF ( TRIM(name) == "pt"   )  p_3d_sec => pt_2
2388    IF ( TRIM(name) == "q"    )  p_3d_sec => q_2
2389    IF ( TRIM(name) == "qc"   )  p_3d_sec => qc_2
2390    IF ( TRIM(name) == "qr"   )  p_3d_sec => qr_2
2391    IF ( TRIM(name) == "nr"   )  p_3d_sec => nr_2
2392    IF ( TRIM(name) == "nc"   )  p_3d_sec => nc_2
2393    IF ( TRIM(name) == "s"    )  p_3d_sec => s_2
2394    IF ( TRIM(name) == "diss" )  p_3d_sec => diss_2
2395    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d_sec => spec_conc_2(:,:,:,n)
2396    IF ( INDEX( TRIM(name), "an_" )   /= 0 )  p_3d_sec => nconc_2(:,:,:,n)
2397    IF ( INDEX( TRIM(name), "am_" )   /= 0 )  p_3d_sec => mconc_2(:,:,:,n)
2398    IF ( INDEX( TRIM(name), "sg_" )   /= 0  .AND.  .NOT. salsa_gases_from_chem ) &
2399                                 p_3d_sec => gconc_2(:,:,:,n)
2400
2401    IF ( ASSOCIATED( p_3d ) )  THEN
2402       CALL pmc_s_set_dataarray( child_id, p_3d, nz_child, nz, array_2 = p_3d_sec )
2403    ELSEIF  ( ASSOCIATED( p_2d ) )  THEN
2404       CALL pmc_s_set_dataarray( child_id, p_2d )
2405    ELSEIF  ( ASSOCIATED( i_2d ) )  THEN
2406       CALL pmc_s_set_dataarray( child_id, i_2d )
2407    ELSE
2408!
2409!--    Give only one message for the root domain
2410       IF ( pmc_is_rootmodel()  .AND.  myid == 0 )  THEN
2411          message_string = 'pointer for array "' // TRIM( name ) // '" can''t be associated'
2412          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
2413       ELSE
2414!
2415!--       Avoid others to continue
2416          CALL MPI_BARRIER( comm2d, ierr )
2417       ENDIF
2418       
2419    ENDIF
2420   
2421#endif
2422   
2423 END SUBROUTINE pmci_set_array_pointer
2424
2425
2426     
2427 INTEGER FUNCTION get_number_of_children()
2428
2429    IMPLICIT NONE
2430
2431   
2432#if defined( __parallel )
2433    get_number_of_children = SIZE( pmc_parent_for_child ) - 1
2434#else
2435    get_number_of_children = 0
2436#endif
2437
2438    RETURN
2439
2440 END FUNCTION get_number_of_children
2441
2442
2443 
2444 INTEGER FUNCTION get_childid( id_index )
2445
2446    IMPLICIT NONE
2447
2448    INTEGER, INTENT(IN) ::  id_index   !<
2449
2450   
2451#if defined( __parallel )
2452    get_childid = pmc_parent_for_child(id_index)
2453#else
2454    get_childid = 0
2455#endif
2456
2457    RETURN
2458
2459 END FUNCTION get_childid
2460
2461
2462
2463 SUBROUTINE get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b, sy_coord, sy_coord_b,   &
2464      ny_coord, ny_coord_b, uz_coord, uz_coord_b )
2465   
2466    IMPLICIT NONE
2467
2468    INTEGER,INTENT(IN)   ::  m                     !<
2469
2470    REAL(wp),INTENT(OUT) ::  lx_coord, lx_coord_b  !<
2471    REAL(wp),INTENT(OUT) ::  rx_coord, rx_coord_b  !<
2472    REAL(wp),INTENT(OUT) ::  ny_coord, ny_coord_b  !<
2473    REAL(wp),INTENT(OUT) ::  sy_coord, sy_coord_b  !<
2474    REAL(wp),INTENT(OUT) ::  uz_coord, uz_coord_b  !<
2475
2476   
2477#if defined( __parallel )
2478   
2479    lx_coord = childgrid(m)%lx_coord
2480    rx_coord = childgrid(m)%rx_coord
2481    sy_coord = childgrid(m)%sy_coord
2482    ny_coord = childgrid(m)%ny_coord
2483    uz_coord = childgrid(m)%uz_coord
2484   
2485    lx_coord_b = childgrid(m)%lx_coord_b
2486    rx_coord_b = childgrid(m)%rx_coord_b
2487    sy_coord_b = childgrid(m)%sy_coord_b
2488    ny_coord_b = childgrid(m)%ny_coord_b
2489    uz_coord_b = childgrid(m)%uz_coord_b
2490   
2491#endif
2492   
2493 END SUBROUTINE get_child_edges
2494
2495
2496
2497 SUBROUTINE  get_child_gridspacing( m, dx, dy,dz )
2498
2499    IMPLICIT NONE
2500   
2501    INTEGER, INTENT(IN)             ::  m      !<
2502
2503    REAL(wp), INTENT(OUT)           ::  dx,dy  !<
2504
2505    REAL(wp), INTENT(OUT), OPTIONAL ::  dz     !<
2506
2507
2508#if defined( __parallel )
2509   
2510    dx = childgrid(m)%dx
2511    dy = childgrid(m)%dy
2512    IF ( PRESENT( dz ) )  THEN
2513       dz = childgrid(m)%dz
2514    ENDIF
2515   
2516#endif
2517   
2518 END SUBROUTINE get_child_gridspacing
2519
2520
2521
2522 SUBROUTINE pmci_create_childs_parent_grid_arrays( name, is, ie, js, je, nzc, n  )
2523
2524    IMPLICIT NONE
2525
2526    INTEGER(iwp), INTENT(IN) ::  ie      !<  RENAME ie, is, je, js?
2527    INTEGER(iwp), INTENT(IN) ::  is      !<
2528    INTEGER(iwp), INTENT(IN) ::  je      !<
2529    INTEGER(iwp), INTENT(IN) ::  js      !<
2530    INTEGER(iwp), INTENT(IN) ::  nzc     !<  nzc is pg%nz, but note that pg%nz is not the original nz of parent,
2531                                         !<  but the highest parent-grid level needed for nesting.
2532    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species / salsa variables
2533   
2534    CHARACTER(LEN=*), INTENT(IN) ::  name    !<
2535
2536#if defined( __parallel )
2537!       
2538!-- Local variables:
2539    INTEGER(iwp) ::  ierr    !<
2540   
2541    INTEGER(idp), POINTER,DIMENSION(:,:)   ::  i_2d    !<
2542   
2543    REAL(wp), POINTER,DIMENSION(:,:)       ::  p_2d    !<
2544    REAL(wp), POINTER,DIMENSION(:,:,:)     ::  p_3d    !<
2545   
2546    NULLIFY( p_3d )
2547    NULLIFY( p_2d )
2548    NULLIFY( i_2d )
2549!
2550!-- List of array names, which can be coupled
2551    IF ( TRIM( name ) == "u" )  THEN
2552       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1,js:je,is:ie) )
2553       p_3d => uc
2554    ELSEIF ( TRIM( name ) == "v" )  THEN
2555       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1,js:je,is:ie) )
2556       p_3d => vc
2557    ELSEIF ( TRIM( name ) == "w" )  THEN
2558       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1,js:je,is:ie) )
2559       p_3d => wc
2560    ELSEIF ( TRIM( name ) == "e" )  THEN
2561       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1,js:je,is:ie) )
2562       p_3d => ec
2563    ELSEIF ( TRIM( name ) == "diss" )  THEN
2564       IF ( .NOT. ALLOCATED( dissc ) )  ALLOCATE( dissc(0:nzc+1,js:je,is:ie) )
2565       p_3d => dissc
2566    ELSEIF ( TRIM( name ) == "pt")  THEN
2567       IF ( .NOT. ALLOCATED( ptc ) )  ALLOCATE( ptc(0:nzc+1,js:je,is:ie) )
2568       p_3d => ptc
2569    ELSEIF ( TRIM( name ) == "q")  THEN
2570       IF ( .NOT. ALLOCATED( q_c ) )  ALLOCATE( q_c(0:nzc+1,js:je,is:ie) )
2571       p_3d => q_c
2572    ELSEIF ( TRIM( name ) == "qc")  THEN
2573       IF ( .NOT. ALLOCATED( qcc ) )  ALLOCATE( qcc(0:nzc+1,js:je,is:ie) )
2574       p_3d => qcc
2575    ELSEIF ( TRIM( name ) == "qr")  THEN
2576       IF ( .NOT. ALLOCATED( qrc ) )  ALLOCATE( qrc(0:nzc+1,js:je,is:ie) )
2577       p_3d => qrc
2578    ELSEIF ( TRIM( name ) == "nr")  THEN
2579       IF ( .NOT. ALLOCATED( nrc ) )  ALLOCATE( nrc(0:nzc+1,js:je,is:ie) )
2580       p_3d => nrc
2581    ELSEIF ( TRIM( name ) == "nc")  THEN
2582       IF ( .NOT. ALLOCATED( ncc ) )  ALLOCATE( ncc(0:nzc+1,js:je,is:ie) )
2583       p_3d => ncc
2584    ELSEIF ( TRIM( name ) == "s")  THEN
2585       IF ( .NOT. ALLOCATED( sc ) )  ALLOCATE( sc(0:nzc+1,js:je,is:ie) )
2586       p_3d => sc
2587    ELSEIF ( TRIM( name ) == "nr_part") THEN
2588       IF ( .NOT. ALLOCATED( nr_partc ) )  ALLOCATE( nr_partc(js:je,is:ie) )
2589       i_2d => nr_partc
2590    ELSEIF ( TRIM( name ) == "part_adr") THEN
2591       IF ( .NOT. ALLOCATED( part_adrc ) )  ALLOCATE( part_adrc(js:je,is:ie) )
2592       i_2d => part_adrc
2593    ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
2594       IF ( .NOT. ALLOCATED( chem_spec_c ) ) ALLOCATE( chem_spec_c(0:nzc+1,js:je,is:ie,1:nspec) )
2595       p_3d => chem_spec_c(:,:,:,n)
2596    ELSEIF ( TRIM( name(1:3) ) == "an_" )  THEN
2597       IF ( .NOT. ALLOCATED( aerosol_number_c ) )                              &
2598          ALLOCATE( aerosol_number_c(0:nzc+1,js:je,is:ie,1:nbins_aerosol) )
2599       p_3d => aerosol_number_c(:,:,:,n)
2600    ELSEIF ( TRIM( name(1:3) ) == "am_" )  THEN
2601       IF ( .NOT. ALLOCATED( aerosol_mass_c ) )                                &
2602          ALLOCATE( aerosol_mass_c(0:nzc+1,js:je,is:ie,1:(nbins_aerosol*ncomponents_mass) ) )
2603       p_3d => aerosol_mass_c(:,:,:,n)
2604    ELSEIF ( TRIM( name(1:3) ) == "sg_"  .AND.  .NOT. salsa_gases_from_chem )  &
2605    THEN
2606       IF ( .NOT. ALLOCATED( salsa_gas_c ) )                                   &
2607          ALLOCATE( salsa_gas_c(0:nzc+1,js:je,is:ie,1:ngases_salsa) )
2608       p_3d => salsa_gas_c(:,:,:,n)
2609    !ELSEIF (trim(name) == "z0") then
2610       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
2611       !p_2d => z0c
2612    ENDIF
2613
2614    IF ( ASSOCIATED( p_3d ) )  THEN
2615       CALL pmc_c_set_dataarray( p_3d )
2616    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
2617       CALL pmc_c_set_dataarray( p_2d )
2618    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
2619       CALL pmc_c_set_dataarray( i_2d )
2620    ELSE
2621!
2622!--    Give only one message for the first child domain.
2623       IF ( cpl_id == 2  .AND.  myid == 0 )  THEN
2624          message_string = 'pointer for array "' // TRIM( name ) //            &
2625               '" can''t be associated'
2626          CALL message( 'pmci_create_childs_parent_grid_arrays', 'PA0170', 3, 2, 0, 6, 0 )
2627       ELSE
2628!
2629!--          Prevent others from continuing in case the abort is to come.
2630          CALL MPI_BARRIER( comm2d, ierr )
2631       ENDIF
2632
2633    ENDIF
2634
2635#endif
2636 END SUBROUTINE pmci_create_childs_parent_grid_arrays
2637
2638
2639 SUBROUTINE pmci_parent_initialize
2640
2641!
2642!-- Send data for the children in order to let them create initial
2643!-- conditions by interpolating the parent-domain fields.
2644#if defined( __parallel )
2645    IMPLICIT NONE
2646
2647    INTEGER(iwp) ::  child_id    !<
2648    INTEGER(iwp) ::  m           !<
2649    REAL(wp) ::  waittime        !<
2650
2651
2652    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
2653       child_id = pmc_parent_for_child(m)
2654       CALL pmc_s_fillbuffer( child_id, waittime=waittime )
2655    ENDDO
2656
2657#endif
2658 END SUBROUTINE pmci_parent_initialize
2659
2660
2661
2662 SUBROUTINE pmci_child_initialize
2663
2664!
2665!-- Create initial conditions for the current child domain by interpolating
2666!-- the parent-domain fields.
2667#if defined( __parallel )
2668    IMPLICIT NONE
2669
2670    INTEGER(iwp) ::  ic         !< Child-grid index in x-direction
2671    INTEGER(iwp) ::  jc         !< Child-grid index in y-direction
2672    INTEGER(iwp) ::  kc         !< Child-grid index in z-direction
2673    INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
2674    INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
2675    INTEGER(iwp) ::  lg         !< Running index for salsa gases
2676    INTEGER(iwp) ::  n          !< Running index for chemical species
2677    REAL(wp) ::  waittime       !< Waiting time
2678
2679!
2680!-- Root model is never anyone's child
2681    IF ( .NOT.  pmc_is_rootmodel() )  THEN
2682!
2683!--    Get data from the parent
2684       CALL pmc_c_getbuffer( waittime = waittime )
2685!
2686!--    The interpolation.
2687       CALL pmci_interp_1sto_all ( u, uc, kcto, iflu, ifuu, jflo, jfuo, kflo, kfuo, 'u' )
2688       CALL pmci_interp_1sto_all ( v, vc, kcto, iflo, ifuo, jflv, jfuv, kflo, kfuo, 'v' )
2689       CALL pmci_interp_1sto_all ( w, wc, kctw, iflo, ifuo, jflo, jfuo, kflw, kfuw, 'w' )
2690
2691       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                              &
2692            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                               &
2693               .NOT. constant_diffusion ) )  THEN
2694          CALL pmci_interp_1sto_all ( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 'e' )
2695       ENDIF
2696
2697       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
2698          CALL pmci_interp_1sto_all ( diss, dissc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2699       ENDIF
2700
2701       IF ( .NOT. neutral )  THEN
2702          CALL pmci_interp_1sto_all ( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2703       ENDIF
2704
2705       IF ( humidity )  THEN
2706
2707          CALL pmci_interp_1sto_all ( q, q_c, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2708
2709          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
2710             CALL pmci_interp_1sto_all ( qc, qcc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' ) 
2711             CALL pmci_interp_1sto_all ( nc, ncc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )   
2712          ENDIF
2713
2714          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
2715             CALL pmci_interp_1sto_all ( qr, qrc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2716             CALL pmci_interp_1sto_all ( nr, nrc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2717          ENDIF
2718
2719       ENDIF
2720
2721       IF ( passive_scalar )  THEN
2722          CALL pmci_interp_1sto_all ( s, sc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2723       ENDIF
2724
2725       IF ( air_chemistry  .AND.  nesting_chem )  THEN
2726          DO  n = 1, nspec
2727             CALL pmci_interp_1sto_all ( chem_species(n)%conc, chem_spec_c(:,:,:,n),                &
2728                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2729          ENDDO
2730       ENDIF
2731
2732       IF ( salsa  .AND.  nesting_salsa )  THEN
2733          DO  lb = 1, nbins_aerosol
2734             CALL pmci_interp_1sto_all ( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),       &
2735                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2736          ENDDO
2737          DO  lc = 1, nbins_aerosol * ncomponents_mass
2738             CALL pmci_interp_1sto_all ( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),           &
2739                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2740          ENDDO
2741          IF ( .NOT. salsa_gases_from_chem )  THEN
2742             DO  lg = 1, ngases_salsa
2743                CALL pmci_interp_1sto_all ( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),              &
2744                                            kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2745             ENDDO
2746          ENDIF
2747       ENDIF
2748
2749       IF ( topography /= 'flat' )  THEN
2750!
2751!--       Inside buildings set velocities back to zero.
2752          DO  ic = nxlg, nxrg
2753             DO  jc = nysg, nyng
2754                DO  kc = nzb, nzt
2755                   u(kc,jc,ic)   = MERGE( u(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 1 ) )
2756                   v(kc,jc,ic)   = MERGE( v(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 2 ) )
2757                   w(kc,jc,ic)   = MERGE( w(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 3 ) )
2758                   u_p(kc,jc,ic) = MERGE( u_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 1 ) )
2759                   v_p(kc,jc,ic) = MERGE( v_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 2 ) )
2760                   w_p(kc,jc,ic) = MERGE( w_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_static_0(kc,jc,ic), 3 ) )
2761                ENDDO
2762             ENDDO
2763          ENDDO
2764       ENDIF
2765    ENDIF
2766
2767
2768 CONTAINS
2769
2770
2771    SUBROUTINE pmci_interp_1sto_all( child_array, parent_array, kct, ifl, ifu, jfl, jfu, kfl, kfu,  &
2772         var )
2773!
2774!--    Interpolation of the internal values for the child-domain initialization
2775       IMPLICIT NONE
2776
2777       INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
2778
2779       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !<  Indicates start index of child cells belonging to certain
2780                                                               !<  parent cell - x direction
2781       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !<  Indicates end index of child cells belonging to certain
2782                                                               !<  parent cell - x direction
2783       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !<  Indicates start index of child cells belonging to certain
2784                                                               !<  parent cell - y direction
2785       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !<  Indicates end index of child cells belonging to certain
2786                                                               !<  parent cell - y direction
2787       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !<  Indicates start index of child cells belonging to certain
2788                                                               !<  parent cell - z direction
2789       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !<  Indicates end index of child cells belonging to certain
2790                                                               !<  parent cell - z direction
2791       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array  !<  Child-grid array
2792       REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN) ::  parent_array        !<  Parent-grid array
2793
2794       CHARACTER(LEN=1), INTENT(IN) ::  var  !<  Variable symbol: 'u', 'v', 'w' or 's'
2795!
2796!--    Local variables:
2797       INTEGER(iwp) ::  ic        !< Running child-grid index in the x-direction
2798       INTEGER(iwp) ::  icb       !< Index pointing to the first redundant ghost point layer behind the actual boundary
2799                                  !< ghost point layer in the x-direction
2800       INTEGER(iwp) ::  icbc      !< Index pointing to the boundary ghost point layer in the x-direction
2801       INTEGER(iwp) ::  icfirst   !< Leftmost child-grid index initialized by the main loops (usually icfirst == icl_init)
2802       INTEGER(iwp) ::  iclast    !< Rightmost child-grid index initialized by the main loops (usually iclast == icr_init)
2803       INTEGER(iwp) ::  icl_init  !< Left child-grid index bound for initialization in the x-direction
2804       INTEGER(iwp) ::  icr_init  !< Right child-grid index bound for initialization in the x-direction
2805       INTEGER(iwp) ::  jc        !< Running child-grid index in the y-direction
2806       INTEGER(iwp) ::  jcb       !< Index pointing to the first redundant ghost point layer behind the actual boundary
2807                                  !< ghost point layer in the y-direction
2808       INTEGER(iwp) ::  jcbc      !< Index pointing to the boundary ghost point layer in the y-direction
2809       INTEGER(iwp) ::  jcfirst   !< Southmost child-grid index initialized by the main loops (usually jcfirst == jcs_init)
2810       INTEGER(iwp) ::  jclast    !< Northmost child-grid index initialized by the main loops (usually jclast == jcn_init)
2811       INTEGER(iwp) ::  jcs_init  !< South child-grid index bound for initialization in the y-direction
2812       INTEGER(iwp) ::  jcn_init  !< North child-grid index bound for initialization in the y-direction
2813       INTEGER(iwp) ::  kc        !< Running child-grid index in the z-direction
2814       INTEGER(iwp) ::  ip        !< Running parent-grid index in the x-direction
2815       INTEGER(iwp) ::  ipl_init  !< Left parent-grid index bound for initialization in the x-direction
2816       INTEGER(iwp) ::  ipr_init  !< Right parent-grid index bound for initialization in the x-direction
2817       INTEGER(iwp) ::  jp        !< Running parent-grid index in the y-direction
2818       INTEGER(iwp) ::  jps_init  !< South parent-grid index bound for initialization in the y-direction
2819       INTEGER(iwp) ::  jpn_init  !< North parent-grid index bound for initialization in the y-direction
2820       INTEGER(iwp) ::  kp        !< Running parent-grid index in the z-direction
2821
2822
2823       ipl_init = ipl
2824       ipr_init = ipr
2825       jps_init = jps
2826       jpn_init = jpn
2827       icl_init = nxl
2828       icr_init = nxr
2829       jcs_init = nys
2830       jcn_init = nyn
2831
2832       icbc = -1
2833       icb  = -2
2834       jcbc = -1
2835       jcb  = -2
2836       IF ( var == 'u' )  THEN
2837          icbc =  0
2838          icb  = -1
2839       ELSE IF ( var == 'v' )  THEN
2840          jcbc =  0
2841          jcb  = -1
2842       ENDIF
2843       
2844       IF ( nesting_mode /= 'vertical' )  THEN
2845          IF ( bc_dirichlet_l )  THEN
2846             ipl_init = ipl + 1
2847             icl_init = nxl - 1
2848!
2849!--          For u, nxl is a ghost node, but not for the other variables
2850             IF ( var == 'u' )  THEN
2851                ipl_init = ipl + 2
2852                icl_init = nxl
2853             ENDIF
2854          ENDIF
2855          IF ( bc_dirichlet_s )  THEN
2856             jps_init = jps + 1
2857             jcs_init = nys - 1 
2858!
2859!--          For v, nys is a ghost node, but not for the other variables
2860             IF ( var == 'v' )  THEN
2861                jps_init = jps + 2
2862                jcs_init = nys
2863             ENDIF
2864          ENDIF
2865          IF ( bc_dirichlet_r )  THEN
2866             ipr_init = ipr - 1
2867             icr_init = nxr + 1 
2868          ENDIF
2869          IF ( bc_dirichlet_n )  THEN
2870             jpn_init = jpn - 1
2871             jcn_init = nyn + 1
2872          ENDIF
2873       ENDIF     
2874
2875       child_array(:,:,:) = 0.0_wp
2876
2877       IF ( var == 'u' )  THEN
2878
2879          icfirst = ifl(ipl_init) 
2880          iclast  = ifl(ipr_init+1) - 1
2881          jcfirst = jfl(jps_init)
2882          jclast  = jfu(jpn_init)
2883          DO  ip = ipl_init, ipr_init
2884             DO  jp = jps_init, jpn_init
2885                DO  kp = 0, kct + 1 
2886
2887                   DO  ic = ifl(ip), ifl(ip+1)-1
2888                      DO  jc = jfl(jp), jfu(jp)
2889                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
2890                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
2891                         ENDDO
2892                      ENDDO
2893                   ENDDO
2894                   
2895                ENDDO
2896             ENDDO
2897          ENDDO
2898
2899       ELSE IF ( var == 'v' )  THEN
2900
2901          icfirst = ifl(ipl_init) 
2902          iclast  = ifu(ipr_init)
2903          jcfirst = jfl(jps_init)
2904          jclast  = jfl(jpn_init+1) - 1
2905          DO  ip = ipl_init, ipr_init
2906             DO  jp = jps_init, jpn_init
2907                DO  kp = 0, kct + 1 
2908
2909                   DO  ic = ifl(ip), ifu(ip)
2910                      DO  jc = jfl(jp), jfl(jp+1)-1
2911                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
2912                            child_array(kc,jc,ic) = parent_array(kp,jp,ip) 
2913                         ENDDO
2914                      ENDDO
2915                   ENDDO
2916
2917                ENDDO
2918             ENDDO
2919          ENDDO
2920
2921       ELSE IF ( var == 'w' )  THEN
2922
2923          icfirst = ifl(ipl_init) 
2924          iclast  = ifu(ipr_init)
2925          jcfirst = jfl(jps_init)
2926          jclast  = jfu(jpn_init)
2927          DO  ip = ipl_init, ipr_init
2928             DO  jp = jps_init, jpn_init
2929                DO  kp = 1, kct + 1 
2930
2931                   DO  ic = ifl(ip), ifu(ip)
2932                      DO  jc = jfl(jp), jfu(jp)
2933!                         
2934!--                      Because the kp-loop for w starts from kp=1 instead of 0
2935                         child_array(nzb,jc,ic) = 0.0_wp
2936                         DO  kc = kfu(kp-1)+1, kfu(kp) 
2937                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
2938                         ENDDO
2939                      ENDDO
2940                   ENDDO
2941                   
2942                ENDDO
2943             ENDDO
2944          ENDDO
2945
2946       ELSE   ! scalars
2947
2948          icfirst = ifl(ipl_init) 
2949          iclast  = ifu(ipr_init)
2950          jcfirst = jfl(jps_init)
2951          jclast  = jfu(jpn_init)
2952          DO  ip = ipl_init, ipr_init
2953             DO  jp = jps_init, jpn_init
2954                DO  kp = 0, kct + 1
2955                                     
2956                   DO  ic = ifl(ip), ifu(ip)
2957                      DO  jc = jfl(jp), jfu(jp)                         
2958                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
2959                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
2960                         ENDDO
2961                      ENDDO
2962                   ENDDO
2963                   
2964                ENDDO
2965             ENDDO
2966          ENDDO
2967
2968       ENDIF  ! var 
2969!
2970!--    If the number of grid points in child subdomain in x- or y-direction
2971!--    (nxr - nxl + 1 and/or nyn - nys + 1) is not integer divisible by the grid spacing
2972!--    ratio in its direction (igsr and/or jgsr), the above loops will return with
2973!--    unfilled gaps in the initial fields. These gaps, if present, are filled here. 
2974       IF ( icfirst > icl_init )  THEN
2975          DO  ic = icl_init, icfirst - 1
2976             child_array(:,:,ic) = child_array(:,:,icfirst)
2977          ENDDO
2978       ENDIF
2979       IF ( iclast < icr_init )  THEN
2980          DO  ic = iclast + 1, icr_init
2981             child_array(:,:,ic) = child_array(:,:,iclast)
2982          ENDDO
2983       ENDIF
2984       IF ( jcfirst > jcs_init )  THEN
2985          DO  jc = jcs_init, jcfirst - 1
2986             child_array(:,jc,:) = child_array(:,jcfirst,:)
2987          ENDDO
2988       ENDIF
2989       IF ( jclast < jcn_init )  THEN
2990          DO  jc = jclast + 1, jcn_init
2991             child_array(:,jc,:) = child_array(:,jclast,:)
2992          ENDDO
2993       ENDIF
2994!
2995!--    Finally, make sure that also the redundant 2nd and 3rd ghost-node layers
2996!--    including the corners are properly filled up.
2997       IF ( nys == 0 )  THEN
2998          DO  jc = -nbgp, jcb  ! jcb = -2 if var == v, else jcb = -1
2999             child_array(0:nzt+1,jc,nxlg:nxrg) = child_array(0:nzt+1,jcbc,nxlg:nxrg)
3000          ENDDO         
3001       ENDIF
3002       IF ( nyn == ny )  THEN
3003          DO  jc = ny+2, ny+nbgp
3004             child_array(0:nzt+1,jc,nxlg:nxrg) = child_array(0:nzt+1,ny+1,nxlg:nxrg)
3005          ENDDO
3006       ENDIF
3007       IF ( nxl == 0 )  THEN
3008          DO  ic = -nbgp, icb  ! icb = -2 if var == u, else icb = -1
3009             child_array(0:nzt+1,nysg:nyng,ic) = child_array(0:nzt+1,nysg:nyng,icbc)
3010          ENDDO         
3011       ENDIF
3012       IF ( nxr == nx )  THEN
3013          DO  ic = nx+2, nx+nbgp
3014             child_array(0:nzt+1,nysg:nyng,ic) = child_array(0:nzt+1,nysg:nyng,nx+1)
3015          ENDDO   
3016       ENDIF
3017
3018    END SUBROUTINE pmci_interp_1sto_all
3019
3020#endif
3021 END SUBROUTINE pmci_child_initialize
3022
3023
3024
3025 SUBROUTINE pmci_check_setting_mismatches
3026!
3027!-- Check for mismatches between settings of root and child variables
3028!-- (e.g., all children have to follow the end_time settings of the root model).
3029!-- The root model overwrites variables in the other models, so these variables
3030!-- only need to be set once in file PARIN.
3031
3032#if defined( __parallel )
3033
3034    USE control_parameters,                                                                         &
3035        ONLY:  dt_restart, end_time, message_string, restart_time, time_restart
3036
3037    IMPLICIT NONE
3038
3039    INTEGER ::  ierr                 !<  MPI error code
3040
3041    REAL(wp) ::  dt_restart_root     !<
3042    REAL(wp) ::  end_time_root       !< 
3043    REAL(wp) ::  restart_time_root   !<
3044    REAL(wp) ::  time_restart_root   !< 
3045
3046!
3047!-- Check the time to be simulated.
3048!-- Here, and in the following, the root process communicates the respective
3049!-- variable to all others, and its value will then be compared with the local
3050!-- values.
3051    IF ( pmc_is_rootmodel() )  end_time_root = end_time
3052    CALL MPI_BCAST( end_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3053
3054    IF ( .NOT. pmc_is_rootmodel() )  THEN
3055       IF ( end_time /= end_time_root )  THEN
3056          WRITE( message_string, * )  'mismatch between root model and ',                           &
3057               'child settings:& end_time(root) = ', end_time_root,                                 &
3058               '& end_time(child) = ', end_time, '& child value is set',                            &
3059               ' to root value'
3060          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6,                      &
3061                        0 )
3062          end_time = end_time_root
3063       ENDIF
3064    ENDIF
3065!
3066!-- Same for restart time
3067    IF ( pmc_is_rootmodel() )  restart_time_root = restart_time
3068    CALL MPI_BCAST( restart_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3069
3070    IF ( .NOT. pmc_is_rootmodel() )  THEN
3071       IF ( restart_time /= restart_time_root )  THEN
3072          WRITE( message_string, * )  'mismatch between root model and ',      &
3073               'child settings: & restart_time(root) = ', restart_time_root,   &
3074               '& restart_time(child) = ', restart_time, '& child ',           &
3075               'value is set to root value'
3076          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3077                        0 )
3078          restart_time = restart_time_root
3079       ENDIF
3080    ENDIF
3081!
3082!-- Same for dt_restart
3083    IF ( pmc_is_rootmodel() )  dt_restart_root = dt_restart
3084    CALL MPI_BCAST( dt_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3085
3086    IF ( .NOT. pmc_is_rootmodel() )  THEN
3087       IF ( dt_restart /= dt_restart_root )  THEN
3088          WRITE( message_string, * )  'mismatch between root model and ',      &
3089               'child settings: & dt_restart(root) = ', dt_restart_root,       &
3090               '& dt_restart(child) = ', dt_restart, '& child ',               &
3091               'value is set to root value'
3092          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3093                        0 )
3094          dt_restart = dt_restart_root
3095       ENDIF
3096    ENDIF
3097!
3098!-- Same for time_restart
3099    IF ( pmc_is_rootmodel() )  time_restart_root = time_restart
3100    CALL MPI_BCAST( time_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3101
3102    IF ( .NOT. pmc_is_rootmodel() )  THEN
3103       IF ( time_restart /= time_restart_root )  THEN
3104          WRITE( message_string, * )  'mismatch between root model and ',      &
3105               'child settings: & time_restart(root) = ', time_restart_root,   &
3106               '& time_restart(child) = ', time_restart, '& child ',           &
3107               'value is set to root value'
3108          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3109                        0 )
3110          time_restart = time_restart_root
3111       ENDIF
3112    ENDIF
3113
3114#endif
3115
3116 END SUBROUTINE pmci_check_setting_mismatches
3117
3118
3119
3120 SUBROUTINE pmci_synchronize
3121
3122#if defined( __parallel )
3123!
3124!-- Unify the time steps for each model and synchronize using
3125!-- MPI_ALLREDUCE with the MPI_MIN operator over all processes using
3126!-- the global communicator MPI_COMM_WORLD.
3127   
3128   IMPLICIT NONE
3129
3130   INTEGER(iwp) ::  ierr  !<  MPI error code
3131   REAL(wp) ::  dtl       !<  Local time step of the current process
3132   REAL(wp) ::  dtg       !<  Global time step defined as the global minimum of dtl of all processes
3133
3134
3135   IF ( debug_output_timestep )  CALL debug_message( 'pmci_synchronize', 'start' )
3136   
3137   dtl = dt_3d
3138   CALL MPI_ALLREDUCE( dtl, dtg, 1, MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
3139   dt_3d  = dtg
3140
3141   IF ( debug_output_timestep )  CALL debug_message( 'pmci_synchronize', 'end' )
3142
3143#endif
3144 END SUBROUTINE pmci_synchronize
3145               
3146
3147
3148 SUBROUTINE pmci_set_swaplevel( swaplevel )
3149
3150!
3151!-- After each Runge-Kutta sub-timestep, alternately set buffer one or buffer
3152!-- two active
3153
3154    IMPLICIT NONE
3155
3156    INTEGER(iwp), INTENT(IN) ::  swaplevel  !< swaplevel (1 or 2) of PALM's timestep
3157
3158    INTEGER(iwp) ::  child_id    !<  Child id of the child number m
3159    INTEGER(iwp) ::  m           !<  Loop index over all children of the current parent
3160
3161#if defined( __parallel )
3162    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
3163       child_id = pmc_parent_for_child(m)
3164       CALL pmc_s_set_active_data_array( child_id, swaplevel )
3165    ENDDO
3166#endif
3167 END SUBROUTINE pmci_set_swaplevel
3168
3169
3170
3171 SUBROUTINE pmci_datatrans( local_nesting_mode )   
3172!
3173!-- This subroutine controls the nesting according to the nestpar
3174!-- parameter nesting_mode (two-way (default) or one-way) and the
3175!-- order of anterpolations according to the nestpar parameter
3176!-- nesting_datatransfer_mode (cascade, overlap or mixed (default)).
3177!-- Although nesting_mode is a variable of this model, pass it as
3178!-- an argument to allow for example to force one-way initialization
3179!-- phase.
3180!-- Note that interpolation ( parent_to_child ) must always be carried
3181!-- out before anterpolation ( child_to_parent ).
3182
3183    IMPLICIT NONE
3184
3185    CHARACTER(LEN=*), INTENT(IN) ::  local_nesting_mode  !<  Nesting mode: 'one-way', 'two-way' or 'vertical'
3186
3187#if defined( __parallel )   
3188
3189    IF ( debug_output_timestep )  CALL debug_message( 'pmci_datatrans', 'start' )
3190
3191    IF ( TRIM( local_nesting_mode ) == 'one-way' )  THEN
3192
3193       CALL pmci_child_datatrans( parent_to_child )
3194       CALL pmci_parent_datatrans( parent_to_child )
3195
3196    ELSE
3197
3198       IF ( nesting_datatransfer_mode == 'cascade' )  THEN
3199
3200          CALL pmci_child_datatrans( parent_to_child )
3201          CALL pmci_parent_datatrans( parent_to_child )
3202
3203          CALL pmci_parent_datatrans( child_to_parent )
3204          CALL pmci_child_datatrans( child_to_parent )
3205
3206       ELSEIF ( nesting_datatransfer_mode == 'overlap')  THEN
3207
3208          CALL pmci_parent_datatrans( parent_to_child )
3209          CALL pmci_child_datatrans( parent_to_child )
3210
3211          CALL pmci_child_datatrans( child_to_parent )
3212          CALL pmci_parent_datatrans( child_to_parent )
3213
3214       ELSEIF ( TRIM( nesting_datatransfer_mode ) == 'mixed' )  THEN
3215
3216          CALL pmci_parent_datatrans( parent_to_child )
3217          CALL pmci_child_datatrans( parent_to_child )
3218
3219          CALL pmci_parent_datatrans( child_to_parent )
3220          CALL pmci_child_datatrans( child_to_parent )
3221
3222       ENDIF
3223
3224    ENDIF
3225
3226    IF ( debug_output_timestep )  CALL debug_message( 'pmci_datatrans', 'end' )
3227
3228#endif
3229 END SUBROUTINE pmci_datatrans
3230
3231
3232
3233 SUBROUTINE pmci_parent_datatrans( direction )
3234   
3235    IMPLICIT NONE
3236
3237    INTEGER(iwp), INTENT(IN) ::  direction   !<  Direction of the data transfer: 'parent_to_child' or 'child_to_parent'
3238
3239#if defined( __parallel )
3240    INTEGER(iwp) ::  child_id    !<  Child id of the child number m
3241    INTEGER(iwp) ::  i           !<  Parent-grid index in x-direction
3242    INTEGER(iwp) ::  j           !<  Parent-grid index in y-direction
3243    INTEGER(iwp) ::  k           !<  Parent-grid index in z-direction
3244    INTEGER(iwp) ::  m           !<  Loop index over all children of the current parent
3245
3246
3247    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
3248       child_id = pmc_parent_for_child(m)
3249       IF ( direction == parent_to_child )  THEN
3250          CALL cpu_log( log_point_s(71), 'pmc parent send', 'start' )
3251          CALL pmc_s_fillbuffer( child_id )
3252          CALL cpu_log( log_point_s(71), 'pmc parent send', 'stop' )
3253       ELSE
3254!
3255!--       Communication from child to parent
3256          CALL cpu_log( log_point_s(72), 'pmc parent recv', 'start' )
3257          child_id = pmc_parent_for_child(m)
3258          CALL pmc_s_getdata_from_buffer( child_id )
3259          CALL cpu_log( log_point_s(72), 'pmc parent recv', 'stop' )
3260!
3261!--       The anterpolated data is now available in u etc
3262          IF ( topography /= 'flat' )  THEN
3263!
3264!--          Inside buildings/topography reset velocities back to zero.
3265!--          Scalars (pt, q, s, km, kh, p, sa, ...) are ignored at
3266!--          present, maybe revise later.
3267!--          Resetting of e is removed as unnecessary since e is not
3268!--          anterpolated, and as incorrect since it overran the default
3269!--          Neumann condition (bc_e_b).
3270             DO   i = nxlg, nxrg
3271                DO   j = nysg, nyng
3272                   DO  k = nzb, nzt+1
3273                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 1 ) )
3274                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 2 ) )
3275                      w(k,j,i) = MERGE( w(k,j,i), 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 3 ) )
3276!
3277!--                 TO_DO: zero setting of temperature within topography creates
3278!--                       wrong results
3279!                   pt(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
3280!                   IF ( humidity  .OR.  passive_scalar )  THEN
3281!                      q(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
3282!                   ENDIF
3283                   ENDDO
3284                ENDDO
3285             ENDDO
3286          ENDIF
3287       ENDIF
3288    ENDDO  ! m
3289
3290#endif
3291 END SUBROUTINE pmci_parent_datatrans
3292
3293
3294
3295 SUBROUTINE pmci_child_datatrans( direction )
3296
3297    IMPLICIT NONE
3298
3299    INTEGER(iwp), INTENT(IN) ::  direction  !< Transfer direction: parent_to_child or child_to_parent
3300
3301#if defined( __parallel )
3302
3303    REAL(wp), DIMENSION(1) ::  dtl          !< Time step size
3304
3305
3306    dtl = dt_3d
3307    IF ( .NOT.  pmc_is_rootmodel() )  THEN
3308
3309       IF ( direction == parent_to_child )  THEN
3310   
3311          CALL cpu_log( log_point_s(73), 'pmc child recv', 'start' )
3312          CALL pmc_c_getbuffer( )
3313          CALL cpu_log( log_point_s(73), 'pmc child recv', 'stop' )
3314
3315          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'start' )
3316          CALL pmci_interpolation
3317          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' )
3318     
3319       ELSE
3320!
3321!--       direction == child_to_parent
3322          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'start' )
3323          CALL pmci_anterpolation
3324          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'stop' )
3325
3326          CALL cpu_log( log_point_s(74), 'pmc child send', 'start' )
3327          CALL pmc_c_putbuffer( )
3328          CALL cpu_log( log_point_s(74), 'pmc child send', 'stop' )
3329
3330       ENDIF
3331    ENDIF
3332
3333 CONTAINS
3334
3335   
3336    SUBROUTINE pmci_interpolation
3337
3338!
3339!--    A wrapper routine for all interpolation actions
3340     
3341       IMPLICIT NONE
3342
3343       INTEGER(iwp) ::  ibgp       !< Index running over the nbgp boundary ghost points in i-direction
3344       INTEGER(iwp) ::  jbgp       !< Index running over the nbgp boundary ghost points in j-direction
3345       INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
3346       INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
3347       INTEGER(iwp) ::  lg         !< Running index for salsa gases
3348       INTEGER(iwp) ::  n          !< Running index for number of chemical species
3349     
3350!
3351!--    In case of vertical nesting no interpolation is needed for the
3352!--    horizontal boundaries
3353       IF ( nesting_mode /= 'vertical' )  THEN
3354!
3355!--       Left border pe:
3356          IF ( bc_dirichlet_l )  THEN
3357
3358             CALL pmci_interp_1sto_lr( u, uc, kcto, jflo, jfuo, kflo, kfuo, 'l', 'u' )
3359             CALL pmci_interp_1sto_lr( v, vc, kcto, jflv, jfuv, kflo, kfuo, 'l', 'v' )
3360             CALL pmci_interp_1sto_lr( w, wc, kctw, jflo, jfuo, kflw, kfuw, 'l', 'w' )
3361
3362             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
3363                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
3364                     .NOT.  constant_diffusion ) )  THEN
3365!                CALL pmci_interp_1sto_lr( e, ec, kcto, jflo, jfuo, kflo, kfuo, 'l', 'e' )
3366!
3367!--             Interpolation of e is replaced by the Neumann condition.
3368                DO  ibgp = -nbgp, -1
3369                   e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,0)
3370                ENDDO
3371
3372             ENDIF
3373
3374             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3375                CALL pmci_interp_1sto_lr( diss, dissc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )
3376             ENDIF
3377
3378             IF ( .NOT. neutral )  THEN
3379                CALL pmci_interp_1sto_lr( pt, ptc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )
3380             ENDIF
3381
3382             IF ( humidity )  THEN
3383
3384                CALL pmci_interp_1sto_lr( q, q_c, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )
3385
3386                IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3387                   CALL pmci_interp_1sto_lr( qc, qcc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' ) 
3388                   CALL pmci_interp_1sto_lr( nc, ncc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )         
3389                ENDIF
3390
3391                IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3392                   CALL pmci_interp_1sto_lr( qr, qrc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' ) 
3393                   CALL pmci_interp_1sto_lr( nr, nrc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )             
3394                ENDIF
3395
3396             ENDIF
3397
3398             IF ( passive_scalar )  THEN
3399                CALL pmci_interp_1sto_lr( s, sc, kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )
3400             ENDIF
3401
3402             IF ( air_chemistry  .AND.  nesting_chem )  THEN
3403                DO  n = 1, nspec
3404                   CALL pmci_interp_1sto_lr( chem_species(n)%conc, chem_spec_c(:,:,:,n),            &
3405                        kcto, jflo, jfuo, kflo, kfuo, 'l', 's' )
3406                ENDDO
3407             ENDIF
3408
3409             IF ( salsa  .AND.  nesting_salsa )  THEN
3410                DO  lb = 1, nbins_aerosol
3411                   CALL pmci_interp_1sto_lr( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
3412                                             kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
3413                ENDDO
3414                DO  lc = 1, nbins_aerosol * ncomponents_mass
3415                   CALL pmci_interp_1sto_lr( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
3416                                             kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
3417                ENDDO
3418                IF ( .NOT. salsa_gases_from_chem )  THEN
3419                   DO  lg = 1, ngases_salsa
3420                      CALL pmci_interp_1sto_lr( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
3421                                                kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
3422                   ENDDO
3423                ENDIF
3424             ENDIF
3425
3426          ENDIF
3427!
3428!--       Right border pe
3429          IF ( bc_dirichlet_r )  THEN
3430             
3431             CALL pmci_interp_1sto_lr( u, uc, kcto, jflo, jfuo, kflo, kfuo, 'r', 'u' )
3432             CALL pmci_interp_1sto_lr( v, vc, kcto, jflv, jfuv, kflo, kfuo, 'r', 'v' )
3433             CALL pmci_interp_1sto_lr( w, wc, kctw, jflo, jfuo, kflw, kfuw, 'r', 'w' )
3434
3435             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
3436                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
3437                     .NOT.  constant_diffusion ) )  THEN
3438!                CALL pmci_interp_1sto_lr( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 'r', 'e' )
3439!
3440!--             Interpolation of e is replaced by the Neumann condition.
3441                DO  ibgp = nx+1, nx+nbgp
3442                   e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,nx)
3443                ENDDO
3444             ENDIF
3445
3446             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3447                CALL pmci_interp_1sto_lr( diss, dissc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3448             ENDIF
3449
3450             IF (  .NOT.  neutral )  THEN
3451                CALL pmci_interp_1sto_lr( pt, ptc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3452             ENDIF
3453
3454             IF ( humidity )  THEN
3455                CALL pmci_interp_1sto_lr( q, q_c, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3456
3457                IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3458                   CALL pmci_interp_1sto_lr( qc, qcc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' ) 
3459                   CALL pmci_interp_1sto_lr( nc, ncc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3460                ENDIF
3461
3462                IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3463                   CALL pmci_interp_1sto_lr( qr, qrc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' ) 
3464                   CALL pmci_interp_1sto_lr( nr, nrc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' ) 
3465                ENDIF
3466
3467             ENDIF
3468
3469             IF ( passive_scalar )  THEN
3470                CALL pmci_interp_1sto_lr( s, sc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3471             ENDIF
3472
3473             IF ( air_chemistry  .AND.  nesting_chem )  THEN
3474                DO  n = 1, nspec
3475                   CALL pmci_interp_1sto_lr( chem_species(n)%conc, chem_spec_c(:,:,:,n),            &
3476                        kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3477                ENDDO
3478             ENDIF
3479
3480             IF ( salsa  .AND.  nesting_salsa )  THEN
3481                DO  lb = 1, nbins_aerosol
3482                   CALL pmci_interp_1sto_lr( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
3483                                             kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3484                ENDDO
3485                DO  lc = 1, nbins_aerosol * ncomponents_mass
3486                   CALL pmci_interp_1sto_lr( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
3487                                             kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3488                ENDDO
3489                IF ( .NOT. salsa_gases_from_chem )  THEN
3490                   DO  lg = 1, ngases_salsa
3491                      CALL pmci_interp_1sto_lr( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
3492                                                kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
3493                   ENDDO
3494                ENDIF
3495             ENDIF
3496
3497          ENDIF
3498!
3499!--       South border pe
3500          IF ( bc_dirichlet_s )  THEN
3501
3502             CALL pmci_interp_1sto_sn( v, vc, kcto, iflo, ifuo, kflo, kfuo, 's', 'v' )
3503             CALL pmci_interp_1sto_sn( w, wc, kctw, iflo, ifuo, kflw, kfuw, 's', 'w' )
3504             CALL pmci_interp_1sto_sn( u, uc, kcto, iflu, ifuu, kflo, kfuo, 's', 'u' )
3505
3506             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
3507                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
3508                     .NOT.  constant_diffusion ) )  THEN
3509!                CALL pmci_interp_1sto_sn( e, ec, kcto, iflo, ifuo, kflo, kfuo, 's', 'e' )
3510!
3511!--             Interpolation of e is replaced by the Neumann condition.
3512                DO  jbgp = -nbgp, -1
3513                   e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,0,nxl:nxr)
3514                ENDDO
3515             ENDIF
3516
3517             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3518                CALL pmci_interp_1sto_sn( diss, dissc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3519             ENDIF
3520
3521             IF (  .NOT.  neutral )  THEN
3522                CALL pmci_interp_1sto_sn( pt, ptc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3523             ENDIF
3524
3525             IF ( humidity )  THEN
3526                CALL pmci_interp_1sto_sn( q, q_c, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3527
3528                IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3529                   CALL pmci_interp_1sto_sn( qc, qcc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3530                   CALL pmci_interp_1sto_sn( nc, ncc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3531                ENDIF
3532
3533                IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3534                   CALL pmci_interp_1sto_sn( qr, qrc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3535                   CALL pmci_interp_1sto_sn( nr, nrc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3536                ENDIF
3537
3538             ENDIF
3539
3540             IF ( passive_scalar )  THEN
3541                CALL pmci_interp_1sto_sn( s, sc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3542             ENDIF
3543
3544             IF ( air_chemistry  .AND.  nesting_chem )  THEN
3545                DO  n = 1, nspec
3546                   CALL pmci_interp_1sto_sn( chem_species(n)%conc, chem_spec_c(:,:,:,n),            &
3547                        kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3548                ENDDO
3549             ENDIF
3550             
3551             IF ( salsa  .AND.  nesting_salsa )  THEN
3552                DO  lb = 1, nbins_aerosol
3553                   CALL pmci_interp_1sto_sn( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
3554                                             kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3555                ENDDO
3556                DO  lc = 1, nbins_aerosol * ncomponents_mass
3557                   CALL pmci_interp_1sto_sn( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
3558                                             kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3559                ENDDO
3560                IF ( .NOT. salsa_gases_from_chem )  THEN
3561                   DO  lg = 1, ngases_salsa
3562                      CALL pmci_interp_1sto_sn( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
3563                                                kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
3564                   ENDDO
3565                ENDIF
3566             ENDIF
3567                       
3568          ENDIF
3569!
3570!--       North border pe
3571          IF ( bc_dirichlet_n )  THEN
3572             
3573             CALL pmci_interp_1sto_sn( v, vc, kcto, iflo, ifuo, kflo, kfuo, 'n', 'v' )
3574             CALL pmci_interp_1sto_sn( w, wc, kctw, iflo, ifuo, kflw, kfuw, 'n', 'w' )
3575             CALL pmci_interp_1sto_sn( u, uc, kcto, iflu, ifuu, kflo, kfuo, 'n', 'u' )
3576
3577             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       & 
3578                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
3579                     .NOT.  constant_diffusion ) )  THEN
3580!                CALL pmci_interp_1sto_sn( e, ec, kcto, iflo, ifuo, kflo, kfuo, 'n', 'e' )
3581!
3582!--             Interpolation of e is replaced by the Neumann condition.
3583                DO  jbgp = ny+1, ny+nbgp
3584                   e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,ny,nxl:nxr)
3585                ENDDO
3586             ENDIF
3587
3588             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3589                CALL pmci_interp_1sto_sn( diss, dissc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3590             ENDIF
3591
3592             IF (  .NOT.  neutral )  THEN
3593                CALL pmci_interp_1sto_sn( pt, ptc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3594             ENDIF
3595
3596             IF ( humidity )  THEN
3597                CALL pmci_interp_1sto_sn( q, q_c, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3598
3599                IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3600                   CALL pmci_interp_1sto_sn( qc, qcc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3601                   CALL pmci_interp_1sto_sn( nc, ncc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3602                ENDIF
3603
3604                IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3605                   CALL pmci_interp_1sto_sn( qr, qrc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3606                   CALL pmci_interp_1sto_sn( nr, nrc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3607                ENDIF
3608
3609             ENDIF
3610
3611             IF ( passive_scalar )  THEN
3612                CALL pmci_interp_1sto_sn( s, sc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3613             ENDIF
3614
3615             IF ( air_chemistry  .AND.  nesting_chem )  THEN
3616                DO  n = 1, nspec
3617                   CALL pmci_interp_1sto_sn( chem_species(n)%conc, chem_spec_c(:,:,:,n),            &
3618                        kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3619                ENDDO
3620             ENDIF
3621             
3622             IF ( salsa  .AND.  nesting_salsa )  THEN
3623                DO  lb = 1, nbins_aerosol
3624                   CALL pmci_interp_1sto_sn( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
3625                                             kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3626                ENDDO
3627                DO  lc = 1, nbins_aerosol * ncomponents_mass
3628                   CALL pmci_interp_1sto_sn( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
3629                                             kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3630                ENDDO
3631                IF ( .NOT. salsa_gases_from_chem )  THEN
3632                   DO  lg = 1, ngases_salsa
3633                      CALL pmci_interp_1sto_sn( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
3634                                                kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
3635                   ENDDO
3636                ENDIF
3637             ENDIF
3638                         
3639          ENDIF
3640
3641       ENDIF       ! IF ( nesting_mode /= 'vertical' )
3642!
3643!--    All PEs are top-border PEs
3644       CALL pmci_interp_1sto_t( w, wc, kctw, iflo, ifuo, jflo, jfuo, 'w' )
3645       CALL pmci_interp_1sto_t( u, uc, kcto, iflu, ifuu, jflo, jfuo, 'u' )
3646       CALL pmci_interp_1sto_t( v, vc, kcto, iflo, ifuo, jflv, jfuv, 'v' )
3647
3648       IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                             &
3649            (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                              &
3650               .NOT.  constant_diffusion ) )  THEN
3651!          CALL pmci_interp_1sto_t( e, ec, kcto, iflo, ifuo, jflo, jfuo, 'e' )
3652!
3653!--       Interpolation of e is replaced by the Neumann condition.
3654          e(nzt+1,nys:nyn,nxl:nxr) = e(nzt,nys:nyn,nxl:nxr)
3655       ENDIF
3656
3657       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3658          CALL pmci_interp_1sto_t( diss, dissc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3659       ENDIF
3660
3661       IF ( .NOT. neutral )  THEN
3662          CALL pmci_interp_1sto_t( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3663       ENDIF
3664
3665       IF ( humidity )  THEN
3666          CALL pmci_interp_1sto_t( q, q_c, kcto, iflo, ifuo, jflo, jfuo, 's' )
3667          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3668             CALL pmci_interp_1sto_t( qc, qcc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3669             CALL pmci_interp_1sto_t( nc, ncc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3670          ENDIF
3671          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3672             CALL pmci_interp_1sto_t( qr, qrc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3673             CALL pmci_interp_1sto_t( nr, nrc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3674          ENDIF
3675       ENDIF
3676
3677       IF ( passive_scalar )  THEN
3678          CALL pmci_interp_1sto_t( s, sc, kcto, iflo, ifuo, jflo, jfuo, 's' )
3679       ENDIF
3680
3681       IF ( air_chemistry  .AND.  nesting_chem )  THEN
3682          DO  n = 1, nspec
3683             CALL pmci_interp_1sto_t( chem_species(n)%conc, chem_spec_c(:,:,:,n),                   &
3684                                      kcto, iflo, ifuo, jflo, jfuo, 's' )
3685          ENDDO
3686       ENDIF
3687       
3688       IF ( salsa  .AND.  nesting_salsa )  THEN
3689          DO  lb = 1, nbins_aerosol
3690             CALL pmci_interp_1sto_t( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),          &
3691                                      kcto, iflo, ifuo, jflo, jfuo, 's' )
3692          ENDDO
3693          DO  lc = 1, nbins_aerosol * ncomponents_mass
3694             CALL pmci_interp_1sto_t( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),              &
3695                                      kcto, iflo, ifuo, jflo, jfuo, 's' )
3696          ENDDO
3697          IF ( .NOT. salsa_gases_from_chem )  THEN
3698             DO  lg = 1, ngases_salsa
3699                CALL pmci_interp_1sto_t( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),                 &
3700                                         kcto, iflo, ifuo, jflo, jfuo, 's' )
3701             ENDDO
3702          ENDIF
3703       ENDIF
3704
3705   END SUBROUTINE pmci_interpolation
3706
3707
3708
3709   SUBROUTINE pmci_anterpolation
3710
3711!
3712!--   A wrapper routine for all anterpolation actions.
3713!--   Note that TKE is not anterpolated.
3714      IMPLICIT NONE
3715      INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
3716      INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
3717      INTEGER(iwp) ::  lg         !< Running index for salsa gases
3718      INTEGER(iwp) ::  n          !< Running index for number of chemical species
3719
3720     
3721      CALL pmci_anterp_tophat( u,  uc,  kcto, iflu, ifuu, jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' )
3722      CALL pmci_anterp_tophat( v,  vc,  kcto, iflo, ifuo, jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' )
3723      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' )
3724!
3725!--   Anterpolation of TKE and dissipation rate if parent and child are in
3726!--   RANS mode.
3727      IF ( rans_mode_parent  .AND.  rans_mode )  THEN
3728         CALL pmci_anterp_tophat( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' )
3729!
3730!--      Anterpolation of dissipation rate only if TKE-e closure is applied.
3731         IF ( rans_tke_e )  THEN
3732            CALL pmci_anterp_tophat( diss, dissc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo,         &
3733                 ijkfc_s, 'diss' )
3734         ENDIF
3735
3736      ENDIF
3737
3738      IF ( .NOT. neutral )  THEN
3739         CALL pmci_anterp_tophat( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 'pt' )
3740      ENDIF
3741
3742      IF ( humidity )  THEN
3743
3744         CALL pmci_anterp_tophat( q, q_c, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 'q' )
3745
3746         IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
3747
3748            CALL pmci_anterp_tophat( qc, qcc, kcto, iflo, ifuo, jflo, jfuo,                         &
3749                                     kflo, kfuo, ijkfc_s, 'qc' )
3750           
3751            CALL pmci_anterp_tophat( nc, ncc, kcto, iflo, ifuo, jflo, jfuo,                         &
3752                                     kflo, kfuo, ijkfc_s, 'nc' )
3753
3754         ENDIF
3755
3756         IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
3757
3758            CALL pmci_anterp_tophat( qr, qrc, kcto, iflo, ifuo, jflo, jfuo,                         &
3759                                     kflo, kfuo, ijkfc_s, 'qr' )
3760
3761            CALL pmci_anterp_tophat( nr, nrc, kcto, iflo, ifuo, jflo, jfuo,                         &
3762                                     kflo, kfuo, ijkfc_s, 'nr' )
3763
3764         ENDIF
3765
3766      ENDIF
3767
3768      IF ( passive_scalar )  THEN
3769         CALL pmci_anterp_tophat( s, sc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
3770      ENDIF
3771
3772      IF ( air_chemistry  .AND.  nesting_chem )  THEN
3773         DO  n = 1, nspec
3774            CALL pmci_anterp_tophat( chem_species(n)%conc, chem_spec_c(:,:,:,n),                    &
3775                                     kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
3776         ENDDO
3777      ENDIF
3778
3779      IF ( salsa  .AND.  nesting_salsa )  THEN
3780         DO  lb = 1, nbins_aerosol
3781            CALL pmci_anterp_tophat( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),           &
3782                                     kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
3783         ENDDO
3784         DO  lc = 1, nbins_aerosol * ncomponents_mass
3785            CALL pmci_anterp_tophat( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),               &
3786                                     kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
3787         ENDDO
3788         IF ( .NOT. salsa_gases_from_chem )  THEN
3789            DO  lg = 1, ngases_salsa
3790               CALL pmci_anterp_tophat( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),                  &
3791                                        kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
3792            ENDDO
3793         ENDIF
3794      ENDIF
3795
3796   END SUBROUTINE pmci_anterpolation
3797
3798
3799
3800   SUBROUTINE pmci_interp_1sto_lr( child_array, parent_array, kct, jfl, jfu, kfl, kfu, edge, var )
3801!
3802!--   Interpolation of ghost-node values used as the child-domain boundary
3803!--   conditions. This subroutine handles the left and right boundaries.
3804      IMPLICIT NONE
3805
3806      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
3807     
3808      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
3809                                                              !< parent cell - y direction
3810      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
3811                                                              !< parent cell - y direction
3812      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
3813                                                              !< parent cell - z direction
3814      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
3815                                                              !< parent cell - z direction
3816
3817      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
3818      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
3819
3820      CHARACTER(LEN=1), INTENT(IN) ::  edge                   !< Edge symbol: 'l' or 'r'
3821      CHARACTER(LEN=1), INTENT(IN) ::  var                    !< Variable symbol: 'u', 'v', 'w' or 's'     
3822!
3823!--   Local variables:
3824      INTEGER(iwp) ::  icb      !< Fixed child-grid index in the x-direction pointing to the node just behind the
3825                                !< boundary-value node
3826      INTEGER(iwp) ::  icbc     !< Fixed child-grid index in the x-direction pointing to the boundary-value nodes
3827      INTEGER(iwp) ::  icbgp    !< Index running over the redundant boundary ghost points in the x-direction
3828      INTEGER(iwp) ::  ierr     !< MPI error code
3829      INTEGER(iwp) ::  ipbeg    !< Parent-grid index in the x-direction pointing to the starting point of workarr_lr
3830                                !< in the parent-grid array
3831      INTEGER(iwp) ::  ipw      !< Reduced parent-grid index in the x-direction for workarr_lr pointing to
3832                                !< the boundary ghost node
3833      INTEGER(iwp) ::  ipwp     !< Reduced parent-grid index in the x-direction for workarr_lr pointing to
3834                                !< the first prognostic node
3835      INTEGER(iwp) ::  jc       !< Running child-grid index in the y-direction
3836      INTEGER(iwp) ::  jp       !< Running parent-grid index in the y-direction
3837      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction
3838      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction     
3839     
3840      REAL(wp) ::  cb           !< Interpolation coefficient for the boundary ghost node 
3841      REAL(wp) ::  cp           !< Interpolation coefficient for the first prognostic node
3842      REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
3843      REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
3844!
3845!--   Check which edge is to be handled
3846      IF ( edge == 'l' )  THEN
3847!
3848!--      For u, nxl is a ghost node, but not for the other variables
3849         IF ( var == 'u' )  THEN
3850            icbc  = nxl   
3851            icb   = icbc - 1
3852            ipw   = 2
3853            ipwp  = ipw        ! This is redundant when var == 'u'
3854            ipbeg = ipl
3855         ELSE
3856            icbc  = nxl - 1
3857            icb   = icbc - 1
3858            ipw   = 1
3859            ipwp  = 2
3860            ipbeg = ipl
3861         ENDIF       
3862      ELSEIF ( edge == 'r' )  THEN
3863         IF ( var == 'u' )  THEN
3864            icbc  = nxr + 1 
3865            icb   = icbc + 1
3866            ipw   = 1
3867            ipwp  = ipw        ! This is redundant when var == 'u'           
3868            ipbeg = ipr - 2
3869         ELSE
3870            icbc  = nxr + 1
3871            icb   = icbc + 1
3872            ipw   = 1
3873            ipwp  = 0
3874            ipbeg = ipr - 2
3875         ENDIF         
3876      ENDIF
3877!
3878!--   Interpolation coefficients
3879      IF ( interpolation_scheme_lrsn == 1 )  THEN
3880         cb = 1.0_wp  ! 1st-order upwind
3881      ELSE IF ( interpolation_scheme_lrsn == 2 )  THEN
3882         cb = 0.5_wp  ! 2nd-order central
3883      ELSE
3884         cb = 0.5_wp  ! 2nd-order central (default)
3885      ENDIF         
3886      cp = 1.0_wp - cb
3887!
3888!--   Substitute the necessary parent-grid data to the work array workarr_lr.
3889      workarr_lr = 0.0_wp     
3890      IF ( pdims(2) > 1 )  THEN
3891         
3892         IF ( bc_dirichlet_s )  THEN
3893            workarr_lr(0:pg%nz+1,jpsw:jpnw-1,0:2) = parent_array(0:pg%nz+1,jpsw:jpnw-1,ipbeg:ipbeg+2)
3894         ELSE IF ( bc_dirichlet_n )  THEN
3895            workarr_lr(0:pg%nz+1,jpsw+1:jpnw,0:2) = parent_array(0:pg%nz+1,jpsw+1:jpnw,ipbeg:ipbeg+2)
3896         ELSE
3897            workarr_lr(0:pg%nz+1,jpsw+1:jpnw-1,0:2)                                                 &
3898                 = parent_array(0:pg%nz+1,jpsw+1:jpnw-1,ipbeg:ipbeg+2)
3899         ENDIF
3900!
3901!--      South-north exchange if more than one PE subdomain in the y-direction.
3902!--      Note that in case of 3-D nesting the south (psouth == MPI_PROC_NULL)
3903!--      and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged
3904!--      because the nest domain is not cyclic.
3905!--      From south to north
3906         CALL MPI_SENDRECV( workarr_lr(0,jpsw+1,0), 1, workarr_lr_exchange_type, psouth,  0,        &
3907                            workarr_lr(0,jpnw,0), 1, workarr_lr_exchange_type, pnorth,  0, comm2d,  &
3908                            status, ierr )                             
3909!                                                                             
3910!--      From north to south                                                 
3911         CALL MPI_SENDRECV( workarr_lr(0,jpnw-1,0), 1, workarr_lr_exchange_type, pnorth,  1,        &
3912                            workarr_lr(0,jpsw,0), 1, workarr_lr_exchange_type, psouth,  1, comm2d,  &
3913                            status, ierr )                               
3914
3915      ELSE                                                                   
3916         workarr_lr(0:pg%nz+1,jpsw:jpnw,0:2) = parent_array(0:pg%nz+1,jpsw:jpnw,ipbeg:ipbeg+2)           
3917      ENDIF
3918
3919      IF ( var == 'u' )  THEN
3920
3921         DO  jp = jpsw, jpnw
3922            DO  kp = 0, kct
3923               
3924               DO  jc = jfl(jp), jfu(jp)
3925                  DO  kc = kfl(kp), kfu(kp)
3926                     child_array(kc,jc,icbc) = workarr_lr(kp,jp,ipw)
3927                  ENDDO
3928               ENDDO
3929
3930            ENDDO
3931         ENDDO
3932
3933      ELSE IF ( var == 'v' )  THEN
3934         
3935         DO  jp = jpsw, jpnw-1
3936            DO  kp = 0, kct
3937!
3938!--            First interpolate to the flux point
3939               c_interp_1 = cb * workarr_lr(kp,jp,ipw)   + cp * workarr_lr(kp,jp,ipwp)
3940               c_interp_2 = cb * workarr_lr(kp,jp+1,ipw) + cp * workarr_lr(kp,jp+1,ipwp)
3941!
3942!--            Use averages of the neighbouring matching grid-line values
3943               DO  jc = jfl(jp), jfl(jp+1)
3944                  child_array(kfl(kp):kfu(kp),jc,icbc) = 0.5_wp * ( c_interp_1 + c_interp_2 )
3945               ENDDO
3946!
3947!--            Then set the values along the matching grid-lines 
3948               IF  ( MOD( jfl(jp), jgsr ) == 0 )  THEN
3949                  child_array(kfl(kp):kfu(kp),jfl(jp),icbc) = c_interp_1
3950               ENDIF
3951            ENDDO
3952         ENDDO
3953!
3954!--      Finally, set the values along the last matching grid-line 
3955         IF ( MOD( jfl(jpnw), jgsr ) == 0 )  THEN
3956            DO  kp = 0, kct
3957               c_interp_1 = cb * workarr_lr(kp,jpnw,ipw) + cp * workarr_lr(kp,jpnw,ipwp)
3958               child_array(kfl(kp):kfu(kp),jfl(jpnw),icbc) = c_interp_1
3959            ENDDO
3960         ENDIF
3961!
3962!--      A gap may still remain in some cases if the subdomain size is not
3963!--      divisible by the grid-spacing ratio. In such a case, fill the
3964!--      gap. Note however, this operation may produce some additional
3965!--      momentum conservation error.
3966         IF ( jfl(jpnw) < nyn )  THEN
3967            DO  kp = 0, kct
3968               DO  jc = jfl(jpnw) + 1, nyn
3969                  child_array(kfl(kp):kfu(kp),jc,icbc) = child_array(kfl(kp):kfu(kp),jfl(jpnw),icbc)
3970               ENDDO
3971            ENDDO
3972         ENDIF
3973
3974      ELSE IF ( var == 'w' )  THEN
3975
3976         DO  jp = jpsw, jpnw
3977            DO  kp = 0, kct + 1   ! It is important to go up to kct+1 
3978!
3979!--            Interpolate to the flux point
3980               c_interp_1 = cb * workarr_lr(kp,jp,ipw) + cp * workarr_lr(kp,jp,ipwp)
3981!
3982!--            First substitute only the matching-node values
3983               child_array(kfu(kp),jfl(jp):jfu(jp),icbc) = c_interp_1
3984                 
3985            ENDDO
3986         ENDDO
3987
3988         DO  jp = jpsw, jpnw
3989            DO  kp = 1, kct + 1   ! It is important to go up to kct+1 
3990!
3991!--            Then fill up the nodes in between with the averages                 
3992               DO  kc = kfu(kp-1) + 1, kfu(kp) - 1 
3993                  child_array(kc,jfl(jp):jfu(jp),icbc) =                                            &
3994                       0.5_wp * ( child_array(kfu(kp-1),jfl(jp):jfu(jp),icbc)                       &
3995                       + child_array(kfu(kp),jfl(jp):jfu(jp),icbc) )
3996               ENDDO
3997                 
3998            ENDDO
3999         ENDDO
4000
4001      ELSE   ! any scalar
4002         
4003         DO  jp = jpsw, jpnw
4004            DO  kp = 0, kct
4005!
4006!--            Interpolate to the flux point
4007               c_interp_1 = cb * workarr_lr(kp,jp,ipw) + cp * workarr_lr(kp,jp,ipwp)
4008               DO  jc = jfl(jp), jfu(jp)
4009                  DO  kc = kfl(kp), kfu(kp)
4010                     child_array(kc,jc,icbc) = c_interp_1
4011                  ENDDO
4012               ENDDO
4013
4014            ENDDO
4015         ENDDO
4016
4017      ENDIF  ! var
4018!
4019!--   Fill up also the redundant 2nd and 3rd ghost-node layers
4020      IF ( edge == 'l' )  THEN
4021         DO  icbgp = -nbgp, icb
4022            child_array(0:nzt+1,nysg:nyng,icbgp) = child_array(0:nzt+1,nysg:nyng,icbc)
4023         ENDDO
4024      ELSEIF ( edge == 'r' )  THEN
4025         DO  icbgp = icb, nx+nbgp
4026            child_array(0:nzt+1,nysg:nyng,icbgp) = child_array(0:nzt+1,nysg:nyng,icbc)
4027         ENDDO
4028      ENDIF
4029
4030   END SUBROUTINE pmci_interp_1sto_lr
4031
4032
4033
4034   SUBROUTINE pmci_interp_1sto_sn( child_array, parent_array, kct, ifl, ifu, kfl, kfu, edge, var )
4035!
4036!--   Interpolation of ghost-node values used as the child-domain boundary
4037!--   conditions. This subroutine handles the south and north boundaries.
4038      IMPLICIT NONE
4039
4040      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
4041     
4042      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
4043                                                              !< parent cell - x direction
4044      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
4045                                                              !< parent cell - x direction
4046      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
4047                                                              !< parent cell - z direction
4048      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
4049                                                              !< parent cell - z direction
4050     
4051      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
4052      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
4053
4054      CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 's' or 'n'
4055      CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
4056!
4057!--   Local variables:     
4058      INTEGER(iwp) ::  ic       !< Running child-grid index in the x-direction
4059      INTEGER(iwp) ::  ierr     !< MPI error code
4060      INTEGER(iwp) ::  ip       !< Running parent-grid index in the x-direction
4061      INTEGER(iwp) ::  jcb      !< Fixed child-grid index in the y-direction pointing to the node just behind the
4062                                !< boundary-value node
4063      INTEGER(iwp) ::  jcbc     !< Fixed child-grid index in the y-direction pointing to the boundary-value nodes
4064      INTEGER(iwp) ::  jcbgp    !< Index running over the redundant boundary ghost points in y-direction
4065      INTEGER(iwp) ::  jpbeg    !< Parent-grid index in the y-direction pointing to the starting point of workarr_sn
4066                                !< in the parent-grid array
4067      INTEGER(iwp) ::  jpw      !< Reduced parent-grid index in the y-direction for workarr_sn pointing to
4068                                !< the boundary ghost node
4069      INTEGER(iwp) ::  jpwp     !< Reduced parent-grid index in the y-direction for workarr_sn pointing to
4070                                !< the first prognostic node
4071      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction     
4072      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction
4073      REAL(wp) ::  cb           !< Interpolation coefficient for the boundary ghost node 
4074      REAL(wp) ::  cp           !< Interpolation coefficient for the first prognostic node
4075      REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
4076      REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
4077
4078     
4079!
4080!--   Check which edge is to be handled: south or north
4081      IF ( edge == 's' )  THEN
4082!
4083!--      For v, nys is a ghost node, but not for the other variables
4084         IF ( var == 'v' )  THEN
4085            jcbc  = nys
4086            jcb   = jcbc - 1
4087            jpw   = 2
4088            jpwp  = 2         ! This is redundant when var == 'v'
4089            jpbeg = jps
4090         ELSE
4091            jcbc  = nys - 1
4092            jcb   = jcbc - 1
4093            jpw   = 1
4094            jpwp  = 2
4095            jpbeg = jps
4096         ENDIF
4097      ELSEIF ( edge == 'n' )  THEN
4098         IF ( var == 'v' )  THEN
4099            jcbc  = nyn + 1
4100            jcb   = jcbc + 1
4101            jpw   = 1
4102            jpwp  = 0         ! This is redundant when var == 'v'     
4103            jpbeg = jpn - 2
4104         ELSE
4105            jcbc  = nyn + 1
4106            jcb   = jcbc + 1
4107            jpw   = 1
4108            jpwp  = 0
4109            jpbeg = jpn - 2
4110         ENDIF
4111      ENDIF
4112!
4113!--   Interpolation coefficients
4114      IF ( interpolation_scheme_lrsn == 1 )  THEN
4115         cb = 1.0_wp  ! 1st-order upwind
4116      ELSE IF ( interpolation_scheme_lrsn == 2 )  THEN
4117         cb = 0.5_wp  ! 2nd-order central
4118      ELSE
4119         cb = 0.5_wp  ! 2nd-order central (default)
4120      ENDIF         
4121      cp = 1.0_wp - cb
4122!
4123!--   Substitute the necessary parent-grid data to the work array workarr_sn.
4124      workarr_sn = 0.0_wp     
4125      IF ( pdims(1) > 1 )  THEN
4126
4127         IF ( bc_dirichlet_l )  THEN
4128            workarr_sn(0:pg%nz+1,0:2,iplw:iprw-1) = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw:iprw-1)
4129         ELSE IF ( bc_dirichlet_r )  THEN
4130            workarr_sn(0:pg%nz+1,0:2,iplw+1:iprw) = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw+1:iprw)
4131         ELSE
4132            workarr_sn(0:pg%nz+1,0:2,iplw+1:iprw-1)                                                 &
4133                 = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw+1:iprw-1)
4134         ENDIF
4135!
4136!--      Left-right exchange if more than one PE subdomain in the x-direction.
4137!--      Note that in case of 3-D nesting the left (pleft == MPI_PROC_NULL) and
4138!--      right (pright == MPI_PROC_NULL) boundaries are not exchanged because
4139!--      the nest domain is not cyclic.
4140!--      From left to right
4141         CALL MPI_SENDRECV( workarr_sn(0,0,iplw+1), 1, workarr_sn_exchange_type, pleft,   0,        &
4142                            workarr_sn(0,0,iprw), 1, workarr_sn_exchange_type, pright, 0, comm2d,   &
4143                            status, ierr )
4144!                                                                           
4145!--      From right to left                                                 
4146         CALL MPI_SENDRECV( workarr_sn(0,0,iprw-1), 1, workarr_sn_exchange_type, pright,  1,        &
4147                            workarr_sn(0,0,iplw), 1, workarr_sn_exchange_type, pleft, 1, comm2d,    &
4148                            status, ierr )
4149
4150      ELSE
4151         workarr_sn(0:pg%nz+1,0:2,iplw:iprw) = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw:iprw)
4152      ENDIF
4153
4154      IF ( var == 'v' )  THEN
4155
4156         DO  ip = iplw, iprw
4157            DO  kp = 0, kct
4158               
4159               DO  ic = ifl(ip), ifu(ip)
4160                  DO  kc = kfl(kp), kfu(kp)
4161                     child_array(kc,jcbc,ic) = workarr_sn(kp,jpw,ip)
4162                  ENDDO
4163               ENDDO
4164
4165            ENDDO
4166         ENDDO
4167
4168      ELSE IF ( var == 'u' )  THEN
4169         
4170         DO  ip = iplw, iprw - 1
4171            DO  kp = 0, kct
4172!
4173!--            First interpolate to the flux point
4174               c_interp_1 = cb * workarr_sn(kp,jpw,ip)   + cp * workarr_sn(kp,jpwp,ip)
4175               c_interp_2 = cb * workarr_sn(kp,jpw,ip+1) + cp * workarr_sn(kp,jpwp,ip+1)
4176!
4177!--            Use averages of the neighbouring matching grid-line values
4178               DO  ic = ifl(ip), ifl(ip+1)
4179                  child_array(kfl(kp):kfu(kp),jcbc,ic) = 0.5_wp * ( c_interp_1 + c_interp_2 )
4180               ENDDO
4181!
4182!--            Then set the values along the matching grid-lines 
4183               IF ( MOD( ifl(ip), igsr ) == 0 )  THEN
4184                  child_array(kfl(kp):kfu(kp),jcbc,ifl(ip)) = c_interp_1
4185               ENDIF
4186
4187            ENDDO
4188         ENDDO
4189!
4190!--      Finally, set the values along the last matching grid-line 
4191         IF ( MOD( ifl(iprw), igsr ) == 0 )  THEN
4192            DO  kp = 0, kct
4193               c_interp_1 = cb * workarr_sn(kp,jpw,iprw) + cp * workarr_sn(kp,jpwp,iprw)
4194               child_array(kfl(kp):kfu(kp),jcbc,ifl(iprw)) = c_interp_1
4195            ENDDO
4196         ENDIF
4197!
4198!--      A gap may still remain in some cases if the subdomain size is not
4199!--      divisible by the grid-spacing ratio. In such a case, fill the
4200!--      gap. Note however, this operation may produce some additional
4201!--      momentum conservation error.
4202         IF ( ifl(iprw) < nxr )  THEN
4203            DO  kp = 0, kct
4204               DO  ic = ifl(iprw) + 1, nxr
4205                  child_array(kfl(kp):kfu(kp),jcbc,ic) = child_array(kfl(kp):kfu(kp),jcbc,ifl(iprw))
4206               ENDDO
4207            ENDDO
4208         ENDIF
4209
4210      ELSE IF ( var == 'w' )  THEN
4211
4212         DO  ip = iplw, iprw
4213            DO  kp = 0, kct + 1   ! It is important to go up to kct+1 
4214!
4215!--            Interpolate to the flux point
4216               c_interp_1 = cb * workarr_sn(kp,jpw,ip) + cp * workarr_sn(kp,jpwp,ip)
4217!
4218!--            First substitute only the matching-node values
4219               child_array(kfu(kp),jcbc,ifl(ip):ifu(ip)) = c_interp_1
4220
4221            ENDDO
4222         ENDDO
4223
4224         DO  ip = iplw, iprw
4225            DO  kp = 1, kct + 1   ! It is important to go up to kct + 1 
4226!
4227!--            Then fill up the nodes in between with the averages
4228               DO  kc = kfu(kp-1) + 1, kfu(kp) - 1 
4229                  child_array(kc,jcbc,ifl(ip):ifu(ip)) =                                            &
4230                       0.5_wp * ( child_array(kfu(kp-1),jcbc,ifl(ip):ifu(ip))                       &
4231                       + child_array(kfu(kp),jcbc,ifl(ip):ifu(ip)) )
4232               ENDDO
4233
4234            ENDDO
4235         ENDDO
4236
4237      ELSE   ! Any scalar
4238         
4239         DO  ip = iplw, iprw
4240            DO  kp = 0, kct
4241!
4242!--            Interpolate to the flux point
4243               c_interp_1 = cb * workarr_sn(kp,jpw,ip) + cp * workarr_sn(kp,jpwp,ip)
4244               DO  ic = ifl(ip), ifu(ip)
4245                  DO  kc = kfl(kp), kfu(kp)
4246                     child_array(kc,jcbc,ic) = c_interp_1
4247                  ENDDO
4248               ENDDO
4249
4250            ENDDO
4251         ENDDO
4252
4253      ENDIF  ! var
4254!
4255!--   Fill up also the redundant 2nd and 3rd ghost-node layers
4256      IF ( edge == 's' )  THEN
4257         DO  jcbgp = -nbgp, jcb
4258            child_array(0:nzt+1,jcbgp,nxlg:nxrg) = child_array(0:nzt+1,jcbc,nxlg:nxrg)
4259         ENDDO
4260      ELSEIF ( edge == 'n' )  THEN
4261         DO  jcbgp = jcb, ny+nbgp
4262            child_array(0:nzt+1,jcbgp,nxlg:nxrg) = child_array(0:nzt+1,jcbc,nxlg:nxrg)
4263         ENDDO
4264      ENDIF
4265
4266   END SUBROUTINE pmci_interp_1sto_sn
4267
4268
4269
4270   SUBROUTINE pmci_interp_1sto_t( child_array, parent_array, kct, ifl, ifu, jfl, jfu, var )
4271!
4272!--   Interpolation of ghost-node values used as the child-domain boundary
4273!--   conditions. This subroutine handles the top boundary.
4274      IMPLICIT NONE
4275
4276      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
4277     
4278      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
4279                                                              !< parent cell - x direction
4280      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
4281                                                              !< parent cell - x direction
4282      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
4283                                                              !< parent cell - y direction
4284      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
4285                                                              !< parent cell - y direction
4286
4287      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
4288      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
4289
4290      CHARACTER(LEN=1), INTENT(IN) ::  var                    !< Variable symbol: 'u', 'v', 'w' or 's'
4291!
4292!--   Local variables:     
4293      INTEGER(iwp) ::  ic          !< Running child-grid index in the x-direction
4294      INTEGER(iwp) ::  ierr        !< MPI error code
4295      INTEGER(iwp) ::  iplc        !< Lower parent-grid index limit in the x-direction for copying parent-grid
4296                                   !< array data to workarr_t
4297      INTEGER(iwp) ::  iprc        !< Upper parent-grid index limit in the x-direction for copying parent-grid
4298                                   !< array data to workarr_t
4299      INTEGER(iwp) ::  jc          !< Running child-grid index in the y-direction
4300      INTEGER(iwp) ::  jpsc        !< Lower parent-grid index limit in the y-direction for copying parent-grid
4301                                   !< array data to workarr_t
4302      INTEGER(iwp) ::  jpnc        !< Upper parent-grid-index limit in the y-direction for copying parent-grid
4303                                   !< array data to workarr_t
4304      INTEGER(iwp) ::  kc          !< Vertical child-grid index fixed to the boundary-value level
4305      INTEGER(iwp) ::  ip          !< Running parent-grid index in the x-direction
4306      INTEGER(iwp) ::  jp          !< Running parent-grid index in the y-direction
4307      INTEGER(iwp) ::  kpw         !< Reduced parent-grid index in the z-direction for workarr_t pointing to
4308                                   !< the boundary ghost node
4309      REAL(wp)     ::  c31         !< Interpolation coefficient for the 3rd-order WS scheme
4310      REAL(wp)     ::  c32         !< Interpolation coefficient for the 3rd-order WS scheme
4311      REAL(wp)     ::  c33         !< Interpolation coefficient for the 3rd-order WS scheme
4312      REAL(wp)     ::  c_interp_1  !< Value interpolated to the flux point in z direction from the parent-grid data
4313      REAL(wp)     ::  c_interp_2  !< Auxiliary value interpolated to the flux point in z direction from the parent-grid data
4314
4315
4316      IF ( var == 'w' )  THEN
4317         kc = nzt
4318      ELSE
4319         kc = nzt + 1
4320      ENDIF
4321      kpw = 1
4322!
4323!--   Interpolation coefficients
4324      IF ( interpolation_scheme_t == 1 )  THEN
4325         c31 =  0.0_wp           ! 1st-order upwind
4326         c32 =  1.0_wp
4327         c33 =  0.0_wp
4328      ELSE IF ( interpolation_scheme_t == 2 )  THEN
4329         c31 =  0.5_wp           ! 2nd-order central
4330         c32 =  0.5_wp
4331         c33 =  0.0_wp
4332      ELSE           
4333         c31 =  2.0_wp / 6.0_wp  ! 3rd-order WS upwind biased (default)
4334         c32 =  5.0_wp / 6.0_wp
4335         c33 = -1.0_wp / 6.0_wp         
4336      ENDIF         
4337!
4338!--   Substitute the necessary parent-grid data to the work array.
4339!--   Note that the dimension of workarr_t is (0:2,jpsw:jpnw,iplw:iprw),
4340!--   And the jc?w and ic?w-index bounds depend on the location of the PE-
4341!--   subdomain relative to the side boundaries.
4342      iplc = iplw + 1
4343      iprc = iprw - 1     
4344      jpsc = jpsw + 1
4345      jpnc = jpnw - 1
4346      IF ( bc_dirichlet_l )  THEN
4347         iplc = iplw
4348      ENDIF
4349      IF ( bc_dirichlet_r )  THEN
4350         iprc = iprw
4351      ENDIF
4352      IF ( bc_dirichlet_s )  THEN
4353         jpsc = jpsw
4354      ENDIF
4355      IF ( bc_dirichlet_n )  THEN
4356         jpnc = jpnw
4357      ENDIF
4358      workarr_t = 0.0_wp
4359      workarr_t(0:2,jpsc:jpnc,iplc:iprc) = parent_array(kct:kct+2,jpsc:jpnc,iplc:iprc)
4360!
4361!--   Left-right exchange if more than one PE subdomain in the x-direction.
4362!--   Note that in case of 3-D nesting the left and right boundaries are
4363!--   not exchanged because the nest domain is not cyclic.
4364      IF ( pdims(1) > 1 )  THEN
4365!
4366!--      From left to right
4367         CALL MPI_SENDRECV( workarr_t(0,jpsw,iplw+1), 1, workarr_t_exchange_type_y, pleft, 0,       &
4368                            workarr_t(0,jpsw,iprw), 1, workarr_t_exchange_type_y, pright, 0,        &
4369                            comm2d, status, ierr )
4370!                                                                             
4371!--      From right to left                                                   
4372         CALL MPI_SENDRECV( workarr_t(0,jpsw,iprw-1), 1, workarr_t_exchange_type_y, pright, 1,      &
4373                            workarr_t(0,jpsw,iplw), 1, workarr_t_exchange_type_y, pleft,  1,        &
4374                            comm2d, status, ierr )                                           
4375      ENDIF                                                                   
4376!                                                                             
4377!--   South-north exchange if more than one PE subdomain in the y-direction.   
4378!--   Note that in case of 3-D nesting the south and north boundaries are     
4379!--   not exchanged because the nest domain is not cyclic.                     
4380      IF ( pdims(2) > 1 )  THEN                                               
4381!                                                                             
4382!--      From south to north                                                   
4383         CALL MPI_SENDRECV( workarr_t(0,jpsw+1,iplw), 1, workarr_t_exchange_type_x, psouth, 2,      &
4384                            workarr_t(0,jpnw,iplw), 1, workarr_t_exchange_type_x, pnorth, 2,        &
4385                            comm2d, status, ierr )                                           
4386!                                                                             
4387!--      From north to south                                                   
4388         CALL MPI_SENDRECV( workarr_t(0,jpnw-1,iplw), 1, workarr_t_exchange_type_x, pnorth, 3,      &
4389                            workarr_t(0,jpsw,iplw), 1, workarr_t_exchange_type_x, psouth, 3,        &
4390                            comm2d, status, ierr )
4391      ENDIF
4392
4393      IF  ( var == 'w' )  THEN
4394         DO  ip = iplw, iprw
4395            DO  jp = jpsw, jpnw
4396 
4397               DO  ic = ifl(ip), ifu(ip)
4398                  DO  jc = jfl(jp), jfu(jp)
4399                     child_array(kc,jc,ic) = workarr_t(kpw,jp,ip)
4400                  ENDDO
4401               ENDDO
4402
4403            ENDDO
4404         ENDDO
4405
4406      ELSE IF  ( var == 'u' )  THEN
4407
4408         DO  ip = iplw, iprw - 1
4409            DO  jp = jpsw, jpnw
4410!
4411!--            First interpolate to the flux point using the 3rd-order WS scheme
4412               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip)   + c32 * workarr_t(kpw,jp,ip)             &
4413                          + c33 * workarr_t(kpw+1,jp,ip)
4414               c_interp_2 = c31 * workarr_t(kpw-1,jp,ip+1) + c32 * workarr_t(kpw,jp,ip+1)           &
4415                          + c33 * workarr_t(kpw+1,jp,ip+1)
4416!
4417!--            Use averages of the neighbouring matching grid-line values
4418               DO  ic = ifl(ip), ifl(ip+1)
4419                  child_array(kc,jfl(jp):jfu(jp),ic) = 0.5_wp * ( c_interp_1 + c_interp_2 )
4420               ENDDO
4421!
4422!--            Then set the values along the matching grid-lines 
4423               IF ( MOD( ifl(ip), igsr ) == 0 )  THEN
4424!
4425!--               First interpolate to the flux point using the 3rd-order WS scheme
4426                  c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)            &
4427                             + c33 * workarr_t(kpw+1,jp,ip)                 
4428                  child_array(kc,jfl(jp):jfu(jp),ifl(ip)) = c_interp_1
4429               ENDIF
4430
4431            ENDDO
4432         ENDDO
4433!
4434!--      Finally, set the values along the last matching grid-line 
4435         IF  ( MOD( ifl(iprw), igsr ) == 0 )  THEN
4436            DO  jp = jpsw, jpnw
4437!
4438!--            First interpolate to the flux point using the 3rd-order WS scheme
4439               c_interp_1 = c31 * workarr_t(kpw-1,jp,iprw) + c32 * workarr_t(kpw,jp,iprw)           &
4440                          + c33 * workarr_t(kpw+1,jp,iprw)
4441               child_array(kc,jfl(jp):jfu(jp),ifl(iprw)) = c_interp_1
4442            ENDDO
4443         ENDIF
4444!
4445!--      A gap may still remain in some cases if the subdomain size is not
4446!--      divisible by the grid-spacing ratio. In such a case, fill the
4447!--      gap. Note however, this operation may produce some additional
4448!--      momentum conservation error.
4449         IF ( ifl(iprw) < nxr )  THEN
4450            DO  jp = jpsw, jpnw
4451               DO  ic = ifl(iprw) + 1, nxr
4452                  child_array(kc,jfl(jp):jfu(jp),ic) = child_array(kc,jfl(jp):jfu(jp),ifl(iprw))
4453               ENDDO
4454            ENDDO
4455         ENDIF
4456
4457      ELSE IF  ( var == 'v' )  THEN
4458
4459         DO  ip = iplw, iprw
4460            DO  jp = jpsw, jpnw-1
4461!
4462!--            First interpolate to the flux point using the 3rd-order WS scheme
4463               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip)   + c32 * workarr_t(kpw,jp,ip)             &
4464                          + c33 * workarr_t(kpw+1,jp,ip)
4465               c_interp_2 = c31 * workarr_t(kpw-1,jp+1,ip) + c32 * workarr_t(kpw,jp+1,ip)           &
4466                          + c33 * workarr_t(kpw+1,jp+1,ip)
4467!
4468!--            Use averages of the neighbouring matching grid-line values
4469               DO  jc = jfl(jp), jfl(jp+1)         
4470                  child_array(kc,jc,ifl(ip):ifu(ip)) = 0.5_wp * ( c_interp_1 + c_interp_2 )
4471               ENDDO
4472!
4473!--            Then set the values along the matching grid-lines 
4474               IF ( MOD( jfl(jp), jgsr ) == 0 )  THEN
4475                  c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)            &
4476                             + c33 * workarr_t(kpw+1,jp,ip)
4477                  child_array(kc,jfl(jp),ifl(ip):ifu(ip)) = c_interp_1
4478               ENDIF
4479               
4480            ENDDO
4481
4482         ENDDO
4483!
4484!--      Finally, set the values along the last matching grid-line
4485         IF ( MOD( jfl(jpnw), jgsr ) == 0 )  THEN
4486            DO  ip = iplw, iprw
4487!
4488!--            First interpolate to the flux point using the 3rd-order WS scheme
4489               c_interp_1 = c31 * workarr_t(kpw-1,jpnw,ip) + c32 * workarr_t(kpw,jpnw,ip)           &
4490                          + c33 * workarr_t(kpw+1,jpnw,ip)
4491               child_array(kc,jfl(jpnw),ifl(ip):ifu(ip)) = c_interp_1
4492            ENDDO
4493         ENDIF
4494!
4495!--      A gap may still remain in some cases if the subdomain size is not
4496!--      divisible by the grid-spacing ratio. In such a case, fill the
4497!--      gap. Note however, this operation may produce some additional
4498!--      momentum conservation error.
4499         IF  ( jfl(jpnw) < nyn )  THEN
4500            DO  ip = iplw, iprw
4501               DO  jc = jfl(jpnw)+1, nyn
4502                  child_array(kc,jc,ifl(ip):ifu(ip)) = child_array(kc,jfl(jpnw),ifl(ip):ifu(ip))
4503               ENDDO
4504            ENDDO
4505         ENDIF
4506
4507      ELSE  ! any scalar variable
4508
4509         DO  ip = iplw, iprw
4510            DO  jp = jpsw, jpnw
4511!
4512!--            First interpolate to the flux point using the 3rd-order WS scheme
4513               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)               &
4514                          + c33 * workarr_t(kpw+1,jp,ip)
4515               DO  ic = ifl(ip), ifu(ip)
4516                  DO  jc = jfl(jp), jfu(jp)
4517                     child_array(kc,jc,ic) = c_interp_1
4518                  ENDDO
4519               ENDDO
4520
4521            ENDDO
4522         ENDDO
4523
4524      ENDIF  ! var
4525!
4526!--   Just fill up the redundant second ghost-node layer in case of var == w.
4527      IF ( var == 'w' )  THEN
4528         child_array(nzt+1,:,:) = child_array(nzt,:,:)
4529      ENDIF
4530
4531   END SUBROUTINE pmci_interp_1sto_t
4532
4533
4534
4535   SUBROUTINE pmci_anterp_tophat( child_array, parent_array, kct, ifl, ifu, jfl, jfu, kfl, kfu,     &
4536                                  ijkfc, var )
4537!
4538!--   Anterpolation of internal-node values to be used as the parent-domain
4539!--   values. This subroutine is based on the first-order numerical
4540!--   integration of the child-grid values contained within the anterpolation
4541!--   cell (Clark & Farley, Journal of the Atmospheric Sciences 41(3), 1984).
4542
4543      IMPLICIT NONE
4544
4545      INTEGER(iwp), INTENT(IN) ::  kct  !< Top boundary index for anterpolation along z
4546     
4547      INTEGER(iwp), DIMENSION(0:pg%nz+1,jpsa:jpna,ipla:ipra), INTENT(IN) ::  ijkfc  !< number of child grid points contributing
4548                                                                                    !< to a parent grid box
4549      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
4550                                                              !< parent cell - x direction
4551      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
4552                                                              !< parent cell - x direction
4553      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
4554                                                              !< parent cell - y direction
4555      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
4556                                                              !< parent cell - y direction
4557      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
4558                                                              !< parent cell - z direction
4559      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
4560                                                              !< parent cell - z direction
4561
4562      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  child_array   !< Child-grid array
4563      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(INOUT)  ::  parent_array  !< Parent-grid array
4564
4565      CHARACTER(LEN=*), INTENT(IN) ::  var                   !< Variable symbol: 'u', 'v', 'w' or 's'
4566!
4567!--   Local variables: 
4568      INTEGER(iwp) ::  ic              !< Running index x-direction - child grid
4569      INTEGER(iwp) ::  ipl_anterp      !< Left boundary index for anterpolation along x
4570      INTEGER(iwp) ::  ipr_anterp      !< Right boundary index for anterpolation along x
4571      INTEGER(iwp) ::  jc              !< Running index y-direction - child grid
4572      INTEGER(iwp) ::  jpn_anterp      !< North boundary index for anterpolation along y
4573      INTEGER(iwp) ::  jps_anterp      !< South boundary index for anterpolation along y
4574      INTEGER(iwp) ::  kc              !< Running index z-direction - child grid     
4575      INTEGER(iwp) ::  kpb_anterp = 0  !< Bottom boundary index for anterpolation along z
4576      INTEGER(iwp) ::  kpt_anterp      !< Top boundary index for anterpolation along z
4577      INTEGER(iwp) ::  ip              !< Running index x-direction - parent grid
4578      INTEGER(iwp) ::  jp              !< Running index y-direction - parent grid
4579      INTEGER(iwp) ::  kp              !< Running index z-direction - parent grid
4580      INTEGER(iwp) ::  var_flag        !< bit number used to flag topography on respective grid
4581
4582      REAL(wp) ::  cellsum       !< sum of respective child cells belonging to parent cell
4583
4584!
4585!--   Define the index bounds ipl_anterp, ipr_anterp, jps_anterp and jpn_anterp.
4586!--   Note that kcb_anterp is simply zero and kct_anterp depends on kct which enters
4587!--   here as a parameter and it is determined in pmci_define_index_mapping.
4588!--   Note that the grid points directly used also for interpolation (from parent to
4589!--   child) are always excluded from anterpolation, e.g. anterpolation is maximally
4590!--   only from 0:kct-1, since kct is directly used for interpolation. Similar restriction is
4591!--   applied to the lateral boundaries as well. An additional buffer is
4592!--   also applied (default value for anterpolation_buffer_width = 2) in order
4593!--   to avoid unphysical accumulation of kinetic energy.
4594      ipl_anterp = ipl
4595      ipr_anterp = ipr
4596      jps_anterp = jps
4597      jpn_anterp = jpn
4598      kpb_anterp = 0
4599      kpt_anterp = kct - 1 - anterpolation_buffer_width
4600
4601      IF ( nesting_mode /= 'vertical' )  THEN
4602!
4603!--      Set the anterpolation buffers on the lateral boundaries
4604         ipl_anterp = MAX( ipl, iplg + 3 + anterpolation_buffer_width )
4605         ipr_anterp = MIN( ipr, iprg - 3 - anterpolation_buffer_width )
4606         jps_anterp = MAX( jps, jpsg + 3 + anterpolation_buffer_width )
4607         jpn_anterp = MIN( jpn, jpng - 3 - anterpolation_buffer_width )
4608         
4609      ENDIF
4610!
4611!--   Set masking bit for topography flags
4612      IF ( var == 'u' )  THEN
4613         var_flag = 1 
4614      ELSEIF ( var == 'v' )  THEN
4615         var_flag = 2 
4616      ELSEIF ( var == 'w' )  THEN
4617         var_flag = 3
4618      ELSE
4619         var_flag = 0
4620      ENDIF
4621!
4622!--   Note that ip, jp, and kp are parent-grid indices and ic,jc, and kc
4623!--   are child-grid indices.
4624      DO  ip = ipl_anterp, ipr_anterp
4625         DO  jp = jps_anterp, jpn_anterp
4626!
4627!--         For simplicity anterpolate within buildings and under elevated
4628!--         terrain too
4629            DO  kp = kpb_anterp, kpt_anterp
4630               cellsum = 0.0_wp
4631               DO  ic = ifl(ip), ifu(ip)
4632                  DO  jc = jfl(jp), jfu(jp)
4633                     DO  kc = kfl(kp), kfu(kp)
4634                        cellsum = cellsum + MERGE( child_array(kc,jc,ic), 0.0_wp,                   &
4635                             BTEST( wall_flags_static_0(kc,jc,ic), var_flag ) )
4636                     ENDDO
4637                  ENDDO
4638               ENDDO
4639!
4640!--            In case all child grid points are inside topography, i.e.
4641!--            ijkfc and cellsum are zero, also parent solution would have
4642!--            zero values at that grid point, which may cause problems in
4643!--            particular for the temperature. Therefore, in case cellsum is
4644!--            zero, keep the parent solution at this point.
4645               IF ( ijkfc(kp,jp,ip) /= 0 )  THEN
4646                  parent_array(kp,jp,ip) = cellsum / REAL( ijkfc(kp,jp,ip), KIND=wp )
4647               ENDIF
4648
4649            ENDDO
4650         ENDDO
4651      ENDDO
4652
4653   END SUBROUTINE pmci_anterp_tophat
4654
4655#endif
4656
4657 END SUBROUTINE pmci_child_datatrans
4658
4659! Description:
4660! ------------
4661!> Set boundary conditions for the prognostic quantities after interpolation
4662!> and anterpolation at upward- and downward facing surfaces. 
4663!> @todo: add Dirichlet boundary conditions for pot. temperature, humdidity and
4664!> passive scalar.
4665!------------------------------------------------------------------------------!
4666 SUBROUTINE pmci_boundary_conds
4667
4668#if defined( __parallel )
4669    USE chem_modules,                                                          &
4670        ONLY:  ibc_cs_b
4671
4672    USE control_parameters,                                                    &
4673        ONLY:  ibc_pt_b, ibc_q_b, ibc_s_b, ibc_uv_b
4674
4675    USE salsa_mod,                                                             &
4676        ONLY:  ibc_salsa_b
4677
4678    USE surface_mod,                                                           &
4679        ONLY:  bc_h
4680
4681    IMPLICIT NONE
4682
4683    INTEGER(iwp) ::  ic  !< Index along x-direction
4684    INTEGER(iwp) ::  jc  !< Index along y-direction
4685    INTEGER(iwp) ::  kc  !< Index along z-direction
4686    INTEGER(iwp) ::  lb  !< Running index for aerosol size bins
4687    INTEGER(iwp) ::  lc  !< Running index for aerosol mass bins
4688    INTEGER(iwp) ::  lg  !< Running index for salsa gases
4689    INTEGER(iwp) ::  m   !< Running index for surface type
4690    INTEGER(iwp) ::  n   !< Running index for number of chemical species
4691   
4692
4693    IF ( debug_output_timestep )  CALL debug_message( 'pmci_boundary_conds', 'start' )
4694!
4695!-- Set Dirichlet boundary conditions for horizontal velocity components
4696    IF ( ibc_uv_b == 0 )  THEN
4697!
4698!--    Upward-facing surfaces
4699       DO  m = 1, bc_h(0)%ns
4700          ic = bc_h(0)%i(m)           
4701          jc = bc_h(0)%j(m)
4702          kc = bc_h(0)%k(m)
4703          u(kc-1,jc,ic) = 0.0_wp
4704          v(kc-1,jc,ic) = 0.0_wp
4705       ENDDO
4706!
4707!--    Downward-facing surfaces
4708       DO  m = 1, bc_h(1)%ns
4709          ic = bc_h(1)%i(m)           
4710          jc = bc_h(1)%j(m)
4711          kc = bc_h(1)%k(m)
4712          u(kc+1,jc,ic) = 0.0_wp
4713          v(kc+1,jc,ic) = 0.0_wp
4714       ENDDO
4715    ENDIF
4716!
4717!-- Set Dirichlet boundary conditions for vertical velocity component
4718!-- Upward-facing surfaces
4719    DO  m = 1, bc_h(0)%ns
4720       ic = bc_h(0)%i(m)           
4721       jc = bc_h(0)%j(m)
4722       kc = bc_h(0)%k(m)
4723       w(kc-1,jc,ic) = 0.0_wp
4724    ENDDO
4725!
4726!-- Downward-facing surfaces
4727    DO  m = 1, bc_h(1)%ns
4728       ic = bc_h(1)%i(m)           
4729       jc = bc_h(1)%j(m)
4730       kc = bc_h(1)%k(m)
4731       w(kc+1,jc,ic) = 0.0_wp
4732    ENDDO
4733!
4734!-- Set Neumann boundary conditions for potential temperature
4735    IF ( .NOT. neutral )  THEN
4736       IF ( ibc_pt_b == 1 )  THEN
4737          DO  m = 1, bc_h(0)%ns
4738             ic = bc_h(0)%i(m)           
4739             jc = bc_h(0)%j(m)
4740             kc = bc_h(0)%k(m)
4741             pt(kc-1,jc,ic) = pt(kc,jc,ic)
4742          ENDDO
4743          DO  m = 1, bc_h(1)%ns
4744             ic = bc_h(1)%i(m)           
4745             jc = bc_h(1)%j(m)
4746             kc = bc_h(1)%k(m)
4747             pt(kc+1,jc,ic) = pt(kc,jc,ic)
4748          ENDDO   
4749       ENDIF
4750    ENDIF
4751!
4752!-- Set Neumann boundary conditions for humidity and cloud-physical quantities
4753    IF ( humidity )  THEN
4754       IF ( ibc_q_b == 1 )  THEN
4755          DO  m = 1, bc_h(0)%ns
4756             ic = bc_h(0)%i(m)           
4757             jc = bc_h(0)%j(m)
4758             kc = bc_h(0)%k(m)
4759             q(kc-1,jc,ic) = q(kc,jc,ic)
4760          ENDDO 
4761          DO  m = 1, bc_h(1)%ns
4762             ic = bc_h(1)%i(m)           
4763             jc = bc_h(1)%j(m)
4764             kc = bc_h(1)%k(m)
4765             q(kc+1,jc,ic) = q(kc,jc,ic)
4766          ENDDO 
4767       ENDIF
4768       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
4769          DO  m = 1, bc_h(0)%ns
4770             ic = bc_h(0)%i(m)           
4771             jc = bc_h(0)%j(m)
4772             kc = bc_h(0)%k(m)
4773             nc(kc-1,jc,ic) = 0.0_wp
4774             qc(kc-1,jc,ic) = 0.0_wp
4775          ENDDO 
4776          DO  m = 1, bc_h(1)%ns
4777             ic = bc_h(1)%i(m)           
4778             jc = bc_h(1)%j(m)
4779             kc = bc_h(1)%k(m)
4780
4781             nc(kc+1,jc,ic) = 0.0_wp
4782             qc(kc+1,jc,ic) = 0.0_wp
4783          ENDDO 
4784       ENDIF
4785
4786       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
4787          DO  m = 1, bc_h(0)%ns
4788             ic = bc_h(0)%i(m)           
4789             jc = bc_h(0)%j(m)
4790             kc = bc_h(0)%k(m)
4791             nr(kc-1,jc,ic) = 0.0_wp
4792             qr(kc-1,jc,ic) = 0.0_wp
4793          ENDDO 
4794          DO  m = 1, bc_h(1)%ns
4795             ic = bc_h(1)%i(m)           
4796             jc = bc_h(1)%j(m)
4797             kc = bc_h(1)%k(m)
4798             nr(kc+1,jc,ic) = 0.0_wp
4799             qr(kc+1,jc,ic) = 0.0_wp
4800          ENDDO 
4801       ENDIF
4802
4803    ENDIF
4804!
4805!-- Set Neumann boundary conditions for passive scalar
4806    IF ( passive_scalar )  THEN
4807       IF ( ibc_s_b == 1 )  THEN
4808          DO  m = 1, bc_h(0)%ns
4809             ic = bc_h(0)%i(m)           
4810             jc = bc_h(0)%j(m)
4811             kc = bc_h(0)%k(m)
4812             s(kc-1,jc,ic) = s(kc,jc,ic)
4813          ENDDO
4814          DO  m = 1, bc_h(1)%ns
4815             ic = bc_h(1)%i(m)           
4816             jc = bc_h(1)%j(m)
4817             kc = bc_h(1)%k(m)
4818             s(kc+1,jc,ic) = s(kc,jc,ic)
4819          ENDDO 
4820       ENDIF
4821    ENDIF
4822!
4823!-- Set Neumann boundary conditions for chemical species
4824    IF ( air_chemistry  .AND.  nesting_chem )  THEN
4825       IF ( ibc_cs_b == 1 )  THEN
4826          DO  n = 1, nspec
4827             DO  m = 1, bc_h(0)%ns
4828                ic = bc_h(0)%i(m)           
4829                jc = bc_h(0)%j(m)
4830                kc = bc_h(0)%k(m)
4831                chem_species(n)%conc(kc-1,jc,ic) = chem_species(n)%conc(kc,jc,ic)
4832             ENDDO
4833             DO  m = 1, bc_h(1)%ns
4834                ic = bc_h(1)%i(m)           
4835                jc = bc_h(1)%j(m)
4836                kc = bc_h(1)%k(m)
4837                chem_species(n)%conc(kc+1,jc,ic) = chem_species(n)%conc(kc,jc,ic)
4838             ENDDO
4839          ENDDO
4840       ENDIF
4841    ENDIF
4842!
4843!-- Set Neumann boundary conditions for aerosols and salsa gases
4844    IF ( salsa  .AND.  nesting_salsa )  THEN
4845       IF ( ibc_salsa_b == 1 )  THEN
4846          DO  m = 1, bc_h(0)%ns
4847             ic = bc_h(0)%i(m)
4848             jc = bc_h(0)%j(m)
4849             kc = bc_h(0)%k(m)
4850             DO  lb = 1, nbins_aerosol
4851                aerosol_number(lb)%conc(kc-1,jc,ic) = aerosol_number(lb)%conc(kc,jc,ic)
4852             ENDDO
4853             DO  lc = 1, nbins_aerosol * ncomponents_mass
4854                aerosol_mass(lc)%conc(kc-1,jc,ic) = aerosol_mass(lc)%conc(kc,jc,ic)
4855             ENDDO
4856             IF ( .NOT. salsa_gases_from_chem )  THEN
4857                DO  lg = 1, ngases_salsa
4858                   salsa_gas(lg)%conc(kc-1,jc,ic) = salsa_gas(lg)%conc(kc,jc,ic)
4859                ENDDO
4860             ENDIF
4861          ENDDO
4862          DO  m = 1, bc_h(1)%ns
4863             ic = bc_h(1)%i(m)
4864             jc = bc_h(1)%j(m)
4865             kc = bc_h(1)%k(m)
4866             DO  lb = 1, nbins_aerosol
4867                aerosol_number(lb)%conc(kc+1,jc,ic) = aerosol_number(lb)%conc(kc,jc,ic)
4868             ENDDO
4869             DO  lc = 1, nbins_aerosol * ncomponents_mass
4870                aerosol_mass(lc)%conc(kc+1,jc,ic) = aerosol_mass(lc)%conc(kc,jc,ic)
4871             ENDDO
4872             IF ( .NOT. salsa_gases_from_chem )  THEN
4873                DO  lg = 1, ngases_salsa
4874                   salsa_gas(lg)%conc(kc+1,jc,ic) = salsa_gas(lg)%conc(kc,jc,ic)
4875                ENDDO
4876             ENDIF
4877          ENDDO
4878       ENDIF
4879    ENDIF   
4880
4881    IF ( debug_output_timestep )  CALL debug_message( 'pmci_boundary_conds', 'end' )
4882
4883#endif
4884 END SUBROUTINE pmci_boundary_conds
4885
4886
4887 
4888 SUBROUTINE pmci_ensure_nest_mass_conservation
4889
4890!
4891!-- Adjust the volume-flow rate through the nested boundaries so that the net volume
4892!-- flow through all boundaries of the current nest domain becomes zero.
4893    IMPLICIT NONE
4894
4895    INTEGER(iwp) ::  i                        !< Running index in the x-direction
4896    INTEGER(iwp) ::  ierr                     !< MPI error code
4897    INTEGER(iwp) ::  j                        !< Running index in the y-direction
4898    INTEGER(iwp) ::  k                        !< Running index in the z-direction
4899    INTEGER(iwp) ::  n                        !< Running index over the boundary faces: l, r, s, n and t
4900
4901    REAL(wp) ::  dxdy                         !< Surface area of grid cell top face
4902    REAL(wp) ::  innor                        !< Inner normal vector of the grid cell face
4903    REAL(wp) ::  sub_sum                      !< Intermediate sum for reducing the loss of signifigant digits in 2-D summations
4904    REAL(wp) ::  u_corr_left                  !< Correction added to the left boundary value of u
4905    REAL(wp) ::  u_corr_right                 !< Correction added to the right boundary value of u
4906    REAL(wp) ::  v_corr_south                 !< Correction added to the south boundary value of v
4907    REAL(wp) ::  v_corr_north                 !< Correction added to the north boundary value of v
4908    REAL(wp) ::  volume_flux_integral         !< Surface integral of volume flux over the domain boundaries
4909    REAL(wp) ::  volume_flux_local            !< Surface integral of volume flux over the subdomain boundary face
4910    REAL(wp) ::  w_corr_top                   !< Correction added to the top boundary value of w
4911
4912    REAL(wp), DIMENSION(5) ::  volume_flux    !< Surface integral of volume flux over each boundary face of the domain
4913
4914   
4915!
4916!-- Sum up the volume flow through the left boundary
4917    volume_flux(1) = 0.0_wp
4918    volume_flux_local = 0.0_wp
4919    IF ( bc_dirichlet_l )  THEN
4920       i = 0
4921       innor = dy
4922       DO   j = nys, nyn
4923          sub_sum = 0.0_wp
4924          DO   k = nzb+1, nzt
4925             sub_sum = sub_sum + innor * u(k,j,i) * dzw(k)                                          &
4926                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 1 ) )
4927          ENDDO
4928          volume_flux_local = volume_flux_local + sub_sum
4929       ENDDO
4930    ENDIF
4931
4932#if defined( __parallel )
4933    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
4934    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(1), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
4935#else
4936    volume_flux(1) = volume_flux_local
4937#endif
4938!
4939!-- Sum up the volume flow through the right boundary
4940    volume_flux(2) = 0.0_wp
4941    volume_flux_local = 0.0_wp
4942    IF ( bc_dirichlet_r )  THEN
4943       i = nx + 1
4944       innor = -dy
4945       DO   j = nys, nyn
4946          sub_sum = 0.0_wp
4947          DO   k = nzb+1, nzt
4948             sub_sum = sub_sum + innor * u(k,j,i) * dzw(k)                                          &
4949                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 1 ) )
4950          ENDDO
4951          volume_flux_local = volume_flux_local + sub_sum
4952       ENDDO
4953    ENDIF
4954
4955#if defined( __parallel )
4956    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
4957    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(2), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
4958#else
4959    volume_flux(2) = volume_flux_local
4960#endif
4961!
4962!-- Sum up the volume flow through the south boundary
4963    volume_flux(3) = 0.0_wp   
4964    volume_flux_local = 0.0_wp
4965    IF ( bc_dirichlet_s )  THEN
4966       j = 0
4967       innor = dx
4968       DO   i = nxl, nxr
4969          sub_sum = 0.0_wp
4970          DO   k = nzb+1, nzt
4971             sub_sum = sub_sum + innor * v(k,j,i) * dzw(k)                                          &
4972                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 2 ) )
4973          ENDDO
4974          volume_flux_local = volume_flux_local + sub_sum
4975       ENDDO
4976    ENDIF
4977   
4978#if defined( __parallel )
4979    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
4980    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(3), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
4981#else
4982    volume_flux(3) = volume_flux_local
4983#endif
4984!   
4985!-- Sum up the volume flow through the north boundary
4986    volume_flux(4) = 0.0_wp
4987    volume_flux_local = 0.0_wp
4988    IF ( bc_dirichlet_n )  THEN
4989       j = ny + 1
4990       innor = -dx
4991       DO  i = nxl, nxr
4992          sub_sum = 0.0_wp
4993          DO  k = nzb+1, nzt
4994             sub_sum = sub_sum + innor * v(k,j,i) * dzw(k)                                          &
4995                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 2 ) )
4996          ENDDO
4997          volume_flux_local = volume_flux_local + sub_sum
4998       ENDDO
4999    ENDIF
5000
5001#if defined( __parallel )
5002    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
5003    CALL MPI_ALLREDUCE( volume_flux_local, volume_flux(4), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
5004#else
5005    volume_flux(4) = volume_flux_local
5006#endif
5007!
5008!-- Sum up the volume flow through the top boundary
5009</