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

Last change on this file since 3976 was 3976, checked in by hellstea, 6 years ago

Child initialization extended to the redundant ghost points behind the nest boundaries

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