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

Last change on this file since 2809 was 2809, checked in by schwenkel, 6 years ago

Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE

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