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

Last change on this file since 2801 was 2801, checked in by thiele, 5 years ago

Introduce particle transfer in nested models

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