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

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

In anterpolation, exclude grid points used for interpolation from parent to child

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