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

Last change on this file since 2868 was 2868, checked in by hellstea, 7 years ago

Local conditional Neumann conditions for one-way coupling removed

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