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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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