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

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

Add missing if statements for call of pmc_set_dataarray_name for TKE and dissipation; minor bugfix for nesting of chemical species

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