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

Last change on this file since 3945 was 3945, checked in by raasch, 5 years ago

messed document changes for r3932 cleaned up

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