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

Last change on this file since 2803 was 2803, checked in by thiele, 4 years ago

preprocessor directive for c_sizeof in case of gfortran added

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