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

Last change on this file since 3943 was 3943, checked in by maronga, 2 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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