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

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

New checks added for nested setups

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