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

Last change on this file since 3824 was 3822, checked in by hellstea, 5 years ago

Temporary increase of the vertical dimension of the parent-grid arrays and workarrc_t is cancelled as unnecessary

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