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

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

bugfix in child initialization in pmc_interface_mod

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