Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |