Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |