Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |