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

Last change on this file since 2317 was 2317, checked in by suehring, 8 years ago

get topograpyh top index via function call

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