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

Last change on this file since 4026 was 4026, checked in by suehring, 2 years ago

Improve albedo initialization in land-surface model by setting always some default albedo types, in order to be independent on any user settings; Mask topography at boundary grid cells in nesting mass conservation

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