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

Last change on this file since 4385 was 4385, checked in by hellstea, 17 months ago

Nesting-related error messages PA0425 and PA0426 made more specific

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