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

Last change on this file since 4346 was 4346, checked in by motisi, 23 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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