Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |