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

Last change on this file since 4260 was 4260, checked in by hellstea, 2 years ago

Rest of the grid-line matching tests in pmc_interface_mod changed to round-off-error tolerant forms

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