Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |