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

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

pmc bugfix and gfortran workaround

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