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

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

Consider multiple pure-vertical nest domains in overlap check

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