Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |