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

Last change on this file since 2795 was 2795, checked in by hellstea, 5 years ago

Bugfix in nesting anterpolation relaxation functions

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