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

Last change on this file since 2701 was 2701, checked in by suehring, 5 years ago

changes from last commit documented

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