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

Last change on this file since 2841 was 2841, checked in by knoop, 4 years ago

Bugfix: wrong placement of include 'mpif.h' corrected

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