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

Last change on this file since 2174 was 2174, checked in by maronga, 5 years ago

nesting extended for cloud physics quantities, bugfixes in nesting, change of default most_method

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