Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |