Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |