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

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

Checks and error messages improved and extended. Number of variables renamed etc

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