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

Last change on this file since 3792 was 3792, checked in by hellstea, 6 years ago

Interpolations improved. Large number of obsolete subroutines removed. All unused variables removed.

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