Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |