Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |