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

Last change on this file since 3484 was 3484, checked in by hellstea, 4 years ago

Nesting interpolation made mass-conservative

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