Home

Context Navigation

source: palm/trunk/SOURCE/pmc_interface.f90 @ 1798

Last change on this file since 1798 was 1798, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 152.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |