Home

Context Navigation

pmc_interface_mod.f90 @ 2230

Last change on this file since 2230 was 2230, checked in by hellstea, 8 years ago

last commit documented

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/trunk/SOURCE/pmc_interface_mod.f90	merged	eligible
/palm/branches/fricke/SOURCE/pmc_interface_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90	296-409
/palm/branches/suehring/pmc_interface_mod.f90	423-666

File size: 175.5 KB

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

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

Context Navigation

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

Download in other formats: