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

Last change on this file since 2020 was 2020, checked in by hellstea, 7 years ago

last commit documented

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