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

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

unused variables removed, missing working precision added, missing preprocessor directives added, bugfix concerning allocation of t_surf_wall_v in nopointer case, declaration statements rearranged to avoid compile time errors, mpi_abort arguments replaced to avoid compile errors

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