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

Last change on this file since 2903 was 2903, checked in by hellstea, 4 years ago

Nesting-related calls to pmci_ensure_nest_mass_conservation and pres after the nest initialization are removed as they may create unwanted initial perturbation in some cases

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