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

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

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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