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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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