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

Last change on this file since 3708 was 3708, checked in by hellstea, 3 years ago

Checks for parent / child grid line matching introduced

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