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

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

last commit documented

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