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

Last change on this file since 2292 was 2292, checked in by schwenkel, 7 years ago

implementation of new bulk microphysics scheme

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