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

Last change on this file since 2318 was 2318, checked in by suehring, 4 years ago

get topograpyhy top index via function call

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