Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |