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

Last change on this file since 3876 was 3876, checked in by knoop, 2 years ago

Moved "photolysis_scheme", "chem_species" and "phot_frequen" to chem_modules

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