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

Last change on this file since 4598 was 4508, checked in by raasch, 12 months ago

salsa decycling replaced by explicit setting of lateral boundary conditions

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