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

Last change on this file since 2311 was 2311, checked in by suehring, 7 years ago

Set bottom boundary conditions after nesting interpolation and anterpolation

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