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

Last change on this file since 2967 was 2967, checked in by raasch, 3 years ago

bugfix: missing parallel cpp-directives added

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