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

Last change on this file since 2853 was 2853, checked in by suehring, 7 years ago

Bugfix in init log-law correction

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