Home

Context Navigation

source: palm/trunk/SOURCE/pmc_interface.f90 @ 1805

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |