Home

Context Navigation

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

Last change on this file since 1764 was 1764, checked in by raasch, 8 years ago
update of the nested domain system + some bugfixes
Property svn:keywords set to `Id`
File size: 144.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |