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

Last change on this file since 4249 was 4249, checked in by hellstea, 22 months ago

Several grid-line matching tests changed to a round-off-error tolerant form

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