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

Last change on this file since 3888 was 3888, checked in by hellstea, 3 years ago

bugfix in pmc_handle_communicator_mod

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