Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |