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

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

some cleaning up

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