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

Last change on this file since 3833 was 3833, checked in by forkel, 5 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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