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

Last change on this file since 3875 was 3864, checked in by monakurppa, 5 years ago

major changes in salsa: data input, format and performance

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