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

Last change on this file since 4413 was 4413, checked in by hellstea, 16 months ago

All the USE-statements within subroutines moved to the module declaration section

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