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

Last change on this file since 2804 was 2804, checked in by thiele, 4 years ago

preprocessor directive for c_sizeof in case of gfortran added

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