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

Last change on this file since 4273 was 4273, checked in by monakurppa, 6 years ago

Add logical switched nesting_chem and nesting_offline_chem

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