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

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

Adjustable anterpolation buffer introduced on all nest boundaries

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