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

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

Bugfix, wrong coarse-grid index in init_tkefactor used.

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