Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |