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

Last change on this file since 3804 was 3804, checked in by hellstea, 3 years ago

Last commit documented

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