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

Last change on this file since 2350 was 2350, checked in by kanani, 4 years ago

bugfixes and workarounds for nopointer version

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