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

Last change on this file since 3794 was 3794, checked in by raasch, 3 years ago

unused variables removed

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