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

Last change on this file since 3183 was 3183, checked in by suehring, 3 years ago

last commit documented

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