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

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

Adjustments according new topography and surface-modelling concept implemented

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