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

Last change on this file since 2669 was 2669, checked in by raasch, 4 years ago

file attributes and activation strings in .palm.iofiles revised, file extensions for nesting, masked output, wind turbine data, etc. changed

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