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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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