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

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

Major update of pmc_interface_mod

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