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

Last change on this file since 3083 was 3083, checked in by gronemeier, 7 years ago

merge with branch rans

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