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

Last change on this file since 3592 was 3592, checked in by suehring, 5 years ago

Determine number of coupled arrays dynamically

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