Home

Context Navigation

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

Last change on this file since 1891 was 1883, checked in by hellstea, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 151.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |