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

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

Bugfixes in computation of the interpolation loglaw-correction parameters

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