Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |