Home

Context Navigation

pmc_interface.f90 @ 1834

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

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

Context Navigation

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

Download in other formats: