Home

Context Navigation

pmc_interface.f90 @ 1791

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

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

Download in other formats:

�