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

Last change on this file since 3182 was 3182, checked in by suehring, 4 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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