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

Last change on this file since 3987 was 3987, checked in by kanani, 2 years ago

clean up location, debug and error messages

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