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

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

restructure/add location/debug messages

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