Home

Context Navigation

pmc_interface_mod.f90 @ 2003

Last change on this file since 2003 was 2003, checked in by suehring, 8 years ago

Humidity and passive scalar also separated in nesting mode

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/trunk/SOURCE/pmc_interface_mod.f90	merged	eligible
/palm/branches/fricke/SOURCE/pmc_interface_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90	296-409
/palm/branches/suehring/pmc_interface_mod.f90	423-666

File size: 161.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: