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

Last change on this file since 4010 was 4010, checked in by hellstea, 2 years ago

Nest mass conservation correction included

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