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

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

Child-domain size checks extended

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