Home

Context Navigation

source: palm/trunk/SOURCE/pmc_interface_mod.f90 @ 1926

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |