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

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

A bug fixed in lateral boundary interpolations

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