Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |