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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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