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

Last change on this file since 4029 was 4029, checked in by raasch, 2 years ago

bugfix: decycling of chemistry species after nesting data transfer, exchange of ghost points and boundary conditions separated for chemical species and SALSA module, nest_chemistry option removed, netcdf variable NF90_NOFILL is used as argument instead of 1 in calls to NF90_DEF_VAR_FILL

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