Home

Context Navigation

source: palm/trunk/SOURCE/vertical_nesting_mod.f90 @ 4444

Last change on this file since 4444 was 4444, checked in by raasch, 5 years ago
bugfix: cpp-directives for serial mode added
Property svn:keywords set to `Id`
File size: 149.8 KB

Line
1	!> @file vertical_nesting_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2020 Leibniz Universitaet Hannover
18	! Copyright 2017-2018 Karlsruhe Institute of Technology
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: vertical_nesting_mod.f90 4444 2020-03-05 15:59:50Z raasch $
28	! bugfix: cpp-directives for serial mode added
29	!
30	! 4360 2020-01-07 11:25:50Z suehring
31	! Corrected "Former revisions" section
32	!
33	! 4102 2019-07-17 16:00:03Z suehring
34	! - Slightly revise setting of boundary conditions at horizontal walls, use
35	! data-structure offset index instead of pre-calculate it for each facing
36	!
37	! 4101 2019-07-17 15:14:26Z gronemeier
38	! remove old_dt
39	!
40	! 3802 2019-03-17 13:33:42Z raasch
41	! unused subroutines commented out
42	!
43	! 3655 2019-01-07 16:51:22Z knoop
44	! unused variables removed
45	!
46	! 2365 2017-08-21 14:59:59Z kanani
47	! Initial revision (SadiqHuq)
48	!
49	!
50	! Description:
51	! ------------
52	!> Module for vertical nesting.
53	!>
54	!> Definition of parameters and variables for vertical nesting
55	!> The horizontal extent of the parent (Coarse Grid) and the child (Fine Grid)
56	!> have to be identical. The vertical extent of the FG should be smaller than CG.
57	!> Only integer nesting ratio supported. Odd nesting ratio preferred
58	!> The code follows MPI-1 standards. The available processors are split into
59	!> two groups using MPI_COMM_SPLIT. Exchange of data from CG to FG is called
60	!> interpolation. FG initialization by interpolation is done once at the start.
61	!> FG boundary conditions are set by interpolated at every timestep.
62	!> Exchange of data from CG to FG is called anterpolation, the two-way interaction
63	!> occurs at every timestep.
64	!> vnest_start_time set in PARIN allows delayed start of the coupling
65	!> after spin-up of the CG
66	!>
67	!> @todo Replace dz(1) appropriatly to account for grid stretching
68	!> @todo Ensure that code can be compiled for serial and parallel mode. Please
69	!> check the placement of the directive "__parallel".
70	!> @todo Add descriptions for all declared variables/parameters, one declaration
71	!> statement per variable
72	!> @todo Add a descriptive header above each subroutine (see land_surface_model)
73	!> @todo FORTRAN language statements must not be used as names for variables
74	!> (e.g. if). Please rename it.
75	!> @todo Revise code according to PALM Coding Standard
76	!------------------------------------------------------------------------------!
77	MODULE vertical_nesting_mod
78
79	USE kinds
80
81	IMPLICIT NONE
82
83	LOGICAL :: vnested = .FALSE. !> set to true if palmrun
84	!> provides specific information via stdin
85	LOGICAL :: vnest_init = .FALSE. !> set to true when FG is initialized
86	REAL(wp) :: vnest_start_time = 9999999.9 !> simulated time when FG should be
87	!> initialized. Should be
88	!> identical in PARIN & PARIN_N
89
90	#if defined( __parallel )
91
92	INTEGER(iwp),DIMENSION(3,2) :: bdims = 0 !> sub-domain grid topology of current PE
93	INTEGER(iwp),DIMENSION(3,2) :: bdims_rem = 0 !> sub-domain grid topology of partner PE
94	INTEGER(iwp) :: cg_nprocs !> no. of PE in CG. Set by palmrun -Y
95	INTEGER(iwp) :: fg_nprocs !> no. of PE in FG. Set by palmrun -Y
96	INTEGER(iwp) :: TYPE_VNEST_BC !> derived contiguous data type for interpolation
97	INTEGER(iwp) :: TYPE_VNEST_ANTER !> derived contiguous data type for anterpolation
98	INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE :: c2f_dims_cg !> One CG PE sends data to multiple FG PEs
99	!> list of grid-topology of partners
100	INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE :: f2c_dims_cg !> One CG PE receives data from multiple FG PEs
101	!> list of grid-topology of partners
102	INTEGER(iwp),DIMENSION(:),ALLOCATABLE :: c2f_dims_fg !> One FG PE sends data to multiple CG PE
103	!> list of grid-topology of partner
104	INTEGER(iwp),DIMENSION(:),ALLOCATABLE :: f2c_dims_fg !> One FG PE sends data to only one CG PE
105	!> list of grid-topology of partner
106
107	INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE :: f_rnk_lst !> list storing rank of FG PE denoted by pdims
108	INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE :: c_rnk_lst !> list storing rank of CG PE denoted by pdims
109	INTEGER(iwp),DIMENSION(3) :: cfratio !> Nesting ratio in x,y and z-directions
110
111	INTEGER(iwp) :: nxc !> no. of CG grid points in x-direction
112	INTEGER(iwp) :: nxf !> no. of FG grid points in x-direction
113	INTEGER(iwp) :: nyc !> no. of CG grid points in y-direction
114	INTEGER(iwp) :: nyf !> no. of FG grid points in y-direction
115	INTEGER(iwp) :: nzc !> no. of CG grid points in z-direction
116	INTEGER(iwp) :: nzf !> no. of FG grid points in z-direction
117	INTEGER(iwp) :: ngp_c !> no. of CG grid points in one vertical level
118	INTEGER(iwp) :: ngp_f !> no. of FG grid points in one vertical level
119
120	INTEGER(iwp) :: n_cell_c !> total no. of CG grid points in a PE
121	INTEGER(iwp),DIMENSION(2) :: pdims_partner !> processor topology of partner PE
122	INTEGER(iwp) :: target_idex !> temporary variable
123	INTEGER(iwp),DIMENSION(2) :: offset !> temporary variable
124	INTEGER(iwp),DIMENSION(2) :: map_coord !> temporary variable
125
126	REAL(wp) :: dxc !> CG grid pacing in x-direction
127	REAL(wp) :: dyc !> FG grid pacing in x-direction
128	REAL(wp) :: dxf !> CG grid pacing in y-direction
129	REAL(wp) :: dyf !> FG grid pacing in y-direction
130	REAL(wp) :: dzc !> CG grid pacing in z-direction
131	REAL(wp) :: dzf !> FG grid pacing in z-direction
132	REAL(wp) :: dtc !> dt calculated for CG
133	REAL(wp) :: dtf !> dt calculated for FG
134
135	REAL(wp), DIMENSION(:), ALLOCATABLE :: zuc !> CG vertical u-levels
136	REAL(wp), DIMENSION(:), ALLOCATABLE :: zuf !> FG vertical u-levels
137	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwc !> CG vertical w-levels
138	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwf !> FG vertical w-levels
139	REAL(wp), DIMENSION(:,:,:), POINTER :: interpol3d !> pointers to simplify function calls
140	REAL(wp), DIMENSION(:,:,:), POINTER :: anterpol3d !> pointers to simplify function calls
141
142
143	REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: work3d !> temporary array for exchange of 3D data
144	REAL(wp),DIMENSION(:,:), ALLOCATABLE :: work2dusws !> temporary array for exchange of 2D data
145	REAL(wp),DIMENSION(:,:), ALLOCATABLE :: work2dvsws !> temporary array for exchange of 2D data
146	REAL(wp),DIMENSION(:,:), ALLOCATABLE :: work2dts !> temporary array for exchange of 2D data
147	REAL(wp),DIMENSION(:,:), ALLOCATABLE :: work2dus !> temporary array for exchange of 2D data
148
149	SAVE
150
151	!-- Public functions
152	PUBLIC vnest_init_fine, vnest_boundary_conds, vnest_anterpolate, &
153	vnest_boundary_conds_khkm, vnest_anterpolate_e, &
154	vnest_init_pegrid_rank, vnest_init_pegrid_domain, vnest_init_grid, &
155	vnest_timestep_sync, vnest_deallocate
156
157	!-- Public constants and variables
158	PUBLIC vnested, vnest_init, vnest_start_time
159
160	PRIVATE bdims, bdims_rem, &
161	work3d, work2dusws, work2dvsws, work2dts, work2dus, &
162	dxc, dyc, dxf, dyf, dzc, dzf, dtc, dtf, &
163	zuc, zuf, zwc, zwf, interpol3d, anterpol3d, &
164	cg_nprocs, fg_nprocs, &
165	c2f_dims_cg, c2f_dims_fg, f2c_dims_cg, f2c_dims_fg, &
166	f_rnk_lst, c_rnk_lst, cfratio, pdims_partner, &
167	nxc, nxf, nyc, nyf, nzc, nzf, &
168	ngp_c, ngp_f, target_idex, n_cell_c, &
169	offset, map_coord, TYPE_VNEST_BC, TYPE_VNEST_ANTER
170
171	INTERFACE vnest_anterpolate
172	MODULE PROCEDURE vnest_anterpolate
173	END INTERFACE vnest_anterpolate
174
175	INTERFACE vnest_anterpolate_e
176	MODULE PROCEDURE vnest_anterpolate_e
177	END INTERFACE vnest_anterpolate_e
178
179	INTERFACE vnest_boundary_conds
180	MODULE PROCEDURE vnest_boundary_conds
181	END INTERFACE vnest_boundary_conds
182
183	INTERFACE vnest_boundary_conds_khkm
184	MODULE PROCEDURE vnest_boundary_conds_khkm
185	END INTERFACE vnest_boundary_conds_khkm
186
187	INTERFACE vnest_check_parameters
188	MODULE PROCEDURE vnest_check_parameters
189	END INTERFACE vnest_check_parameters
190
191	INTERFACE vnest_deallocate
192	MODULE PROCEDURE vnest_deallocate
193	END INTERFACE vnest_deallocate
194
195	INTERFACE vnest_init_fine
196	MODULE PROCEDURE vnest_init_fine
197	END INTERFACE vnest_init_fine
198
199	INTERFACE vnest_init_grid
200	MODULE PROCEDURE vnest_init_grid
201	END INTERFACE vnest_init_grid
202
203	INTERFACE vnest_init_pegrid_domain
204	MODULE PROCEDURE vnest_init_pegrid_domain
205	END INTERFACE vnest_init_pegrid_domain
206
207	INTERFACE vnest_init_pegrid_rank
208	MODULE PROCEDURE vnest_init_pegrid_rank
209	END INTERFACE vnest_init_pegrid_rank
210
211	INTERFACE vnest_timestep_sync
212	MODULE PROCEDURE vnest_timestep_sync
213	END INTERFACE vnest_timestep_sync
214
215	CONTAINS
216
217
218
219	SUBROUTINE vnest_init_fine
220	#if defined( __parallel )
221
222	!--------------------------------------------------------------------------------!
223	! Description:
224	! ------------
225	! At the specified vnest_start_time initialize the Fine Grid based on the coarse
226	! grid values
227	!------------------------------------------------------------------------------!
228
229
230	USE arrays_3d
231	USE control_parameters
232	USE grid_variables
233	USE indices
234	USE interfaces
235	USE pegrid
236	USE turbulence_closure_mod, &
237	ONLY : tcm_diffusivities
238
239
240	IMPLICIT NONE
241
242	REAL(wp) :: time_since_reference_point_rem
243	INTEGER(iwp) :: i
244	INTEGER(iwp) :: j
245	INTEGER(iwp) :: iif
246	INTEGER(iwp) :: jjf
247	INTEGER(iwp) :: kkf
248
249
250	if (myid ==0 ) print *, ' TIME TO INIT FINE from COARSE', simulated_time
251
252	!
253	!-- In case of model termination initiated by the remote model
254	!-- (terminate_coupled_remote > 0), initiate termination of the local model.
255	!-- The rest of the coupler must then be skipped because it would cause an MPI
256	!-- intercomminucation hang.
257	!-- If necessary, the coupler will be called at the beginning of the next
258	!-- restart run.
259
260	IF ( myid == 0) THEN
261	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &
262	target_id, 0, &
263	terminate_coupled_remote, 1, MPI_INTEGER, &
264	target_id, 0, &
265	comm_inter, status, ierr )
266	ENDIF
267	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
268	ierr )
269
270	IF ( terminate_coupled_remote > 0 ) THEN
271	WRITE( message_string, * ) 'remote model "', &
272	TRIM( coupling_mode_remote ), &
273	'" terminated', &
274	'&with terminate_coupled_remote = ', &
275	terminate_coupled_remote, &
276	'&local model "', TRIM( coupling_mode ), &
277	'" has', &
278	'&terminate_coupled = ', &
279	terminate_coupled
280	CALL message( 'vnest_init_fine', 'PA0310', 1, 2, 0, 6, 0 )
281	RETURN
282	ENDIF
283
284
285	!
286	!-- Exchange the current simulated time between the models,
287	!-- currently just for total_2ding
288	IF ( myid == 0 ) THEN
289
290	CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
291	11, comm_inter, ierr )
292	CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
293	target_id, 11, comm_inter, status, ierr )
294
295	ENDIF
296
297	CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
298	ierr )
299
300
301	IF ( coupling_mode == 'vnested_crse' ) THEN
302	!-- Send data to fine grid for initialization
303
304	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
305	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
306
307	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
308	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
309	map_coord(1) = i+offset(1)
310	map_coord(2) = j+offset(2)
311
312	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
313
314	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 10, &
315	comm_inter,status, ierr )
316
317	bdims (1,1) = bdims_rem (1,1) / cfratio(1)
318	bdims (1,2) = bdims_rem (1,2) / cfratio(1)
319	bdims (2,1) = bdims_rem (2,1) / cfratio(2)
320	bdims (2,2) = bdims_rem (2,2) / cfratio(2)
321	bdims (3,1) = bdims_rem (3,1)
322	bdims (3,2) = bdims_rem (3,2) / cfratio(3)
323
324
325	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 9, &
326	comm_inter, ierr )
327
328
329	n_cell_c = (bdims(1,2)-bdims(1,1)+3) * &
330	(bdims(2,2)-bdims(2,1)+3) * &
331	(bdims(3,2)-bdims(3,1)+3)
332
333	CALL MPI_SEND( u( bdims(3,1):bdims(3,2)+2, &
334	bdims(2,1)-1:bdims(2,2)+1, &
335	bdims(1,1)-1:bdims(1,2)+1),&
336	n_cell_c, MPI_REAL, target_idex, &
337	101, comm_inter, ierr)
338
339	CALL MPI_SEND( v( bdims(3,1):bdims(3,2)+2, &
340	bdims(2,1)-1:bdims(2,2)+1, &
341	bdims(1,1)-1:bdims(1,2)+1),&
342	n_cell_c, MPI_REAL, target_idex, &
343	102, comm_inter, ierr)
344
345	CALL MPI_SEND( w( bdims(3,1):bdims(3,2)+2, &
346	bdims(2,1)-1:bdims(2,2)+1, &
347	bdims(1,1)-1:bdims(1,2)+1),&
348	n_cell_c, MPI_REAL, target_idex, &
349	103, comm_inter, ierr)
350
351	CALL MPI_SEND( pt(bdims(3,1):bdims(3,2)+2, &
352	bdims(2,1)-1:bdims(2,2)+1, &
353	bdims(1,1)-1:bdims(1,2)+1),&
354	n_cell_c, MPI_REAL, target_idex, &
355	105, comm_inter, ierr)
356
357	IF ( humidity ) THEN
358	CALL MPI_SEND( q(bdims(3,1):bdims(3,2)+2, &
359	bdims(2,1)-1:bdims(2,2)+1, &
360	bdims(1,1)-1:bdims(1,2)+1),&
361	n_cell_c, MPI_REAL, target_idex, &
362	116, comm_inter, ierr)
363	ENDIF
364
365	CALL MPI_SEND( e( bdims(3,1):bdims(3,2)+2, &
366	bdims(2,1)-1:bdims(2,2)+1, &
367	bdims(1,1)-1:bdims(1,2)+1),&
368	n_cell_c, MPI_REAL, target_idex, &
369	104, comm_inter, ierr)
370
371	CALL MPI_SEND(kh( bdims(3,1):bdims(3,2)+2, &
372	bdims(2,1)-1:bdims(2,2)+1, &
373	bdims(1,1)-1:bdims(1,2)+1),&
374	n_cell_c, MPI_REAL, target_idex, &
375	106, comm_inter, ierr)
376
377	CALL MPI_SEND(km( bdims(3,1):bdims(3,2)+2, &
378	bdims(2,1)-1:bdims(2,2)+1, &
379	bdims(1,1)-1:bdims(1,2)+1),&
380	n_cell_c, MPI_REAL, target_idex, &
381	107, comm_inter, ierr)
382
383	!-- Send Surface fluxes
384	IF ( use_surface_fluxes ) THEN
385
386	n_cell_c = (bdims(1,2)-bdims(1,1)+3) * &
387	(bdims(2,2)-bdims(2,1)+3)
388
389	!
390	!-- shf and z0 for CG / FG need to initialized in input file or user_code
391	!-- TODO
392	!-- initialization of usws, vsws, ts and us not vital to vnest FG
393	!-- variables are not compatible with the new surface layer module
394	!
395	! CALL MPI_SEND(surf_def_h(0)%usws( bdims(2,1)-1:bdims(2,2)+1, &
396	! bdims(1,1)-1:bdims(1,2)+1),&
397	! n_cell_c, MPI_REAL, target_idex, &
398	! 110, comm_inter, ierr )
399	!
400	! CALL MPI_SEND(surf_def_h(0)%vsws( bdims(2,1)-1:bdims(2,2)+1, &
401	! bdims(1,1)-1:bdims(1,2)+1),&
402	! n_cell_c, MPI_REAL, target_idex, &
403	! 111, comm_inter, ierr )
404	!
405	! CALL MPI_SEND(ts ( bdims(2,1)-1:bdims(2,2)+1, &
406	! bdims(1,1)-1:bdims(1,2)+1),&
407	! n_cell_c, MPI_REAL, target_idex, &
408	! 112, comm_inter, ierr )
409	!
410	! CALL MPI_SEND(us ( bdims(2,1)-1:bdims(2,2)+1, &
411	! bdims(1,1)-1:bdims(1,2)+1),&
412	! n_cell_c, MPI_REAL, target_idex, &
413	! 113, comm_inter, ierr )
414	!
415	ENDIF
416
417
418
419
420	end do
421	end do
422
423	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
424	!-- Receive data from coarse grid for initialization
425
426	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
427	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
428	map_coord(1) = offset(1)
429	map_coord(2) = offset(2)
430	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
431
432	bdims (1,1) = nxl
433	bdims (1,2) = nxr
434	bdims (2,1) = nys
435	bdims (2,2) = nyn
436	bdims (3,1) = nzb
437	bdims (3,2) = nzt
438
439	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 10, &
440	comm_inter, ierr )
441
442	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 9, &
443	comm_inter,status, ierr )
444
445	n_cell_c = (bdims_rem(1,2)-bdims_rem(1,1)+3) * &
446	(bdims_rem(2,2)-bdims_rem(2,1)+3) * &
447	(bdims_rem(3,2)-bdims_rem(3,1)+3)
448
449	ALLOCATE( work3d ( bdims_rem(3,1) :bdims_rem(3,2)+2, &
450	bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
451	bdims_rem(1,1)-1:bdims_rem(1,2)+1))
452
453
454	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 101, &
455	comm_inter,status, ierr )
456	interpol3d => u
457	call interpolate_to_fine_u
458
459	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 102, &
460	comm_inter,status, ierr )
461	interpol3d => v
462	call interpolate_to_fine_v
463
464	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 103, &
465	comm_inter,status, ierr )
466	interpol3d => w
467	call interpolate_to_fine_w
468
469	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 105, &
470	comm_inter,status, ierr )
471	interpol3d => pt
472	call interpolate_to_fine_s
473
474	IF ( humidity ) THEN
475	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 116, &
476	comm_inter,status, ierr )
477	interpol3d => q
478	call interpolate_to_fine_s
479	ENDIF
480
481	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 104, &
482	comm_inter,status, ierr )
483	interpol3d => e
484	call interpolate_to_fine_s
485
486	!-- kh,km no target attribute, use of pointer not possible
487	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 106, &
488	comm_inter,status, ierr )
489	call interpolate_to_fine_kh
490
491	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 107, &
492	comm_inter,status, ierr )
493	call interpolate_to_fine_km
494
495	DEALLOCATE( work3d )
496	NULLIFY ( interpol3d )
497
498	!-- Recv Surface Fluxes
499	IF ( use_surface_fluxes ) THEN
500	n_cell_c = (bdims_rem(1,2)-bdims_rem(1,1)+3) * &
501	(bdims_rem(2,2)-bdims_rem(2,1)+3)
502
503	ALLOCATE( work2dusws ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
504	bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
505	ALLOCATE( work2dvsws ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
506	bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
507	ALLOCATE( work2dts ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
508	bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
509	ALLOCATE( work2dus ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
510	bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
511
512	!
513	!-- shf and z0 for CG / FG need to initialized in input file or user_code
514	!-- TODO
515	!-- initialization of usws, vsws, ts and us not vital to vnest FG
516	!-- variables are not compatible with the new surface layer module
517	!
518	! CALL MPI_RECV( work2dusws,n_cell_c, MPI_REAL, target_idex, 110, &
519	! comm_inter,status, ierr )
520	!
521	! CALL MPI_RECV( work2dvsws,n_cell_c, MPI_REAL, target_idex, 111, &
522	! comm_inter,status, ierr )
523	!
524	! CALL MPI_RECV( work2dts ,n_cell_c, MPI_REAL, target_idex, 112, &
525	! comm_inter,status, ierr )
526	!
527	! CALL MPI_RECV( work2dus ,n_cell_c, MPI_REAL, target_idex, 113, &
528	! comm_inter,status, ierr )
529	!
530	! CALL interpolate_to_fine_flux ( 108 )
531
532	DEALLOCATE( work2dusws )
533	DEALLOCATE( work2dvsws )
534	DEALLOCATE( work2dts )
535	DEALLOCATE( work2dus )
536	ENDIF
537
538	IF ( .NOT. constant_diffusion ) THEN
539	DO kkf = bdims(3,1)+1,bdims(3,2)+1
540	DO jjf = bdims(2,1),bdims(2,2)
541	DO iif = bdims(1,1),bdims(1,2)
542
543	IF ( e(kkf,jjf,iif) < 0.0 ) THEN
544	e(kkf,jjf,iif) = 1E-15_wp
545	END IF
546
547	END DO
548	END DO
549	END DO
550	ENDIF
551
552	w(nzt+1,:,:) = w(nzt,:,:)
553
554	CALL exchange_horiz( u, nbgp )
555	CALL exchange_horiz( v, nbgp )
556	CALL exchange_horiz( w, nbgp )
557	CALL exchange_horiz( pt, nbgp )
558	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e, nbgp )
559	IF ( humidity ) CALL exchange_horiz( q, nbgp )
560
561	!
562	!-- Velocity boundary conditions at the bottom boundary
563	IF ( ibc_uv_b == 0 ) THEN
564	u(nzb,:,:) = 0.0_wp
565	v(nzb,:,:) = 0.0_wp
566	ELSE
567	u(nzb,:,:) = u(nzb+1,:,:)
568	v(nzb,:,:) = v(nzb+1,:,:)
569	END IF
570
571
572	w(nzb,:,:) = 0.0_wp
573
574	!
575	!-- Temperature boundary conditions at the bottom boundary
576	IF ( ibc_pt_b /= 0 ) THEN
577	pt(nzb,:,:) = pt(nzb+1,:,:)
578	END IF
579
580	!
581	!-- Bottom boundary condition for the turbulent kinetic energy
582	!-- Generally a Neumann condition with de/dz=0 is assumed
583	IF ( .NOT. constant_diffusion ) THEN
584	e(nzb,:,:) = e(nzb+1,:,:)
585	END IF
586
587	!
588	!-- Bottom boundary condition for turbulent diffusion coefficients
589	km(nzb,:,:) = km(nzb+1,:,:)
590	kh(nzb,:,:) = kh(nzb+1,:,:)
591
592	!diffusivities required
593	IF ( .NOT. humidity ) THEN
594	CALL tcm_diffusivities( pt, pt_reference )
595	ELSE
596	CALL tcm_diffusivities( vpt, pt_reference )
597	ENDIF
598
599
600	!
601	!-- Reset Fine Grid top Boundary Condition
602	!-- At the top of the FG, the scalars always follow Dirichlet condition
603
604	ibc_pt_t = 0
605
606	!-- Initialize old time levels
607	pt_p = pt; u_p = u; v_p = v; w_p = w
608	IF ( .NOT. constant_diffusion ) e_p = e
609	IF ( humidity ) THEN
610	ibc_q_t = 0
611	q_p = q
612	ENDIF
613
614	ENDIF
615
616
617	if (myid==0) print , '* Fine Initalized ** simulated_time:', simulated_time
618
619	CONTAINS
620
621	SUBROUTINE interpolate_to_fine_w
622
623	USE arrays_3d
624	USE control_parameters
625	USE grid_variables
626	USE indices
627	USE pegrid
628
629
630	IMPLICIT NONE
631
632	INTEGER(iwp) :: i
633	INTEGER(iwp) :: j
634	INTEGER(iwp) :: k
635	INTEGER(iwp) :: iif
636	INTEGER(iwp) :: jjf
637	INTEGER(iwp) :: kkf
638	INTEGER(iwp) :: nzbottom
639	INTEGER(iwp) :: nztop
640	INTEGER(iwp) :: bottomx
641	INTEGER(iwp) :: bottomy
642	INTEGER(iwp) :: bottomz
643	INTEGER(iwp) :: topx
644	INTEGER(iwp) :: topy
645	INTEGER(iwp) :: topz
646	REAL(wp) :: eps
647	REAL(wp) :: alpha
648	REAL(wp) :: eminus
649	REAL(wp) :: edot
650	REAL(wp) :: eplus
651	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprs
652	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: wprf
653
654
655	nzbottom = bdims_rem (3,1)
656	nztop = bdims_rem (3,2)
657
658	ALLOCATE( wprf(nzbottom:nztop, bdims_rem(2,1)-1: bdims_rem(2,2)+1,nxl:nxr) )
659	ALLOCATE( wprs(nzbottom:nztop,nys:nyn,nxl:nxr) )
660
661
662	!
663	!-- Initialisation of the velocity component w
664	!
665	!-- Interpolation in x-direction
666	DO k = nzbottom, nztop
667	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
668	DO i = bdims_rem(1,1),bdims_rem(1,2)
669
670	bottomx = (nxf+1)/(nxc+1) * i
671	topx = (nxf+1)/(nxc+1) * (i+1) - 1
672
673	DO iif = bottomx, topx
674
675	eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
676	alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
677	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
678	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
679	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
680
681	wprf(k,j,iif) = eminus * work3d(k,j,i-1) &
682	+ edot * work3d(k,j,i) &
683	+ eplus * work3d(k,j,i+1)
684	END DO
685
686	END DO
687	END DO
688	END DO
689
690	!
691	!-- Interpolation in y-direction (quadratic, Clark and Farley)
692	DO k = nzbottom, nztop
693	DO j = bdims_rem(2,1), bdims_rem(2,2)
694
695	bottomy = (nyf+1)/(nyc+1) * j
696	topy = (nyf+1)/(nyc+1) * (j+1) - 1
697
698	DO iif = nxl, nxr
699	DO jjf = bottomy, topy
700
701	eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
702	alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
703	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
704	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
705	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
706
707	wprs(k,jjf,iif) = eminus * wprf(k,j-1,iif) &
708	+ edot * wprf(k,j,iif) &
709	+ eplus * wprf(k,j+1,iif)
710
711	END DO
712	END DO
713
714	END DO
715	END DO
716
717	!
718	!-- Interpolation in z-direction (linear)
719
720	DO k = nzbottom, nztop-1
721
722	bottomz = (dzc/dzf) * k
723	topz = (dzc/dzf) * (k+1) - 1
724
725	DO jjf = nys, nyn
726	DO iif = nxl, nxr
727	DO kkf = bottomz, topz
728
729	w(kkf,jjf,iif) = wprs(k,jjf,iif) + ( zwf(kkf) - zwc(k) ) &
730	* ( wprs(k+1,jjf,iif) - wprs(k,jjf,iif) ) / dzc
731
732	END DO
733	END DO
734	END DO
735
736	END DO
737
738	DO jjf = nys, nyn
739	DO iif = nxl, nxr
740
741	w(nzt,jjf,iif) = wprs(nztop,jjf,iif)
742
743	END DO
744	END DO
745	!
746	! w(nzb:nzt+1,nys:nyn,nxl:nxr) = 0
747
748	DEALLOCATE( wprf, wprs )
749
750	END SUBROUTINE interpolate_to_fine_w
751
752	SUBROUTINE interpolate_to_fine_u
753
754
755	USE arrays_3d
756	USE control_parameters
757	USE grid_variables
758	USE indices
759	USE pegrid
760
761
762	IMPLICIT NONE
763
764	INTEGER(iwp) :: i
765	INTEGER(iwp) :: j
766	INTEGER(iwp) :: k
767	INTEGER(iwp) :: iif
768	INTEGER(iwp) :: jjf
769	INTEGER(iwp) :: kkf
770	INTEGER(iwp) :: nzbottom
771	INTEGER(iwp) :: nztop
772	INTEGER(iwp) :: bottomx
773	INTEGER(iwp) :: bottomy
774	INTEGER(iwp) :: bottomz
775	INTEGER(iwp) :: topx
776	INTEGER(iwp) :: topy
777	INTEGER(iwp) :: topz
778	REAL(wp) :: eps
779	REAL(wp) :: alpha
780	REAL(wp) :: eminus
781	REAL(wp) :: edot
782	REAL(wp) :: eplus
783	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprf
784	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprs
785
786
787
788	nzbottom = bdims_rem (3,1)
789	nztop = bdims_rem (3,2)
790
791	ALLOCATE( uprf(nzbottom:nztop+2,nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
792	ALLOCATE( uprs(nzb+1:nzt+1,nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
793
794	!
795	!-- Initialisation of the velocity component uf
796
797	!
798	!-- Interpolation in y-direction (quadratic, Clark and Farley)
799
800	DO k = nzbottom, nztop+2
801	DO j = bdims_rem(2,1), bdims_rem(2,2)
802
803	bottomy = (nyf+1)/(nyc+1) * j
804	topy = (nyf+1)/(nyc+1) * (j+1) - 1
805
806	DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
807	DO jjf = bottomy, topy
808
809	eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
810	alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
811	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
812	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
813	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
814
815	uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
816	+ edot * work3d(k,j,i) &
817	+ eplus * work3d(k,j+1,i)
818
819	END DO
820	END DO
821
822	END DO
823	END DO
824
825	!
826	!-- Interpolation in z-direction (quadratic, Clark and Farley)
827
828	DO k = nzbottom+1, nztop
829
830	bottomz = (dzc/dzf) * (k-1) + 1
831	topz = (dzc/dzf) * k
832
833	DO jjf = nys, nyn
834	DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
835	DO kkf = bottomz, topz
836
837	eps = ( zuf(kkf) - zuc(k) ) / dzc
838	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
839	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
840	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
841	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
842
843	uprs(kkf,jjf,i) = eminus * uprf(k-1,jjf,i) &
844	+ edot * uprf(k,jjf,i) &
845	+ eplus * uprf(k+1,jjf,i)
846
847	END DO
848	END DO
849	END DO
850
851	END DO
852
853	DO jjf = nys, nyn
854	DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
855
856	eps = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
857	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
858	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
859	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
860	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
861
862	uprs(nzt+1,jjf,i) = eminus * uprf(nztop,jjf,i) &
863	+ edot * uprf(nztop+1,jjf,i) &
864	+ eplus * uprf(nztop+2,jjf,i)
865
866	END DO
867	END DO
868
869	!
870	!-- Interpolation in x-direction (linear)
871
872	DO kkf = nzb+1, nzt+1
873	DO jjf = nys, nyn
874	DO i = bdims_rem(1,1), bdims_rem(1,2)
875
876	bottomx = (nxf+1)/(nxc+1) * i
877	topx = (nxf+1)/(nxc+1) * (i+1) - 1
878
879	DO iif = bottomx, topx
880	u(kkf,jjf,iif) = uprs(kkf,jjf,i) + ( iif * dxf - i * dxc ) &
881	* ( uprs(kkf,jjf,i+1) - uprs(kkf,jjf,i) ) / dxc
882	END DO
883
884	END DO
885	END DO
886	END DO
887	!
888	!-- Determination of uf at the bottom boundary
889
890
891
892	DEALLOCATE( uprf, uprs )
893
894	END SUBROUTINE interpolate_to_fine_u
895
896
897	SUBROUTINE interpolate_to_fine_v
898
899
900	USE arrays_3d
901	USE control_parameters
902	USE grid_variables
903	USE indices
904	USE pegrid
905
906
907	IMPLICIT NONE
908
909	INTEGER(iwp) :: i
910	INTEGER(iwp) :: j
911	INTEGER(iwp) :: k
912	INTEGER(iwp) :: iif
913	INTEGER(iwp) :: jjf
914	INTEGER(iwp) :: kkf
915	INTEGER(iwp) :: nzbottom
916	INTEGER(iwp) :: nztop
917	INTEGER(iwp) :: bottomx
918	INTEGER(iwp) :: bottomy
919	INTEGER(iwp) :: bottomz
920	INTEGER(iwp) :: topx
921	INTEGER(iwp) :: topy
922	INTEGER(iwp) :: topz
923	REAL(wp) :: eps
924	REAL(wp) :: alpha
925	REAL(wp) :: eminus
926	REAL(wp) :: edot
927	REAL(wp) :: eplus
928	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprs
929	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprf
930
931
932	nzbottom = bdims_rem (3,1)
933	nztop = bdims_rem (3,2)
934
935	ALLOCATE( vprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
936	ALLOCATE( vprs(nzb+1:nzt+1, bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
937	!
938	!-- Initialisation of the velocity component vf
939
940	!
941	!-- Interpolation in x-direction (quadratic, Clark and Farley)
942
943	DO k = nzbottom, nztop+2
944	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
945	DO i = bdims_rem(1,1), bdims_rem(1,2)
946
947	bottomx = (nxf+1)/(nxc+1) * i
948	topx = (nxf+1)/(nxc+1) * (i+1) - 1
949
950	DO iif = bottomx, topx
951
952	eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
953	alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
954	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
955	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
956	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
957
958	vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
959	+ edot * work3d(k,j,i) &
960	+ eplus * work3d(k,j,i+1)
961
962	END DO
963
964	END DO
965	END DO
966	END DO
967
968	!
969	!-- Interpolation in z-direction (quadratic, Clark and Farley)
970
971	DO k = nzbottom+1, nztop
972
973	bottomz = (dzc/dzf) * (k-1) + 1
974	topz = (dzc/dzf) * k
975
976	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
977	DO iif = nxl, nxr
978	DO kkf = bottomz, topz
979
980	eps = ( zuf(kkf) - zuc(k) ) / dzc
981	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
982	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
983	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
984	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
985
986	vprs(kkf,j,iif) = eminus * vprf(k-1,j,iif) &
987	+ edot * vprf(k,j,iif) &
988	+ eplus * vprf(k+1,j,iif)
989
990	END DO
991	END DO
992	END DO
993
994	END DO
995
996	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
997	DO iif = nxl, nxr
998
999	eps = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
1000	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1001	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1002	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1003	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1004
1005	vprs(nzt+1,j,iif) = eminus * vprf(nztop,j,iif) &
1006	+ edot * vprf(nztop+1,j,iif) &
1007	+ eplus * vprf(nztop+2,j,iif)
1008
1009	END DO
1010	END DO
1011
1012	!
1013	!-- Interpolation in y-direction (linear)
1014
1015	DO kkf = nzb+1, nzt+1
1016	DO j = bdims_rem(2,1), bdims_rem(2,2)
1017
1018	bottomy = (nyf+1)/(nyc+1) * j
1019	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1020
1021	DO iif = nxl, nxr
1022	DO jjf = bottomy, topy
1023	v (kkf,jjf,iif) = vprs(kkf,j,iif) + ( jjf * dyf - j * dyc ) &
1024	* ( vprs(kkf,j+1,iif) - vprs(kkf,j,iif) ) / dyc
1025	END DO
1026	END DO
1027
1028	END DO
1029	END DO
1030
1031	!
1032	!-- Determination of vf at the bottom boundary
1033
1034
1035	DEALLOCATE( vprf, vprs )
1036
1037	END SUBROUTINE interpolate_to_fine_v
1038
1039
1040	SUBROUTINE interpolate_to_fine_s
1041
1042
1043	USE arrays_3d
1044	USE control_parameters
1045	USE grid_variables
1046	USE indices
1047	USE pegrid
1048
1049
1050	IMPLICIT NONE
1051
1052	INTEGER(iwp) :: i
1053	INTEGER(iwp) :: j
1054	INTEGER(iwp) :: k
1055	INTEGER(iwp) :: iif
1056	INTEGER(iwp) :: jjf
1057	INTEGER(iwp) :: kkf
1058	INTEGER(iwp) :: nzbottom
1059	INTEGER(iwp) :: nztop
1060	INTEGER(iwp) :: bottomx
1061	INTEGER(iwp) :: bottomy
1062	INTEGER(iwp) :: bottomz
1063	INTEGER(iwp) :: topx
1064	INTEGER(iwp) :: topy
1065	INTEGER(iwp) :: topz
1066	REAL(wp) :: eps
1067	REAL(wp) :: alpha
1068	REAL(wp) :: eminus
1069	REAL(wp) :: edot
1070	REAL(wp) :: eplus
1071	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
1072	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
1073
1074
1075	nzbottom = bdims_rem (3,1)
1076	nztop = bdims_rem (3,2)
1077
1078	ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1079	ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
1080
1081	!
1082	!-- Initialisation of scalar variables
1083
1084	!
1085	!-- Interpolation in x-direction (quadratic, Clark and Farley)
1086
1087	DO k = nzbottom, nztop+2
1088	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
1089	DO i = bdims_rem(1,1), bdims_rem(1,2)
1090
1091	bottomx = (nxf+1)/(nxc+1) * i
1092	topx = (nxf+1)/(nxc+1) * (i+1) - 1
1093
1094	DO iif = bottomx, topx
1095
1096	eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
1097	alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
1098	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1099	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1100	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1101
1102	ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
1103	+ edot * work3d(k,j,i) &
1104	+ eplus * work3d(k,j,i+1)
1105	END DO
1106
1107	END DO
1108	END DO
1109	END DO
1110
1111	!
1112	!-- Interpolation in y-direction (quadratic, Clark and Farley)
1113
1114	DO k = nzbottom, nztop+2
1115	DO j = bdims_rem(2,1), bdims_rem(2,2)
1116
1117	bottomy = (nyf+1)/(nyc+1) * j
1118	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1119
1120	DO iif = nxl, nxr
1121	DO jjf = bottomy, topy
1122
1123	eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
1124	alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
1125	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1126	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1127	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1128
1129	ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1130	+ edot * ptprf(k,j,iif) &
1131	+ eplus * ptprf(k,j+1,iif)
1132
1133	END DO
1134	END DO
1135
1136	END DO
1137	END DO
1138
1139	!
1140	!-- Interpolation in z-direction (quadratic, Clark and Farley)
1141
1142	DO k = nzbottom+1, nztop
1143
1144	bottomz = (dzc/dzf) * (k-1) + 1
1145	topz = (dzc/dzf) * k
1146
1147	DO jjf = nys, nyn
1148	DO iif = nxl, nxr
1149	DO kkf = bottomz, topz
1150
1151	eps = ( zuf(kkf) - zuc(k) ) / dzc
1152	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1153	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1154	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1155	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1156
1157	interpol3d(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1158	+ edot * ptprs(k,jjf,iif) &
1159	+ eplus * ptprs(k+1,jjf,iif)
1160
1161	END DO
1162	END DO
1163	END DO
1164
1165	END DO
1166
1167	DO jjf = nys, nyn
1168	DO iif = nxl, nxr
1169
1170	eps = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
1171	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1172	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1173	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1174	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1175
1176	interpol3d(nzt+1,jjf,iif) = eminus * ptprs(nztop,jjf,iif) &
1177	+ edot * ptprs(nztop+1,jjf,iif) &
1178	+ eplus * ptprs(nztop+2,jjf,iif)
1179
1180	END DO
1181	END DO
1182
1183
1184	DEALLOCATE( ptprf, ptprs )
1185
1186	END SUBROUTINE interpolate_to_fine_s
1187
1188
1189	SUBROUTINE interpolate_to_fine_kh
1190
1191
1192	USE arrays_3d
1193	USE control_parameters
1194	USE grid_variables
1195	USE indices
1196	USE pegrid
1197
1198
1199	IMPLICIT NONE
1200
1201	INTEGER(iwp) :: i
1202	INTEGER(iwp) :: j
1203	INTEGER(iwp) :: k
1204	INTEGER(iwp) :: iif
1205	INTEGER(iwp) :: jjf
1206	INTEGER(iwp) :: kkf
1207	INTEGER(iwp) :: nzbottom
1208	INTEGER(iwp) :: nztop
1209	INTEGER(iwp) :: bottomx
1210	INTEGER(iwp) :: bottomy
1211	INTEGER(iwp) :: bottomz
1212	INTEGER(iwp) :: topx
1213	INTEGER(iwp) :: topy
1214	INTEGER(iwp) :: topz
1215	REAL(wp) :: eps
1216	REAL(wp) :: alpha
1217	REAL(wp) :: eminus
1218	REAL(wp) :: edot
1219	REAL(wp) :: eplus
1220	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
1221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
1222
1223
1224	nzbottom = bdims_rem (3,1)
1225	nztop = bdims_rem (3,2)
1226	! nztop = blk_dim_rem (3,2)+1
1227
1228
1229	ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1230	ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
1231
1232
1233	!
1234	!-- Initialisation of scalar variables
1235
1236	!
1237	!-- Interpolation in x-direction (quadratic, Clark and Farley)
1238
1239	DO k = nzbottom, nztop+2
1240	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
1241	DO i = bdims_rem(1,1), bdims_rem(1,2)
1242
1243	bottomx = (nxf+1)/(nxc+1) * i
1244	topx = (nxf+1)/(nxc+1) * (i+1) - 1
1245
1246	DO iif = bottomx, topx
1247
1248	eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
1249	alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
1250	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1251	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1252	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1253
1254	ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
1255	+ edot * work3d(k,j,i) &
1256	+ eplus * work3d(k,j,i+1)
1257	END DO
1258
1259	END DO
1260	END DO
1261	END DO
1262
1263	!
1264	!-- Interpolation in y-direction (quadratic, Clark and Farley)
1265
1266	DO k = nzbottom, nztop+2
1267	DO j = bdims_rem(2,1), bdims_rem(2,2)
1268
1269	bottomy = (nyf+1)/(nyc+1) * j
1270	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1271
1272	DO iif = nxl, nxr
1273	DO jjf = bottomy, topy
1274
1275	eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
1276	alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
1277	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1278	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1279	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1280
1281	ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1282	+ edot * ptprf(k,j,iif) &
1283	+ eplus * ptprf(k,j+1,iif)
1284
1285	END DO
1286	END DO
1287
1288	END DO
1289	END DO
1290
1291	!
1292	!-- Interpolation in z-direction (quadratic, Clark and Farley)
1293
1294	DO k = nzbottom+1, nztop
1295
1296	bottomz = (dzc/dzf) * (k-1) + 1
1297	topz = (dzc/dzf) * k
1298
1299	DO jjf = nys, nyn
1300	DO iif = nxl, nxr
1301	DO kkf = bottomz, topz
1302
1303	eps = ( zuf(kkf) - zuc(k) ) / dzc
1304	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1305	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1306	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1307	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1308
1309	kh(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1310	+ edot * ptprs(k,jjf,iif) &
1311	+ eplus * ptprs(k+1,jjf,iif)
1312
1313	END DO
1314	END DO
1315	END DO
1316
1317	END DO
1318
1319	DO jjf = nys, nyn
1320	DO iif = nxl, nxr
1321
1322	eps = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
1323	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1324	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1325	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1326	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1327
1328	kh(nzt+1,jjf,iif) = eminus * ptprs(nztop,jjf,iif) &
1329	+ edot * ptprs(nztop+1,jjf,iif) &
1330	+ eplus * ptprs(nztop+2,jjf,iif)
1331
1332	END DO
1333	END DO
1334
1335
1336	DEALLOCATE( ptprf, ptprs )
1337
1338	END SUBROUTINE interpolate_to_fine_kh
1339
1340	SUBROUTINE interpolate_to_fine_km
1341
1342
1343	USE arrays_3d
1344	USE control_parameters
1345	USE grid_variables
1346	USE indices
1347	USE pegrid
1348
1349
1350	IMPLICIT NONE
1351
1352	INTEGER(iwp) :: i
1353	INTEGER(iwp) :: j
1354	INTEGER(iwp) :: k
1355	INTEGER(iwp) :: iif
1356	INTEGER(iwp) :: jjf
1357	INTEGER(iwp) :: kkf
1358	INTEGER(iwp) :: nzbottom
1359	INTEGER(iwp) :: nztop
1360	INTEGER(iwp) :: bottomx
1361	INTEGER(iwp) :: bottomy
1362	INTEGER(iwp) :: bottomz
1363	INTEGER(iwp) :: topx
1364	INTEGER(iwp) :: topy
1365	INTEGER(iwp) :: topz
1366	REAL(wp) :: eps
1367	REAL(wp) :: alpha
1368	REAL(wp) :: eminus
1369	REAL(wp) :: edot
1370	REAL(wp) :: eplus
1371	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
1372	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
1373
1374
1375	nzbottom = bdims_rem (3,1)
1376	nztop = bdims_rem (3,2)
1377	! nztop = blk_dim_rem (3,2)+1
1378
1379
1380	ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1381	ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
1382
1383
1384	!
1385	!-- Initialisation of scalar variables
1386
1387	!
1388	!-- Interpolation in x-direction (quadratic, Clark and Farley)
1389
1390	DO k = nzbottom, nztop+2
1391	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
1392	DO i = bdims_rem(1,1), bdims_rem(1,2)
1393
1394	bottomx = (nxf+1)/(nxc+1) * i
1395	topx = (nxf+1)/(nxc+1) * (i+1) - 1
1396
1397	DO iif = bottomx, topx
1398
1399	eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
1400	alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
1401	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1402	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1403	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1404
1405	ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
1406	+ edot * work3d(k,j,i) &
1407	+ eplus * work3d(k,j,i+1)
1408	END DO
1409
1410	END DO
1411	END DO
1412	END DO
1413
1414	!
1415	!-- Interpolation in y-direction (quadratic, Clark and Farley)
1416
1417	DO k = nzbottom, nztop+2
1418	DO j = bdims_rem(2,1), bdims_rem(2,2)
1419
1420	bottomy = (nyf+1)/(nyc+1) * j
1421	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1422
1423	DO iif = nxl, nxr
1424	DO jjf = bottomy, topy
1425
1426	eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
1427	alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
1428	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1429	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1430	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1431
1432	ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1433	+ edot * ptprf(k,j,iif) &
1434	+ eplus * ptprf(k,j+1,iif)
1435
1436	END DO
1437	END DO
1438
1439	END DO
1440	END DO
1441
1442	!
1443	!-- Interpolation in z-direction (quadratic, Clark and Farley)
1444
1445	DO k = nzbottom+1, nztop
1446
1447	bottomz = (dzc/dzf) * (k-1) + 1
1448	topz = (dzc/dzf) * k
1449
1450	DO jjf = nys, nyn
1451	DO iif = nxl, nxr
1452	DO kkf = bottomz, topz
1453
1454	eps = ( zuf(kkf) - zuc(k) ) / dzc
1455	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1456	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1457	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1458	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1459
1460	km(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1461	+ edot * ptprs(k,jjf,iif) &
1462	+ eplus * ptprs(k+1,jjf,iif)
1463
1464	END DO
1465	END DO
1466	END DO
1467
1468	END DO
1469
1470	DO jjf = nys, nyn
1471	DO iif = nxl, nxr
1472
1473	eps = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
1474	alpha = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
1475	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1476	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1477	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1478
1479	km(nzt+1,jjf,iif) = eminus * ptprs(nztop,jjf,iif) &
1480	+ edot * ptprs(nztop+1,jjf,iif) &
1481	+ eplus * ptprs(nztop+2,jjf,iif)
1482
1483	END DO
1484	END DO
1485
1486
1487	DEALLOCATE( ptprf, ptprs )
1488
1489	END SUBROUTINE interpolate_to_fine_km
1490
1491
1492
1493
1494	! SUBROUTINE interpolate_to_fine_flux
1495	!
1496	!
1497	! USE arrays_3d
1498	! USE control_parameters
1499	! USE grid_variables
1500	! USE indices
1501	! USE pegrid
1502	!
1503	!
1504	! IMPLICIT NONE
1505	!
1506	! INTEGER(iwp) :: i
1507	! INTEGER(iwp) :: j
1508	! INTEGER(iwp) :: iif
1509	! INTEGER(iwp) :: jjf
1510	! INTEGER(iwp) :: bottomx
1511	! INTEGER(iwp) :: bottomy
1512	! INTEGER(iwp) :: topx
1513	! INTEGER(iwp) :: topy
1514	! REAL(wp) :: eps
1515	! REAL(wp) :: alpha
1516	! REAL(wp) :: eminus
1517	! REAL(wp) :: edot
1518	! REAL(wp) :: eplus
1519	! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uswspr
1520	! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vswspr
1521	! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tspr
1522	! REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uspr
1523	!
1524	! ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1525	! ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1526	! ALLOCATE( tspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1527	! ALLOCATE( uspr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1528	!
1529	! !
1530	! !-- Initialisation of scalar variables (2D)
1531	!
1532	! !
1533	! !-- Interpolation in x-direction (quadratic, Clark and Farley)
1534	!
1535	! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
1536	! DO i = bdims_rem(1,1), bdims_rem(1,2)
1537	!
1538	! bottomx = (nxf+1)/(nxc+1) * i
1539	! topx = (nxf+1)/(nxc+1) * (i+1) - 1
1540	!
1541	! DO iif = bottomx, topx
1542	!
1543	! eps = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
1544	! alpha = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
1545	! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1546	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1547	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1548	!
1549	! uswspr(j,iif) = eminus * work2dusws(j,i-1) &
1550	! + edot * work2dusws(j,i) &
1551	! + eplus * work2dusws(j,i+1)
1552	!
1553	! vswspr(j,iif) = eminus * work2dvsws(j,i-1) &
1554	! + edot * work2dvsws(j,i) &
1555	! + eplus * work2dvsws(j,i+1)
1556	!
1557	! tspr(j,iif) = eminus * work2dts(j,i-1) &
1558	! + edot * work2dts(j,i) &
1559	! + eplus * work2dts(j,i+1)
1560	!
1561	! uspr(j,iif) = eminus * work2dus(j,i-1) &
1562	! + edot * work2dus(j,i) &
1563	! + eplus * work2dus(j,i+1)
1564	!
1565	! END DO
1566	!
1567	! END DO
1568	! END DO
1569	!
1570	! !
1571	! !-- Interpolation in y-direction (quadratic, Clark and Farley)
1572	!
1573	! DO j = bdims_rem(2,1), bdims_rem(2,2)
1574	!
1575	! bottomy = (nyf+1)/(nyc+1) * j
1576	! topy = (nyf+1)/(nyc+1) * (j+1) - 1
1577	!
1578	! DO iif = nxl, nxr
1579	! DO jjf = bottomy, topy
1580	!
1581	! eps = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
1582	! alpha = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
1583	! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1584	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1585	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1586	!
1587	!!
1588	!!-- TODO
1589	!-- variables are not compatible with the new surface layer module
1590	!
1591	! surf_def_h(0)%usws(jjf,iif) = eminus * uswspr(j-1,if) &
1592	! + edot * uswspr(j,iif) &
1593	! + eplus * uswspr(j+1,iif)
1594	!
1595	! surf_def_h(0)%vsws(jjf,iif) = eminus * vswspr(j-1,if) &
1596	! + edot * vswspr(j,iif) &
1597	! + eplus * vswspr(j+1,iif)
1598	!
1599	! ts(jjf,iif) = eminus * tspr(j-1,if) &
1600	! + edot * tspr(j,iif) &
1601	! + eplus * tspr(j+1,iif)
1602	!
1603	! us(jjf,iif) = eminus * uspr(j-1,if) &
1604	! + edot * uspr(j,iif) &
1605	! + eplus * uspr(j+1,iif)
1606	!
1607	! END DO
1608	! END DO
1609	!
1610	! END DO
1611	!
1612	!
1613	! DEALLOCATE( uswspr, vswspr )
1614	! DEALLOCATE( tspr, uspr )
1615	!
1616	!
1617	! END SUBROUTINE interpolate_to_fine_flux
1618
1619
1620	#endif
1621	END SUBROUTINE vnest_init_fine
1622
1623	SUBROUTINE vnest_boundary_conds
1624	#if defined( __parallel )
1625	!------------------------------------------------------------------------------!
1626	! Description:
1627	! ------------
1628	! Boundary conditions for the prognostic quantities.
1629	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
1630	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
1631	! handled in routine exchange_horiz. Pressure boundary conditions are
1632	! explicitly set in routines pres, poisfft, poismg and sor.
1633	!------------------------------------------------------------------------------!
1634
1635	USE arrays_3d
1636	USE control_parameters
1637	USE grid_variables
1638	USE indices
1639	USE pegrid
1640
1641
1642	IMPLICIT NONE
1643
1644	INTEGER(iwp) :: i
1645	INTEGER(iwp) :: j
1646	INTEGER(iwp) :: iif
1647	INTEGER(iwp) :: jjf
1648
1649
1650	!
1651	!-- vnest: top boundary conditions
1652
1653	IF ( coupling_mode == 'vnested_crse' ) THEN
1654	!-- Send data to fine grid for TOP BC
1655
1656	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
1657	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
1658
1659	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
1660	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
1661	map_coord(1) = i+offset(1)
1662	map_coord(2) = j+offset(2)
1663
1664	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
1665
1666	bdims (1,1) = c2f_dims_cg (0,map_coord(1),map_coord(2))
1667	bdims (1,2) = c2f_dims_cg (1,map_coord(1),map_coord(2))
1668	bdims (2,1) = c2f_dims_cg (2,map_coord(1),map_coord(2))
1669	bdims (2,2) = c2f_dims_cg (3,map_coord(1),map_coord(2))
1670	bdims (3,1) = c2f_dims_cg (4,map_coord(1),map_coord(2))
1671	bdims (3,2) = c2f_dims_cg (5,map_coord(1),map_coord(2))
1672
1673	n_cell_c = ( (bdims(1,2)-bdims(1,1)) + 3 ) * &
1674	( (bdims(2,2)-bdims(2,1)) + 3 ) * &
1675	( (bdims(3,2)-bdims(3,1)) + 1 )
1676
1677	CALL MPI_SEND(u (bdims(3,1), bdims(2,1)-1, bdims(1,1)-1), &
1678	1, TYPE_VNEST_BC, target_idex, &
1679	201, comm_inter, ierr)
1680
1681	CALL MPI_SEND(v(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
1682	1, TYPE_VNEST_BC, target_idex, &
1683	202, comm_inter, ierr)
1684
1685	CALL MPI_SEND(w(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
1686	1, TYPE_VNEST_BC, target_idex, &
1687	203, comm_inter, ierr)
1688
1689	CALL MPI_SEND(pt(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
1690	1, TYPE_VNEST_BC, target_idex, &
1691	205, comm_inter, ierr)
1692
1693	IF ( humidity ) THEN
1694	CALL MPI_SEND(q(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
1695	1, TYPE_VNEST_BC, target_idex, &
1696	209, comm_inter, ierr)
1697	ENDIF
1698
1699	end do
1700	end do
1701
1702	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
1703	!-- Receive data from coarse grid for TOP BC
1704
1705	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
1706	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
1707	map_coord(1) = offset(1)
1708	map_coord(2) = offset(2)
1709	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
1710
1711	bdims_rem (1,1) = c2f_dims_fg(0)
1712	bdims_rem (1,2) = c2f_dims_fg(1)
1713	bdims_rem (2,1) = c2f_dims_fg(2)
1714	bdims_rem (2,2) = c2f_dims_fg(3)
1715	bdims_rem (3,1) = c2f_dims_fg(4)
1716	bdims_rem (3,2) = c2f_dims_fg(5)
1717
1718	n_cell_c = &
1719	( (bdims_rem(1,2)-bdims_rem(1,1)) + 3 ) * &
1720	( (bdims_rem(2,2)-bdims_rem(2,1)) + 3 ) * &
1721	( (bdims_rem(3,2)-bdims_rem(3,1)) + 1 )
1722
1723	ALLOCATE( work3d ( &
1724	bdims_rem(3,1) :bdims_rem(3,2) , &
1725	bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
1726	bdims_rem(1,1)-1:bdims_rem(1,2)+1))
1727
1728
1729	CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 201, &
1730	comm_inter,status, ierr )
1731	interpol3d => u
1732	call vnest_set_topbc_u
1733
1734	CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 202, &
1735	comm_inter,status, ierr )
1736	interpol3d => v
1737	call vnest_set_topbc_v
1738
1739	CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 203, &
1740	comm_inter,status, ierr )
1741	interpol3d => w
1742	call vnest_set_topbc_w
1743
1744	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 205, &
1745	comm_inter,status, ierr )
1746	interpol3d => pt
1747	call vnest_set_topbc_s
1748
1749	IF ( humidity ) THEN
1750	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 209, &
1751	comm_inter,status, ierr )
1752	interpol3d => q
1753	call vnest_set_topbc_s
1754
1755	CALL exchange_horiz_2d(q (nzt+1,:,:) )
1756	ENDIF
1757
1758	!-- TKE Neumann BC for FG top
1759	DO jjf = nys, nyn
1760	DO iif = nxl, nxr
1761	e(nzt+1,jjf,iif) = e(nzt,jjf,iif)
1762	END DO
1763	END DO
1764
1765	!
1766	!-- w level nzt+1 does not impact results. Only to avoid jumps while
1767	!-- plotting profiles
1768	w(nzt+1,:,:) = w(nzt,:,:)
1769
1770	CALL exchange_horiz_2d(u (nzt+1,:,:) )
1771	CALL exchange_horiz_2d(v (nzt+1,:,:) )
1772	CALL exchange_horiz_2d(pt(nzt+1,:,:) )
1773	CALL exchange_horiz_2d(e (nzt+1,:,:) )
1774	CALL exchange_horiz_2d(w (nzt+1,:,:) )
1775	CALL exchange_horiz_2d(w (nzt ,:,:) )
1776
1777	NULLIFY ( interpol3d )
1778	DEALLOCATE ( work3d )
1779
1780	ENDIF
1781
1782
1783	CONTAINS
1784
1785	SUBROUTINE vnest_set_topbc_w
1786
1787
1788	USE arrays_3d
1789	USE control_parameters
1790	USE grid_variables
1791	USE indices
1792	USE pegrid
1793
1794
1795	IMPLICIT NONE
1796
1797	INTEGER(iwp) :: i
1798	INTEGER(iwp) :: j
1799	INTEGER(iwp) :: iif
1800	INTEGER(iwp) :: jjf
1801	INTEGER(iwp) :: bottomx
1802	INTEGER(iwp) :: bottomy
1803	INTEGER(iwp) :: topx
1804	INTEGER(iwp) :: topy
1805	REAL(wp) :: eps
1806	REAL(wp) :: alpha
1807	REAL(wp) :: eminus
1808	REAL(wp) :: edot
1809	REAL(wp) :: eplus
1810	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wprf
1811
1812
1813	ALLOCATE( wprf(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
1814
1815	!
1816	!-- Determination of a boundary condition for the vertical velocity component w:
1817	!-- In this case only interpolation in x- and y- direction is necessary, as the
1818	!-- boundary w-node of the fine grid coincides with a w-node in the coarse grid.
1819	!-- For both interpolations the scheme of Clark and Farley is used.
1820
1821	!
1822	!-- Interpolation in x-direction
1823
1824	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
1825
1826	DO i = bdims_rem(1,1), bdims_rem(1,2)
1827
1828	bottomx = (nxf+1)/(nxc+1) * i
1829	topx = (nxf+1)/(nxc+1) * (i+1) - 1
1830
1831	DO iif = bottomx, topx
1832
1833	eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
1834	alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
1835	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1836	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1837	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1838	wprf(j,iif) = eminus * work3d(bdims_rem(3,1),j,i-1) &
1839	+ edot * work3d(bdims_rem(3,1),j,i) &
1840	+ eplus * work3d(bdims_rem(3,1),j,i+1)
1841
1842	END DO
1843
1844	END DO
1845
1846	END DO
1847
1848	!
1849	!-- Interpolation in y-direction
1850
1851	DO j = bdims_rem(2,1), bdims_rem(2,2)
1852
1853	bottomy = (nyf+1)/(nyc+1) * j
1854	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1855
1856	DO iif = nxl, nxr
1857
1858	DO jjf = bottomy, topy
1859
1860	eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
1861
1862	alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
1863
1864	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1865
1866	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1867
1868	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1869
1870	w(nzt,jjf,iif) = eminus * wprf(j-1,iif) &
1871	+ edot * wprf(j,iif) &
1872	+ eplus * wprf(j+1,iif)
1873
1874	END DO
1875
1876	END DO
1877
1878	END DO
1879
1880	DEALLOCATE( wprf )
1881
1882	END SUBROUTINE vnest_set_topbc_w
1883
1884
1885	SUBROUTINE vnest_set_topbc_u
1886
1887
1888	USE arrays_3d
1889	USE control_parameters
1890	USE grid_variables
1891	USE indices
1892	USE pegrid
1893
1894
1895	IMPLICIT NONE
1896
1897	INTEGER(iwp) :: i
1898	INTEGER(iwp) :: j
1899	INTEGER(iwp) :: k
1900	INTEGER(iwp) :: iif
1901	INTEGER(iwp) :: jjf
1902	INTEGER(iwp) :: bottomx
1903	INTEGER(iwp) :: bottomy
1904	INTEGER(iwp) :: topx
1905	INTEGER(iwp) :: topy
1906	REAL(wp) :: eps
1907	REAL(wp) :: alpha
1908	REAL(wp) :: eminus
1909	REAL(wp) :: edot
1910	REAL(wp) :: eplus
1911	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: uprf
1912	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uprs
1913
1914	ALLOCATE( uprf(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
1915	ALLOCATE( uprs(nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
1916
1917
1918	!
1919	!-- Interpolation in y-direction
1920
1921	DO k = bdims_rem(3,1), bdims_rem(3,2)
1922	DO j = bdims_rem(2,1), bdims_rem(2,2)
1923
1924	bottomy = (nyf+1)/(nyc+1) * j
1925	topy = (nyf+1)/(nyc+1) * (j+1) - 1
1926
1927	DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
1928	DO jjf = bottomy, topy
1929
1930	eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
1931	alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
1932	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1933	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1934	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1935
1936	uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
1937	+ edot * work3d(k,j,i) &
1938	+ eplus * work3d(k,j+1,i)
1939	END DO
1940	END DO
1941
1942	END DO
1943	END DO
1944
1945	!
1946	!-- Interpolation in z-direction
1947
1948	DO jjf = nys, nyn
1949	DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
1950	eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
1951	alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
1952	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
1953	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
1954	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
1955	uprs(jjf,i) = eminus * uprf(bdims_rem(3,1),jjf,i) &
1956	+ edot * uprf(bdims_rem(3,1)+1,jjf,i) &
1957	+ eplus * uprf(bdims_rem(3,1)+2,jjf,i)
1958	END DO
1959	END DO
1960
1961	!
1962	!-- Interpolation in x-direction
1963
1964	DO jjf = nys, nyn
1965	DO i = bdims_rem(1,1), bdims_rem(1,2)
1966
1967	bottomx = (nxf+1)/(nxc+1) * i
1968	topx = (nxf+1)/(nxc+1) * (i+1) - 1
1969
1970	DO iif = bottomx, topx
1971	u(nzt+1,jjf,iif) = uprs(jjf,i) + ( iif * dxf - i * dxc ) * ( uprs(jjf,i+1) - uprs(jjf,i) ) / dxc
1972	END DO
1973
1974	END DO
1975	END DO
1976
1977
1978
1979	DEALLOCATE ( uprf, uprs )
1980
1981	END SUBROUTINE vnest_set_topbc_u
1982
1983
1984	SUBROUTINE vnest_set_topbc_v
1985
1986
1987	USE arrays_3d
1988	USE control_parameters
1989	USE grid_variables
1990	USE indices
1991	USE pegrid
1992
1993
1994	IMPLICIT NONE
1995
1996	INTEGER(iwp) :: i
1997	INTEGER(iwp) :: j
1998	INTEGER(iwp) :: k
1999	INTEGER(iwp) :: iif
2000	INTEGER(iwp) :: jjf
2001	INTEGER(iwp) :: bottomx
2002	INTEGER(iwp) :: bottomy
2003	INTEGER(iwp) :: topx
2004	INTEGER(iwp) :: topy
2005	REAL(wp) :: eps
2006	REAL(wp) :: alpha
2007	REAL(wp) :: eminus
2008	REAL(wp) :: edot
2009	REAL(wp) :: eplus
2010	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: vprf
2011	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vprs
2012
2013
2014
2015	ALLOCATE( vprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
2016	ALLOCATE( vprs(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
2017	!
2018	!-- Determination of a boundary condition for the horizontal velocity component v:
2019	!-- Interpolation in x- and z-direction is carried out by using the scheme,
2020	!-- which was derived by Clark and Farley (1984). In y-direction a
2021	!-- linear interpolation is carried out.
2022
2023	!
2024	!-- Interpolation in x-direction
2025
2026	DO k = bdims_rem(3,1), bdims_rem(3,2)
2027	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
2028	DO i = bdims_rem(1,1), bdims_rem(1,2)
2029
2030	bottomx = (nxf+1)/(nxc+1) * i
2031	topx = (nxf+1)/(nxc+1) * (i+1) - 1
2032
2033	DO iif = bottomx, topx
2034
2035	eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
2036	alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
2037	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
2038	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2039	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2040	vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
2041	+ edot * work3d(k,j,i) &
2042	+ eplus * work3d(k,j,i+1)
2043	END DO
2044
2045	END DO
2046	END DO
2047	END DO
2048
2049	!
2050	!-- Interpolation in z-direction
2051
2052	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
2053	DO iif = nxl, nxr
2054
2055	eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
2056	alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
2057	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
2058	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2059	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2060	vprs(j,iif) = eminus * vprf(bdims_rem(3,1),j,iif) &
2061	+ edot * vprf(bdims_rem(3,1)+1,j,iif) &
2062	+ eplus * vprf(bdims_rem(3,1)+2,j,iif)
2063
2064	END DO
2065	END DO
2066
2067	!
2068	!-- Interpolation in y-direction
2069
2070	DO j = bdims_rem(2,1), bdims_rem(2,2)
2071	DO iif = nxl, nxr
2072
2073	bottomy = (nyf+1)/(nyc+1) * j
2074	topy = (nyf+1)/(nyc+1) * (j+1) - 1
2075
2076	DO jjf = bottomy, topy
2077
2078	v(nzt+1,jjf,iif) = vprs(j,iif) + ( jjf * dyf - j * dyc ) * ( vprs(j+1,iif) - vprs(j,iif) ) / dyc
2079
2080	END DO
2081	END DO
2082	END DO
2083
2084
2085	DEALLOCATE ( vprf, vprs)
2086
2087
2088
2089	END SUBROUTINE vnest_set_topbc_v
2090
2091
2092	SUBROUTINE vnest_set_topbc_s
2093
2094
2095	USE arrays_3d
2096	USE control_parameters
2097	USE grid_variables
2098	USE indices
2099	USE pegrid
2100
2101
2102	IMPLICIT NONE
2103
2104	INTEGER(iwp) :: i
2105	INTEGER(iwp) :: j
2106	INTEGER(iwp) :: k
2107	INTEGER(iwp) :: iif
2108	INTEGER(iwp) :: jjf
2109	INTEGER(iwp) :: bottomx
2110	INTEGER(iwp) :: bottomy
2111	INTEGER(iwp) :: topx
2112	INTEGER(iwp) :: topy
2113	REAL(wp) :: eps
2114	REAL(wp) :: alpha
2115	REAL(wp) :: eminus
2116	REAL(wp) :: edot
2117	REAL(wp) :: eplus
2118	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
2119	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
2120
2121
2122
2123	ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
2124	ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
2125
2126	!
2127	!-- Determination of a boundary condition for the potential temperature pt:
2128	!-- The scheme derived by Clark and Farley can be used in all three dimensions.
2129
2130	!
2131	!-- Interpolation in x-direction
2132
2133	DO k = bdims_rem(3,1), bdims_rem(3,2)
2134
2135	DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
2136
2137	DO i = bdims_rem(1,1), bdims_rem(1,2)
2138
2139	bottomx = (nxf+1)/(nxc+1) * i
2140	topx = (nxf+1)/(nxc+1) *(i+1) - 1
2141
2142	DO iif = bottomx, topx
2143
2144	eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
2145
2146	alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
2147
2148	eminus = eps * (eps - 1.0 ) / 2.0 + alpha
2149
2150	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2151
2152	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2153
2154	ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
2155	+ edot * work3d(k,j,i) &
2156	+ eplus * work3d(k,j,i+1)
2157	END DO
2158
2159	END DO
2160
2161	END DO
2162
2163	END DO
2164
2165	!
2166	!-- Interpolation in y-direction
2167
2168	DO k = bdims_rem(3,1), bdims_rem(3,2)
2169
2170	DO j = bdims_rem(2,1), bdims_rem(2,2)
2171
2172	bottomy = (nyf+1)/(nyc+1) * j
2173	topy = (nyf+1)/(nyc+1) * (j+1) - 1
2174
2175	DO iif = nxl, nxr
2176
2177	DO jjf = bottomy, topy
2178
2179	eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
2180
2181	alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
2182
2183	eminus = eps * (eps - 1.0) / 2.0 + alpha
2184
2185	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2186
2187	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2188
2189	ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
2190	+ edot * ptprf(k,j,iif) &
2191	+ eplus * ptprf(k,j+1,iif)
2192	END DO
2193
2194	END DO
2195
2196	END DO
2197
2198	END DO
2199
2200	!
2201	!-- Interpolation in z-direction
2202
2203	DO jjf = nys, nyn
2204	DO iif = nxl, nxr
2205
2206	eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
2207
2208	alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
2209
2210	eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
2211
2212	edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2213
2214	eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2215
2216	interpol3d (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) &
2217	+ edot * ptprs(bdims_rem(3,1)+1,jjf,iif) &
2218	+ eplus * ptprs(bdims_rem(3,1)+2,jjf,iif)
2219
2220	END DO
2221	END DO
2222
2223	DEALLOCATE ( ptprf, ptprs )
2224
2225
2226
2227	END SUBROUTINE vnest_set_topbc_s
2228	#endif
2229	END SUBROUTINE vnest_boundary_conds
2230
2231
2232	SUBROUTINE vnest_boundary_conds_khkm
2233	#if defined( __parallel )
2234
2235	!--------------------------------------------------------------------------------!
2236	! Description:
2237	! ------------
2238	! Boundary conditions for the prognostic quantities.
2239	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
2240	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
2241	! handled in routine exchange_horiz. Pressure boundary conditions are
2242	! explicitly set in routines pres, poisfft, poismg and sor.
2243	!------------------------------------------------------------------------------!
2244
2245	USE arrays_3d
2246	USE control_parameters
2247	USE grid_variables
2248	USE indices
2249	USE pegrid
2250
2251
2252	IMPLICIT NONE
2253
2254	INTEGER(iwp) :: i
2255	INTEGER(iwp) :: j
2256	INTEGER(iwp) :: iif
2257	INTEGER(iwp) :: jjf
2258
2259
2260	IF ( coupling_mode == 'vnested_crse' ) THEN
2261	! Send data to fine grid for TOP BC
2262
2263	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
2264	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
2265
2266	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
2267	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
2268	map_coord(1) = i+offset(1)
2269	map_coord(2) = j+offset(2)
2270
2271	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
2272
2273	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 10, &
2274	comm_inter,status, ierr )
2275
2276	bdims (1,1) = bdims_rem (1,1) / cfratio(1)
2277	bdims (1,2) = bdims_rem (1,2) / cfratio(1)
2278	bdims (2,1) = bdims_rem (2,1) / cfratio(2)
2279	bdims (2,2) = bdims_rem (2,2) / cfratio(2)
2280	bdims (3,1) = bdims_rem (3,2) / cfratio(3)
2281	bdims (3,2) = bdims (3,1) + 2
2282
2283	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 9, &
2284	comm_inter, ierr )
2285
2286
2287	n_cell_c = ( (bdims(1,2)-bdims(1,1)) + 3 ) * &
2288	( (bdims(2,2)-bdims(2,1)) + 3 ) * &
2289	( (bdims(3,2)-bdims(3,1)) + 1 )
2290
2291	CALL MPI_SEND(kh(bdims(3,1) :bdims(3,2) , &
2292	bdims(2,1)-1:bdims(2,2)+1, &
2293	bdims(1,1)-1:bdims(1,2)+1),&
2294	n_cell_c, MPI_REAL, target_idex, &
2295	207, comm_inter, ierr)
2296
2297	CALL MPI_SEND(km(bdims(3,1) :bdims(3,2) , &
2298	bdims(2,1)-1:bdims(2,2)+1, &
2299	bdims(1,1)-1:bdims(1,2)+1),&
2300	n_cell_c, MPI_REAL, target_idex, &
2301	208, comm_inter, ierr)
2302
2303
2304
2305	end do
2306	end do
2307
2308	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
2309	! Receive data from coarse grid for TOP BC
2310
2311	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
2312	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
2313	map_coord(1) = offset(1)
2314	map_coord(2) = offset(2)
2315	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
2316
2317	bdims (1,1) = nxl
2318	bdims (1,2) = nxr
2319	bdims (2,1) = nys
2320	bdims (2,2) = nyn
2321	bdims (3,1) = nzb
2322	bdims (3,2) = nzt
2323
2324	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 10, &
2325	comm_inter, ierr )
2326
2327	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 9, &
2328	comm_inter,status, ierr )
2329
2330	n_cell_c = ( (bdims_rem(1,2)-bdims_rem(1,1)) + 3 ) * &
2331	( (bdims_rem(2,2)-bdims_rem(2,1)) + 3 ) * &
2332	( (bdims_rem(3,2)-bdims_rem(3,1)) + 1 )
2333
2334	ALLOCATE( work3d ( bdims_rem(3,1) :bdims_rem(3,2) , &
2335	bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
2336	bdims_rem(1,1)-1:bdims_rem(1,2)+1))
2337
2338
2339	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 207, &
2340	comm_inter,status, ierr )
2341
2342	! Neumann BC for FG kh
2343	DO jjf = nys, nyn
2344	DO iif = nxl, nxr
2345	kh(nzt+1,jjf,iif) = kh(nzt,jjf,iif)
2346	END DO
2347	END DO
2348
2349	CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 208, &
2350	comm_inter,status, ierr )
2351
2352	! Neumann BC for FG kh
2353	DO jjf = nys, nyn
2354	DO iif = nxl, nxr
2355	km(nzt+1,jjf,iif) = km(nzt,jjf,iif)
2356	END DO
2357	END DO
2358
2359
2360	!
2361	!-- The following evaluation can only be performed, if the fine grid is situated below the inversion
2362	!! DO jjf = nys-1, nyn+1
2363	!! DO iif = nxl-1, nxr+1
2364	!!
2365	!! km(nzt+1,jjf,iif) = 0.1 * l_grid(nzt+1) * SQRT( e(nzt+1,jjf,iif) )
2366	!! kh(nzt+1,jjf,iif) = 3.0 * km(nzt+1,jjf,iif)
2367	!!
2368	!! END DO
2369	!! END DO
2370
2371	CALL exchange_horiz_2d(km(nzt+1,:,:) )
2372	CALL exchange_horiz_2d(kh(nzt+1,:,:) )
2373
2374	DEALLOCATE ( work3d )
2375
2376	ENDIF
2377
2378
2379	! CONTAINS
2380	!
2381	! SUBROUTINE vnest_set_topbc_kh
2382	!
2383	!
2384	! USE arrays_3d
2385	! USE control_parameters
2386	! USE grid_variables
2387	! USE indices
2388	! USE pegrid
2389	!
2390	!
2391	! IMPLICIT NONE
2392	!
2393	! INTEGER(iwp) :: i
2394	! INTEGER(iwp) :: j
2395	! INTEGER(iwp) :: k
2396	! INTEGER(iwp) :: iif
2397	! INTEGER(iwp) :: jjf
2398	! INTEGER(iwp) :: bottomx
2399	! INTEGER(iwp) :: bottomy
2400	! INTEGER(iwp) :: topx
2401	! INTEGER(iwp) :: topy
2402	! REAL(wp) :: eps
2403	! REAL(wp) :: alpha
2404	! REAL(wp) :: eminus
2405	! REAL(wp) :: edot
2406	! REAL(wp) :: eplus
2407	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
2408	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
2409	!
2410	!
2411	!
2412	! ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
2413	! ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
2414	!
2415	! !
2416	! !-- Determination of a boundary condition for the potential temperature pt:
2417	! !-- The scheme derived by Clark and Farley can be used in all three dimensions.
2418	!
2419	! !
2420	! !-- Interpolation in x-direction
2421	!
2422	! DO k = bdims_rem(3,1), bdims_rem(3,2)
2423	!
2424	! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
2425	!
2426	! DO i = bdims_rem(1,1), bdims_rem(1,2)
2427	!
2428	! bottomx = (nxf+1)/(nxc+1) * i
2429	! topx = (nxf+1)/(nxc+1) *(i+1) - 1
2430	!
2431	! DO iif = bottomx, topx
2432	!
2433	! eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
2434	!
2435	! alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
2436	!
2437	! eminus = eps * (eps - 1.0 ) / 2.0 + alpha
2438	!
2439	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2440	!
2441	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2442	!
2443	! ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
2444	! + edot * work3d(k,j,i) &
2445	! + eplus * work3d(k,j,i+1)
2446	! END DO
2447	!
2448	! END DO
2449	!
2450	! END DO
2451	!
2452	! END DO
2453	!
2454	! !
2455	! !-- Interpolation in y-direction
2456	!
2457	! DO k = bdims_rem(3,1), bdims_rem(3,2)
2458	!
2459	! DO j = bdims_rem(2,1), bdims_rem(2,2)
2460	!
2461	! bottomy = (nyf+1)/(nyc+1) * j
2462	! topy = (nyf+1)/(nyc+1) * (j+1) - 1
2463	!
2464	! DO iif = nxl, nxr
2465	!
2466	! DO jjf = bottomy, topy
2467	!
2468	! eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
2469	!
2470	! alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
2471	!
2472	! eminus = eps * (eps - 1.0) / 2.0 + alpha
2473	!
2474	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2475	!
2476	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2477	!
2478	! ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
2479	! + edot * ptprf(k,j,iif) &
2480	! + eplus * ptprf(k,j+1,iif)
2481	! END DO
2482	!
2483	! END DO
2484	!
2485	! END DO
2486	!
2487	! END DO
2488	!
2489	! !
2490	! !-- Interpolation in z-direction
2491	!
2492	! DO jjf = nys, nyn
2493	! DO iif = nxl, nxr
2494	!
2495	! eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
2496	!
2497	! alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
2498	!
2499	! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
2500	!
2501	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2502	!
2503	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2504	!
2505	! kh (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) &
2506	! + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) &
2507	! + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif)
2508	!
2509	! END DO
2510	! END DO
2511	!
2512	! DEALLOCATE ( ptprf, ptprs )
2513	!
2514	!
2515	!
2516	! END SUBROUTINE vnest_set_topbc_kh
2517
2518	! SUBROUTINE vnest_set_topbc_km
2519	!
2520	!
2521	! USE arrays_3d
2522	! USE control_parameters
2523	! USE grid_variables
2524	! USE indices
2525	! USE pegrid
2526	!
2527	!
2528	! IMPLICIT NONE
2529	!
2530	! INTEGER(iwp) :: i
2531	! INTEGER(iwp) :: j
2532	! INTEGER(iwp) :: k
2533	! INTEGER(iwp) :: iif
2534	! INTEGER(iwp) :: jjf
2535	! INTEGER(iwp) :: bottomx
2536	! INTEGER(iwp) :: bottomy
2537	! INTEGER(iwp) :: topx
2538	! INTEGER(iwp) :: topy
2539	! REAL(wp) :: eps
2540	! REAL(wp) :: alpha
2541	! REAL(wp) :: eminus
2542	! REAL(wp) :: edot
2543	! REAL(wp) :: eplus
2544	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprf
2545	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptprs
2546	!
2547	!
2548	!
2549	! ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
2550	! ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
2551	!
2552	! !
2553	! !-- Determination of a boundary condition for the potential temperature pt:
2554	! !-- The scheme derived by Clark and Farley can be used in all three dimensions.
2555	!
2556	! !
2557	! !-- Interpolation in x-direction
2558	!
2559	! DO k = bdims_rem(3,1), bdims_rem(3,2)
2560	!
2561	! DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
2562	!
2563	! DO i = bdims_rem(1,1), bdims_rem(1,2)
2564	!
2565	! bottomx = (nxf+1)/(nxc+1) * i
2566	! topx = (nxf+1)/(nxc+1) *(i+1) - 1
2567	!
2568	! DO iif = bottomx, topx
2569	!
2570	! eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
2571	!
2572	! alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
2573	!
2574	! eminus = eps * (eps - 1.0 ) / 2.0 + alpha
2575	!
2576	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2577	!
2578	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2579	!
2580	! ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
2581	! + edot * work3d(k,j,i) &
2582	! + eplus * work3d(k,j,i+1)
2583	! END DO
2584	!
2585	! END DO
2586	!
2587	! END DO
2588	!
2589	! END DO
2590	!
2591	! !
2592	! !-- Interpolation in y-direction
2593	!
2594	! DO k = bdims_rem(3,1), bdims_rem(3,2)
2595	!
2596	! DO j = bdims_rem(2,1), bdims_rem(2,2)
2597	!
2598	! bottomy = (nyf+1)/(nyc+1) * j
2599	! topy = (nyf+1)/(nyc+1) * (j+1) - 1
2600	!
2601	! DO iif = nxl, nxr
2602	!
2603	! DO jjf = bottomy, topy
2604	!
2605	! eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
2606	!
2607	! alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
2608	!
2609	! eminus = eps * (eps - 1.0) / 2.0 + alpha
2610	!
2611	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2612	!
2613	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2614	!
2615	! ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
2616	! + edot * ptprf(k,j,iif) &
2617	! + eplus * ptprf(k,j+1,iif)
2618	! END DO
2619	!
2620	! END DO
2621	!
2622	! END DO
2623	!
2624	! END DO
2625	!
2626	! !
2627	! !-- Interpolation in z-direction
2628	!
2629	! DO jjf = nys, nyn
2630	! DO iif = nxl, nxr
2631	!
2632	! eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
2633	!
2634	! alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
2635	!
2636	! eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
2637	!
2638	! edot = ( 1.0 - eps*2.0 ) - 2.0 alpha
2639	!
2640	! eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
2641	!
2642	! km (nzt+1,jjf,iif) = eminus * ptprs(bdims_rem(3,1),jjf,iif) &
2643	! + edot * ptprs(bdims_rem(3,1)+1,jjf,iif) &
2644	! + eplus * ptprs(bdims_rem(3,1)+2,jjf,iif)
2645	!
2646	! END DO
2647	! END DO
2648	!
2649	! DEALLOCATE ( ptprf, ptprs )
2650	!
2651	!
2652	!
2653	! END SUBROUTINE vnest_set_topbc_km
2654
2655
2656	#endif
2657	END SUBROUTINE vnest_boundary_conds_khkm
2658
2659
2660
2661	SUBROUTINE vnest_anterpolate
2662
2663	#if defined( __parallel )
2664
2665	!--------------------------------------------------------------------------------!
2666	! Description:
2667	! ------------
2668	! Anterpolate data from fine grid to coarse grid.
2669	!------------------------------------------------------------------------------!
2670
2671	USE arrays_3d
2672	USE control_parameters
2673	USE grid_variables
2674	USE indices
2675	USE interfaces
2676	USE pegrid
2677	USE surface_mod, &
2678	ONLY : bc_h
2679
2680
2681	IMPLICIT NONE
2682
2683	REAL(wp) :: time_since_reference_point_rem
2684	INTEGER(iwp) :: i
2685	INTEGER(iwp) :: j
2686	INTEGER(iwp) :: k
2687	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
2688	INTEGER(iwp) :: m !< running index surface elements
2689
2690
2691
2692	!
2693	!-- In case of model termination initiated by the remote model
2694	!-- (terminate_coupled_remote > 0), initiate termination of the local model.
2695	!-- The rest of the coupler must then be skipped because it would cause an MPI
2696	!-- intercomminucation hang.
2697	!-- If necessary, the coupler will be called at the beginning of the next
2698	!-- restart run.
2699
2700	IF ( myid == 0) THEN
2701	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &
2702	target_id, 0, &
2703	terminate_coupled_remote, 1, MPI_INTEGER, &
2704	target_id, 0, &
2705	comm_inter, status, ierr )
2706	ENDIF
2707	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
2708	ierr )
2709
2710	IF ( terminate_coupled_remote > 0 ) THEN
2711	WRITE( message_string, * ) 'remote model "', &
2712	TRIM( coupling_mode_remote ), &
2713	'" terminated', &
2714	'&with terminate_coupled_remote = ', &
2715	terminate_coupled_remote, &
2716	'&local model "', TRIM( coupling_mode ), &
2717	'" has', &
2718	'&terminate_coupled = ', &
2719	terminate_coupled
2720	CALL message( 'vnest_anterpolate', 'PA0310', 1, 2, 0, 6, 0 )
2721	RETURN
2722	ENDIF
2723
2724
2725	!
2726	!-- Exchange the current simulated time between the models
2727
2728	IF ( myid == 0 ) THEN
2729
2730	CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
2731	11, comm_inter, ierr )
2732	CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
2733	target_id, 11, comm_inter, status, ierr )
2734
2735	ENDIF
2736
2737	CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
2738	ierr )
2739
2740	IF ( coupling_mode == 'vnested_crse' ) THEN
2741	! Receive data from fine grid for anterpolation
2742
2743	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
2744	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
2745
2746	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
2747	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
2748	map_coord(1) = i+offset(1)
2749	map_coord(2) = j+offset(2)
2750
2751	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
2752
2753	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 10, &
2754	comm_inter,status, ierr )
2755
2756	bdims (1,1) = bdims_rem (1,1) / cfratio(1)
2757	bdims (1,2) = bdims_rem (1,2) / cfratio(1)
2758	bdims (2,1) = bdims_rem (2,1) / cfratio(2)
2759	bdims (2,2) = bdims_rem (2,2) / cfratio(2)
2760	bdims (3,1) = bdims_rem (3,1)
2761	bdims (3,2) = bdims_rem (3,2) / cfratio(3)
2762
2763	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 9, &
2764	comm_inter, ierr )
2765
2766	n_cell_c = &
2767	(bdims(1,2)-bdims(1,1)+1) * &
2768	(bdims(2,2)-bdims(2,1)+1) * &
2769	(bdims(3,2)-bdims(3,1)+0)
2770
2771	CALL MPI_RECV( u( &
2772	bdims(3,1)+1:bdims(3,2), &
2773	bdims(2,1) :bdims(2,2), &
2774	bdims(1,1) :bdims(1,2)),&
2775	n_cell_c, MPI_REAL, target_idex, 101, &
2776	comm_inter,status, ierr )
2777
2778	CALL MPI_RECV( v( &
2779	bdims(3,1)+1:bdims(3,2), &
2780	bdims(2,1) :bdims(2,2), &
2781	bdims(1,1) :bdims(1,2)),&
2782	n_cell_c, MPI_REAL, target_idex, 102, &
2783	comm_inter,status, ierr )
2784
2785	CALL MPI_RECV(pt( &
2786	bdims(3,1)+1:bdims(3,2), &
2787	bdims(2,1) :bdims(2,2), &
2788	bdims(1,1) :bdims(1,2)),&
2789	n_cell_c, MPI_REAL, target_idex, 105, &
2790	comm_inter,status, ierr )
2791
2792	IF ( humidity ) THEN
2793	CALL MPI_RECV(q( &
2794	bdims(3,1)+1:bdims(3,2), &
2795	bdims(2,1) :bdims(2,2), &
2796	bdims(1,1) :bdims(1,2)),&
2797	n_cell_c, MPI_REAL, target_idex, 106, &
2798	comm_inter,status, ierr )
2799	ENDIF
2800
2801	CALL MPI_RECV( w( &
2802	bdims(3,1) :bdims(3,2)-1, &
2803	bdims(2,1) :bdims(2,2), &
2804	bdims(1,1) :bdims(1,2)), &
2805	n_cell_c, MPI_REAL, target_idex, 103, &
2806	comm_inter,status, ierr )
2807
2808	end do
2809	end do
2810
2811
2812
2813	!
2814	!-- Boundary conditions for the velocity components u and v
2815
2816
2817	IF ( ibc_uv_b == 0 ) THEN
2818	u(nzb,:,:) = 0.0_wp
2819	v(nzb,:,:) = 0.0_wp
2820	ELSE
2821	u(nzb,:,:) = u(nzb+1,:,:)
2822	v(nzb,:,:) = v(nzb+1,:,:)
2823	END IF
2824	!
2825	!-- Boundary conditions for the velocity components w
2826
2827	w(nzb,:,:) = 0.0_wp
2828
2829	!
2830	!-- Temperature at bottom boundary.
2831	!-- Neumann, zero-gradient
2832	IF ( ibc_pt_b == 1 ) THEN
2833	DO l = 0, 1
2834	DO m = 1, bc_h(l)%ns
2835	i = bc_h(l)%i(m)
2836	j = bc_h(l)%j(m)
2837	k = bc_h(l)%k(m)
2838	pt(k+bc_h(l)%koff,j,i) = pt(k,j,i)
2839	ENDDO
2840	ENDDO
2841	ENDIF
2842
2843
2844	CALL exchange_horiz( u, nbgp )
2845	CALL exchange_horiz( v, nbgp )
2846	CALL exchange_horiz( w, nbgp )
2847	CALL exchange_horiz( pt, nbgp )
2848
2849
2850	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
2851	! Send data to coarse grid for anterpolation
2852
2853	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
2854	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
2855	map_coord(1) = offset(1)
2856	map_coord(2) = offset(2)
2857	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
2858
2859	!-- Limit anterpolation level to nzt - z nesting ratio (a pseudo-buffer layer)
2860	bdims (1,1) = nxl
2861	bdims (1,2) = nxr
2862	bdims (2,1) = nys
2863	bdims (2,2) = nyn
2864	bdims (3,1) = nzb
2865	bdims (3,2) = nzt-cfratio(3)
2866
2867	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 10, &
2868	comm_inter, ierr )
2869
2870	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 9, &
2871	comm_inter,status, ierr )
2872
2873
2874	ALLOCATE( work3d ( &
2875	bdims_rem(3,1)+1:bdims_rem(3,2), &
2876	bdims_rem(2,1) :bdims_rem(2,2), &
2877	bdims_rem(1,1) :bdims_rem(1,2)))
2878
2879
2880	anterpol3d => u
2881
2882	CALL anterpolate_to_crse_u
2883	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
2884	101, comm_inter, ierr)
2885
2886	anterpol3d => v
2887
2888	CALL anterpolate_to_crse_v
2889	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
2890	102, comm_inter, ierr)
2891
2892	anterpol3d => pt
2893
2894	CALL anterpolate_to_crse_s
2895	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
2896	105, comm_inter, ierr)
2897
2898
2899	IF ( humidity ) THEN
2900
2901	anterpol3d => q
2902
2903	CALL anterpolate_to_crse_s
2904	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
2905	106, comm_inter, ierr)
2906	ENDIF
2907
2908
2909	DEALLOCATE( work3d )
2910	ALLOCATE( work3d ( bdims_rem(3,1) :bdims_rem(3,2)-1, &
2911	bdims_rem(2,1) :bdims_rem(2,2), &
2912	bdims_rem(1,1) :bdims_rem(1,2)))
2913	anterpol3d => w
2914	CALL anterpolate_to_crse_w
2915	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
2916	103, comm_inter, ierr)
2917
2918	NULLIFY ( anterpol3d )
2919	DEALLOCATE( work3d )
2920
2921	ENDIF
2922
2923
2924
2925	CONTAINS
2926	SUBROUTINE anterpolate_to_crse_u
2927
2928
2929	USE arrays_3d
2930	USE control_parameters
2931	USE grid_variables
2932	USE indices
2933	USE pegrid
2934
2935
2936	IMPLICIT NONE
2937
2938	INTEGER(iwp) :: i
2939	INTEGER(iwp) :: j
2940	INTEGER(iwp) :: k
2941	INTEGER(iwp) :: iif
2942	INTEGER(iwp) :: jjf
2943	INTEGER(iwp) :: kkf
2944	INTEGER(iwp) :: bottomy
2945	INTEGER(iwp) :: bottomz
2946	INTEGER(iwp) :: topy
2947	INTEGER(iwp) :: topz
2948	REAL(wp) :: aweight
2949
2950	!
2951	!-- Anterpolation of the velocity components u
2952	!-- only values in yz-planes that coincide in the fine and
2953	!-- the coarse grid are considered
2954
2955	DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
2956
2957	bottomz = (dzc/dzf) * (k-1) + 1
2958	topz = (dzc/dzf) * k
2959
2960	DO j = bdims_rem(2,1),bdims_rem(2,2)
2961
2962	bottomy = (nyf+1) / (nyc+1) * j
2963	topy = (nyf+1) / (nyc+1) * (j+1) - 1
2964
2965	DO i = bdims_rem(1,1),bdims_rem(1,2)
2966
2967	iif = (nxf+1) / (nxc+1) * i
2968
2969	aweight = 0.0
2970
2971	DO kkf = bottomz, topz
2972	DO jjf = bottomy, topy
2973
2974	aweight = aweight + anterpol3d(kkf,jjf,iif) * &
2975	(dzf/dzc) * (dyf/dyc)
2976
2977	END DO
2978	END DO
2979
2980	work3d(k,j,i) = aweight
2981
2982	END DO
2983
2984	END DO
2985
2986	END DO
2987
2988
2989
2990	END SUBROUTINE anterpolate_to_crse_u
2991
2992
2993	SUBROUTINE anterpolate_to_crse_v
2994
2995
2996	USE arrays_3d
2997	USE control_parameters
2998	USE grid_variables
2999	USE indices
3000	USE pegrid
3001
3002
3003	IMPLICIT NONE
3004
3005	INTEGER(iwp) :: i
3006	INTEGER(iwp) :: j
3007	INTEGER(iwp) :: k
3008	INTEGER(iwp) :: iif
3009	INTEGER(iwp) :: jjf
3010	INTEGER(iwp) :: kkf
3011	INTEGER(iwp) :: bottomx
3012	INTEGER(iwp) :: bottomz
3013	INTEGER(iwp) :: topx
3014	INTEGER(iwp) :: topz
3015	REAL(wp) :: aweight
3016
3017	!
3018	!-- Anterpolation of the velocity components v
3019	!-- only values in xz-planes that coincide in the fine and
3020	!-- the coarse grid are considered
3021
3022	DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
3023
3024	bottomz = (dzc/dzf) * (k-1) + 1
3025	topz = (dzc/dzf) * k
3026
3027	DO j = bdims_rem(2,1), bdims_rem(2,2)
3028
3029	jjf = (nyf+1) / (nyc+1) * j
3030
3031	DO i = bdims_rem(1,1), bdims_rem(1,2)
3032
3033	bottomx = (nxf+1) / (nxc+1) * i
3034	topx = (nxf+1) / (nxc+1) * (i+1) - 1
3035
3036	aweight = 0.0
3037
3038	DO kkf = bottomz, topz
3039	DO iif = bottomx, topx
3040
3041	aweight = aweight + anterpol3d(kkf,jjf,iif) * &
3042	(dzf/dzc) * (dxf/dxc)
3043
3044
3045	END DO
3046	END DO
3047
3048	work3d(k,j,i) = aweight
3049
3050	END DO
3051	END DO
3052	END DO
3053
3054
3055
3056	END SUBROUTINE anterpolate_to_crse_v
3057
3058
3059	SUBROUTINE anterpolate_to_crse_w
3060
3061
3062	USE arrays_3d
3063	USE control_parameters
3064	USE grid_variables
3065	USE indices
3066	USE pegrid
3067
3068
3069	IMPLICIT NONE
3070
3071	INTEGER(iwp) :: i
3072	INTEGER(iwp) :: j
3073	INTEGER(iwp) :: k
3074	INTEGER(iwp) :: iif
3075	INTEGER(iwp) :: jjf
3076	INTEGER(iwp) :: kkf
3077	INTEGER(iwp) :: bottomx
3078	INTEGER(iwp) :: bottomy
3079	INTEGER(iwp) :: topx
3080	INTEGER(iwp) :: topy
3081	REAL(wp) :: aweight
3082
3083	!
3084	!-- Anterpolation of the velocity components w
3085	!-- only values in xy-planes that coincide in the fine and
3086	!-- the coarse grid are considered
3087
3088	DO k = bdims_rem(3,1), bdims_rem(3,2)-1
3089
3090	kkf = cfratio(3) * k
3091
3092	DO j = bdims_rem(2,1), bdims_rem(2,2)
3093
3094	bottomy = (nyf+1) / (nyc+1) * j
3095	topy = (nyf+1) / (nyc+1) * (j+1) - 1
3096
3097	DO i = bdims_rem(1,1), bdims_rem(1,2)
3098
3099	bottomx = (nxf+1) / (nxc+1) * i
3100	topx = (nxf+1) / (nxc+1) * (i+1) - 1
3101
3102	aweight = 0.0
3103
3104	DO jjf = bottomy, topy
3105	DO iif = bottomx, topx
3106
3107	aweight = aweight + anterpol3d (kkf,jjf,iif) * &
3108	(dxf/dxc) * (dyf/dyc)
3109
3110	END DO
3111	END DO
3112
3113	work3d(k,j,i) = aweight
3114
3115	END DO
3116
3117	END DO
3118
3119	END DO
3120
3121
3122	END SUBROUTINE anterpolate_to_crse_w
3123
3124
3125	SUBROUTINE anterpolate_to_crse_s
3126
3127
3128	USE arrays_3d
3129	USE control_parameters
3130	USE grid_variables
3131	USE indices
3132	USE pegrid
3133
3134
3135	IMPLICIT NONE
3136
3137	INTEGER(iwp) :: i
3138	INTEGER(iwp) :: j
3139	INTEGER(iwp) :: k
3140	INTEGER(iwp) :: iif
3141	INTEGER(iwp) :: jjf
3142	INTEGER(iwp) :: kkf
3143	INTEGER(iwp) :: bottomx
3144	INTEGER(iwp) :: bottomy
3145	INTEGER(iwp) :: bottomz
3146	INTEGER(iwp) :: topx
3147	INTEGER(iwp) :: topy
3148	INTEGER(iwp) :: topz
3149	REAL(wp) :: aweight
3150
3151	!
3152	!-- Anterpolation of the potential temperature pt
3153	!-- all fine grid values are considered
3154
3155	DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
3156
3157	bottomz = (dzc/dzf) * (k-1) + 1
3158	topz = (dzc/dzf) * k
3159
3160	DO j = bdims_rem(2,1), bdims_rem(2,2)
3161
3162	bottomy = (nyf+1) / (nyc+1) * j
3163	topy = (nyf+1) / (nyc+1) * (j+1) - 1
3164
3165	DO i = bdims_rem(1,1), bdims_rem(1,2)
3166
3167	bottomx = (nxf+1) / (nxc+1) * i
3168	topx = (nxf+1) / (nxc+1) * (i+1) - 1
3169
3170	aweight = 0.0
3171
3172	DO kkf = bottomz, topz
3173	DO jjf = bottomy, topy
3174	DO iif = bottomx, topx
3175
3176	aweight = aweight + anterpol3d(kkf,jjf,iif) * &
3177	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3178
3179	END DO
3180	END DO
3181	END DO
3182
3183	work3d(k,j,i) = aweight
3184
3185	END DO
3186
3187	END DO
3188
3189	END DO
3190
3191
3192	END SUBROUTINE anterpolate_to_crse_s
3193	#endif
3194	END SUBROUTINE vnest_anterpolate
3195
3196
3197
3198	SUBROUTINE vnest_anterpolate_e
3199	#if defined( __parallel )
3200
3201	!--------------------------------------------------------------------------------!
3202	! Description:
3203	! ------------
3204	! Anterpolate TKE from fine grid to coarse grid.
3205	!------------------------------------------------------------------------------!
3206
3207	USE arrays_3d
3208	USE control_parameters
3209	USE grid_variables
3210	USE indices
3211	USE interfaces
3212	USE pegrid
3213
3214
3215	IMPLICIT NONE
3216
3217	REAL(wp) :: time_since_reference_point_rem
3218	INTEGER(iwp) :: i
3219	INTEGER(iwp) :: j
3220
3221	!
3222	!-- In case of model termination initiated by the remote model
3223	!-- (terminate_coupled_remote > 0), initiate termination of the local model.
3224	!-- The rest of the coupler must then be skipped because it would cause an MPI
3225	!-- intercomminucation hang.
3226	!-- If necessary, the coupler will be called at the beginning of the next
3227	!-- restart run.
3228
3229	IF ( myid == 0) THEN
3230	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &
3231	target_id, 0, &
3232	terminate_coupled_remote, 1, MPI_INTEGER, &
3233	target_id, 0, &
3234	comm_inter, status, ierr )
3235	ENDIF
3236	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
3237	ierr )
3238
3239	IF ( terminate_coupled_remote > 0 ) THEN
3240	WRITE( message_string, * ) 'remote model "', &
3241	TRIM( coupling_mode_remote ), &
3242	'" terminated', &
3243	'&with terminate_coupled_remote = ', &
3244	terminate_coupled_remote, &
3245	'&local model "', TRIM( coupling_mode ), &
3246	'" has', &
3247	'&terminate_coupled = ', &
3248	terminate_coupled
3249	CALL message( 'vnest_anterpolate_e', 'PA0310', 1, 2, 0, 6, 0 )
3250	RETURN
3251	ENDIF
3252
3253
3254	!
3255	!-- Exchange the current simulated time between the models
3256	IF ( myid == 0 ) THEN
3257
3258	CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
3259	11, comm_inter, ierr )
3260	CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
3261	target_id, 11, comm_inter, status, ierr )
3262
3263	ENDIF
3264
3265	CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
3266	ierr )
3267
3268	IF ( coupling_mode == 'vnested_crse' ) THEN
3269	! Receive data from fine grid for anterpolation
3270
3271	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
3272	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
3273
3274	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
3275	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
3276	map_coord(1) = i+offset(1)
3277	map_coord(2) = j+offset(2)
3278
3279	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
3280
3281	bdims (1,1) = f2c_dims_cg (0,map_coord(1),map_coord(2))
3282	bdims (1,2) = f2c_dims_cg (1,map_coord(1),map_coord(2))
3283	bdims (2,1) = f2c_dims_cg (2,map_coord(1),map_coord(2))
3284	bdims (2,2) = f2c_dims_cg (3,map_coord(1),map_coord(2))
3285	bdims (3,1) = f2c_dims_cg (4,map_coord(1),map_coord(2))
3286	bdims (3,2) = f2c_dims_cg (5,map_coord(1),map_coord(2))
3287
3288
3289	n_cell_c = (bdims(1,2)-bdims(1,1)+1) * &
3290	(bdims(2,2)-bdims(2,1)+1) * &
3291	(bdims(3,2)-bdims(3,1)+0)
3292
3293
3294	CALL MPI_RECV( e( bdims(3,1)+1:bdims(3,2), &
3295	bdims(2,1) :bdims(2,2), &
3296	bdims(1,1) :bdims(1,2)),&
3297	n_cell_c, MPI_REAL, target_idex, 104, &
3298	comm_inter,status, ierr )
3299	end do
3300	end do
3301
3302
3303	!
3304	!-- Boundary conditions
3305
3306	IF ( .NOT. constant_diffusion ) THEN
3307	e(nzb,:,:) = e(nzb+1,:,:)
3308	END IF
3309
3310	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e, nbgp )
3311
3312	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
3313	! Send data to coarse grid for anterpolation
3314
3315	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
3316	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
3317	map_coord(1) = offset(1)
3318	map_coord(2) = offset(2)
3319	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
3320
3321	bdims_rem (1,1) = f2c_dims_fg (0)
3322	bdims_rem (1,2) = f2c_dims_fg (1)
3323	bdims_rem (2,1) = f2c_dims_fg (2)
3324	bdims_rem (2,2) = f2c_dims_fg (3)
3325	bdims_rem (3,1) = f2c_dims_fg (4)
3326	bdims_rem (3,2) = f2c_dims_fg (5)
3327
3328	ALLOCATE( work3d ( &
3329	bdims_rem(3,1)+1:bdims_rem(3,2), &
3330	bdims_rem(2,1) :bdims_rem(2,2), &
3331	bdims_rem(1,1) :bdims_rem(1,2)))
3332
3333	anterpol3d => e
3334
3335	CALL anterpolate_to_crse_e
3336
3337	CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex, &
3338	104, comm_inter, ierr)
3339
3340	NULLIFY ( anterpol3d )
3341	DEALLOCATE( work3d )
3342	ENDIF
3343
3344
3345	CONTAINS
3346
3347
3348
3349
3350
3351	SUBROUTINE anterpolate_to_crse_e
3352
3353
3354	USE arrays_3d
3355	USE control_parameters
3356	USE grid_variables
3357	USE indices
3358	USE pegrid
3359
3360
3361	IMPLICIT NONE
3362
3363	INTEGER(iwp) :: i
3364	INTEGER(iwp) :: j
3365	INTEGER(iwp) :: k
3366	INTEGER(iwp) :: iif
3367	INTEGER(iwp) :: jjf
3368	INTEGER(iwp) :: kkf
3369	INTEGER(iwp) :: bottomx
3370	INTEGER(iwp) :: bottomy
3371	INTEGER(iwp) :: bottomz
3372	INTEGER(iwp) :: topx
3373	INTEGER(iwp) :: topy
3374	INTEGER(iwp) :: topz
3375	REAL(wp) :: aweight_a
3376	REAL(wp) :: aweight_b
3377	REAL(wp) :: aweight_c
3378	REAL(wp) :: aweight_d
3379	REAL(wp) :: aweight_e
3380	REAL(wp) :: energ
3381
3382	DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
3383
3384	bottomz = (dzc/dzf) * (k-1) + 1
3385	topz = (dzc/dzf) * k
3386
3387	DO j = bdims_rem(2,1), bdims_rem(2,2)
3388
3389	bottomy = (nyf+1) / (nyc+1) * j
3390	topy = (nyf+1) / (nyc+1) * (j+1) - 1
3391
3392	DO i = bdims_rem(1,1), bdims_rem(1,2)
3393
3394	bottomx = (nxf+1) / (nxc+1) * i
3395	topx = (nxf+1) / (nxc+1) * (i+1) - 1
3396
3397	aweight_a = 0.0
3398	aweight_b = 0.0
3399	aweight_c = 0.0
3400	aweight_d = 0.0
3401	aweight_e = 0.0
3402
3403	DO kkf = bottomz, topz
3404	DO jjf = bottomy, topy
3405	DO iif = bottomx, topx
3406
3407	aweight_a = aweight_a + anterpol3d(kkf,jjf,iif) * &
3408	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3409
3410
3411	energ = ( 0.5 * ( u(kkf,jjf,iif) + u(kkf,jjf,iif+1) ) )**2.0 + &
3412	( 0.5 * ( v(kkf,jjf,iif) + v(kkf,jjf+1,iif) ) )**2.0 + &
3413	( 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) )**2.0
3414
3415	aweight_b = aweight_b + energ * &
3416	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3417
3418	aweight_c = aweight_c + 0.5 * ( u(kkf,jjf,iif) + u(kkf,jjf,iif+1) ) * &
3419	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3420
3421	aweight_d = aweight_d + 0.5 * ( v(kkf,jjf,iif) + v(kkf,jjf+1,iif) ) * &
3422	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3423
3424	aweight_e = aweight_e + 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) * &
3425	(dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3426
3427
3428	END DO
3429	END DO
3430	END DO
3431
3432	work3d(k,j,i) = aweight_a + 0.5 * ( aweight_b - &
3433	aweight_c**2.0 - &
3434	aweight_d**2.0 - &
3435	aweight_e**2.0 )
3436
3437	END DO
3438
3439	END DO
3440
3441	END DO
3442
3443
3444
3445	END SUBROUTINE anterpolate_to_crse_e
3446	#endif
3447	END SUBROUTINE vnest_anterpolate_e
3448
3449	SUBROUTINE vnest_init_pegrid_rank
3450	#if defined( __parallel )
3451	! Domain decomposition and exchange of grid variables between coarse and fine
3452	! Given processor coordinates as index f_rnk_lst(pcoord(1), pcoord(2))
3453	! returns the rank. A single coarse block will have to send data to multiple
3454	! fine blocks. In the coarse grid the pcoords of the remote block is first found and then using
3455	! f_rnk_lst the target_idex is identified.
3456	! blk_dim stores the index limits of a given block. blk_dim_remote is received
3457	! from the asscoiated nest partner.
3458	! cf_ratio(1:3) is the ratio between fine and coarse grid: nxc/nxf, nyc/nyf and
3459	! ceiling(dxc/dxf)
3460
3461
3462	USE control_parameters, &
3463	ONLY: coupling_mode
3464
3465	USE kinds
3466
3467	USE pegrid
3468
3469
3470	IMPLICIT NONE
3471
3472	INTEGER(iwp) :: dest_rnk
3473	INTEGER(iwp) :: i !<
3474
3475	IF (myid == 0) THEN
3476	IF ( coupling_mode == 'vnested_crse') THEN
3477	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 33, comm_inter, &
3478	ierr )
3479	CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, numprocs, 66, &
3480	comm_inter, status, ierr )
3481	ELSEIF ( coupling_mode == 'vnested_fine') THEN
3482	CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 33, &
3483	comm_inter, status, ierr )
3484	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 66, comm_inter, &
3485	ierr )
3486	ENDIF
3487	ENDIF
3488
3489
3490	IF ( coupling_mode == 'vnested_crse') THEN
3491	CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3492	ALLOCATE( c_rnk_lst( 0:(pdims(1)-1) ,0:(pdims(2)-1) ) )
3493	ALLOCATE( f_rnk_lst( 0:(pdims_partner(1)-1) ,0:(pdims_partner(2)-1) ) )
3494	do i=0,numprocs-1
3495	CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
3496	call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
3497	c_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
3498	end do
3499	ELSEIF ( coupling_mode == 'vnested_fine') THEN
3500	CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3501	ALLOCATE( c_rnk_lst( 0:(pdims_partner(1)-1) ,0:(pdims_partner(2)-1) ) )
3502	ALLOCATE( f_rnk_lst( 0:(pdims(1)-1) ,0:(pdims(2)-1) ) )
3503
3504	do i=0,numprocs-1
3505	CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
3506	call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
3507	f_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
3508	enddo
3509	ENDIF
3510
3511
3512	IF ( coupling_mode == 'vnested_crse') THEN
3513	if (myid == 0) then
3514	CALL MPI_SEND( c_rnk_lst, pdims(1)*pdims(2), MPI_INTEGER, numprocs, 0, comm_inter, ierr )
3515	CALL MPI_RECV( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, numprocs, 4, comm_inter,status, ierr )
3516	end if
3517	CALL MPI_BCAST( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
3518	ELSEIF ( coupling_mode == 'vnested_fine') THEN
3519	if (myid == 0) then
3520	CALL MPI_RECV( c_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, 0, comm_inter,status, ierr )
3521	CALL MPI_SEND( f_rnk_lst, pdims(1)*pdims(2), MPI_INTEGER, 0, 4, comm_inter, ierr )
3522	end if
3523	CALL MPI_BCAST( c_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
3524	ENDIF
3525
3526	!-- Reason for MPI error unknown; solved if three lines duplicated
3527	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
3528	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
3529	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
3530
3531
3532	#endif
3533
3534	END SUBROUTINE vnest_init_pegrid_rank
3535
3536
3537	SUBROUTINE vnest_init_pegrid_domain
3538	#if defined( __parallel )
3539
3540	USE control_parameters, &
3541	ONLY: coupling_mode, coupling_topology, dz, &
3542	dz_stretch_level_start, message_string
3543
3544	USE grid_variables, &
3545	ONLY: dx, dy
3546
3547	USE indices, &
3548	ONLY: nbgp, nx, ny, nz, nxl, nxr, nys, nyn, nzb, nzt
3549
3550	USE kinds
3551
3552	USE pegrid
3553
3554	IMPLICIT NONE
3555
3556	INTEGER(iwp) :: i !<
3557	INTEGER(iwp) :: j !<
3558	INTEGER(iwp) :: tempx
3559	INTEGER(iwp) :: tempy
3560	INTEGER(iwp) :: TYPE_INT_YZ
3561	INTEGER(iwp) :: SIZEOFREAL
3562	INTEGER(iwp) :: MTV_X
3563	INTEGER(iwp) :: MTV_Y
3564	INTEGER(iwp) :: MTV_Z
3565	INTEGER(iwp) :: MTV_RX
3566	INTEGER(iwp) :: MTV_RY
3567	INTEGER(iwp) :: MTV_RZ
3568
3569	!
3570	!-- Pass the number of grid points of the coarse model to
3571	!-- the nested model and vice versa
3572	IF ( coupling_mode == 'vnested_crse' ) THEN
3573
3574	nxc = nx
3575	nyc = ny
3576	nzc = nz
3577	dxc = dx
3578	dyc = dy
3579	dzc = dz(1)
3580	cg_nprocs = numprocs
3581
3582	IF ( myid == 0 ) THEN
3583
3584	CALL MPI_SEND( nxc, 1, MPI_INTEGER , numprocs, 1, comm_inter, &
3585	ierr )
3586	CALL MPI_SEND( nyc, 1, MPI_INTEGER , numprocs, 2, comm_inter, &
3587	ierr )
3588	CALL MPI_SEND( nzc, 1, MPI_INTEGER , numprocs, 3, comm_inter, &
3589	ierr )
3590	CALL MPI_SEND( dxc, 1, MPI_REAL , numprocs, 4, comm_inter, &
3591	ierr )
3592	CALL MPI_SEND( dyc, 1, MPI_REAL , numprocs, 5, comm_inter, &
3593	ierr )
3594	CALL MPI_SEND( dzc, 1, MPI_REAL , numprocs, 6, comm_inter, &
3595	ierr )
3596	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7, comm_inter, &
3597	ierr )
3598	CALL MPI_SEND( cg_nprocs, 1, MPI_INTEGER, numprocs, 8, comm_inter, &
3599	ierr )
3600	CALL MPI_RECV( nxf, 1, MPI_INTEGER, numprocs, 21, comm_inter, &
3601	status, ierr )
3602	CALL MPI_RECV( nyf, 1, MPI_INTEGER, numprocs, 22, comm_inter, &
3603	status, ierr )
3604	CALL MPI_RECV( nzf, 1, MPI_INTEGER, numprocs, 23, comm_inter, &
3605	status, ierr )
3606	CALL MPI_RECV( dxf, 1, MPI_REAL, numprocs, 24, comm_inter, &
3607	status, ierr )
3608	CALL MPI_RECV( dyf, 1, MPI_REAL, numprocs, 25, comm_inter, &
3609	status, ierr )
3610	CALL MPI_RECV( dzf, 1, MPI_REAL, numprocs, 26, comm_inter, &
3611	status, ierr )
3612	CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, &
3613	numprocs, 27, comm_inter, status, ierr )
3614	CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER, &
3615	numprocs, 28, comm_inter, status, ierr )
3616	ENDIF
3617
3618	CALL MPI_BCAST( nxf, 1, MPI_INTEGER, 0, comm2d, ierr )
3619	CALL MPI_BCAST( nyf, 1, MPI_INTEGER, 0, comm2d, ierr )
3620	CALL MPI_BCAST( nzf, 1, MPI_INTEGER, 0, comm2d, ierr )
3621	CALL MPI_BCAST( dxf, 1, MPI_REAL, 0, comm2d, ierr )
3622	CALL MPI_BCAST( dyf, 1, MPI_REAL, 0, comm2d, ierr )
3623	CALL MPI_BCAST( dzf, 1, MPI_REAL, 0, comm2d, ierr )
3624	CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3625	CALL MPI_BCAST( fg_nprocs, 1, MPI_INTEGER, 0, comm2d, ierr )
3626
3627	!
3628	!-- Check if stretching is used within the nested domain. ABS(...) is
3629	!-- necessary because of the default value of -9999999.9_wp (negative)
3630	IF ( ABS( dz_stretch_level_start(1) ) <= (nzf+1)*dzf ) THEN
3631	message_string = 'Stretching in the parent domain is '// &
3632	'only allowed above the nested domain'
3633	CALL message( 'vnest_init_pegrid_domain', 'PA0497', 1, 2, 0, 6, 0 )
3634	ENDIF
3635
3636	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
3637
3638	nxf = nx
3639	nyf = ny
3640	nzf = nz
3641	dxf = dx
3642	dyf = dy
3643	dzf = dz(1)
3644	fg_nprocs = numprocs
3645
3646	IF ( myid == 0 ) THEN
3647
3648	CALL MPI_RECV( nxc, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
3649	ierr )
3650	CALL MPI_RECV( nyc, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
3651	ierr )
3652	CALL MPI_RECV( nzc, 1, MPI_INTEGER, 0, 3, comm_inter, status, &
3653	ierr )
3654	CALL MPI_RECV( dxc, 1, MPI_REAL, 0, 4, comm_inter, status, &
3655	ierr )
3656	CALL MPI_RECV( dyc, 1, MPI_REAL, 0, 5, comm_inter, status, &
3657	ierr )
3658	CALL MPI_RECV( dzc, 1, MPI_REAL, 0, 6, comm_inter, status, &
3659	ierr )
3660	CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 7, comm_inter, &
3661	status, ierr )
3662	CALL MPI_RECV( cg_nprocs, 1, MPI_INTEGER, 0, 8, comm_inter, &
3663	status, ierr )
3664	CALL MPI_SEND( nxf, 1, MPI_INTEGER, 0, 21, comm_inter, ierr )
3665	CALL MPI_SEND( nyf, 1, MPI_INTEGER, 0, 22, comm_inter, ierr )
3666	CALL MPI_SEND( nzf, 1, MPI_INTEGER, 0, 23, comm_inter, ierr )
3667	CALL MPI_SEND( dxf, 1, MPI_REAL, 0, 24, comm_inter, ierr )
3668	CALL MPI_SEND( dyf, 1, MPI_REAL, 0, 25, comm_inter, ierr )
3669	CALL MPI_SEND( dzf, 1, MPI_REAL, 0, 26, comm_inter, ierr )
3670	CALL MPI_SEND( pdims,2,MPI_INTEGER, 0, 27, comm_inter, ierr )
3671	CALL MPI_SEND( fg_nprocs,1,MPI_INTEGER, 0, 28, comm_inter, ierr )
3672	ENDIF
3673
3674	CALL MPI_BCAST( nxc, 1, MPI_INTEGER, 0, comm2d, ierr)
3675	CALL MPI_BCAST( nyc, 1, MPI_INTEGER, 0, comm2d, ierr)
3676	CALL MPI_BCAST( nzc, 1, MPI_INTEGER, 0, comm2d, ierr)
3677	CALL MPI_BCAST( dxc, 1, MPI_REAL, 0, comm2d, ierr)
3678	CALL MPI_BCAST( dyc, 1, MPI_REAL, 0, comm2d, ierr)
3679	CALL MPI_BCAST( dzc, 1, MPI_REAL, 0, comm2d, ierr)
3680	CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr)
3681	CALL MPI_BCAST( cg_nprocs, 1, MPI_INTEGER, 0, comm2d, ierr)
3682
3683	ENDIF
3684
3685	ngp_c = ( nxc+1 + 2 * nbgp ) * ( nyc+1 + 2 * nbgp )
3686	ngp_f = ( nxf+1 + 2 * nbgp ) * ( nyf+1 + 2 * nbgp )
3687
3688	IF ( coupling_mode(1:8) == 'vnested_') coupling_topology = 1
3689
3690
3691	!-- Nesting Ratio: For each coarse grid cell how many fine grid cells exist
3692	cfratio(1) = INT ( (nxf+1) / (nxc+1) )
3693	cfratio(2) = INT ( (nyf+1) / (nyc+1) )
3694	cfratio(3) = CEILING ( dzc / dzf )
3695
3696	!-- target_id is used only for exhange of information like simulated_time
3697	!-- which are then MPI_BCAST to other processors in the group
3698	IF ( myid == 0 ) THEN
3699
3700	IF ( TRIM( coupling_mode ) == 'vnested_crse' ) THEN
3701	target_id = numprocs
3702	ELSE IF ( TRIM( coupling_mode ) == 'vnested_fine' ) THEN
3703	target_id = 0
3704	ENDIF
3705
3706	ENDIF
3707
3708	!-- Store partner grid dimenstions and create MPI derived types
3709
3710	IF ( coupling_mode == 'vnested_crse' ) THEN
3711
3712	offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
3713	offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
3714
3715	tempx = ( pdims_partner(1) / pdims(1) ) - 1
3716	tempy = ( pdims_partner(2) / pdims(2) ) - 1
3717	ALLOCATE( c2f_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
3718	ALLOCATE( f2c_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
3719
3720	do j = 0, ( pdims_partner(2) / pdims(2) ) - 1
3721	do i = 0, ( pdims_partner(1) / pdims(1) ) - 1
3722	map_coord(1) = i+offset(1)
3723	map_coord(2) = j+offset(2)
3724
3725	target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
3726
3727	CALL MPI_RECV( bdims_rem, 6, MPI_INTEGER, target_idex, 10, &
3728	comm_inter,status, ierr )
3729
3730	!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
3731	!-- One CG can have multiple FG partners. The 3D array is mapped by partner proc co-ord
3732	c2f_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
3733	c2f_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
3734	c2f_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
3735	c2f_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
3736	c2f_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,2) / cfratio(3)
3737	c2f_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2) / cfratio(3)) + 2
3738
3739	!-- Store the CG dimensions that correspond to the FG partner; needed for anterpolation
3740	f2c_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
3741	f2c_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
3742	f2c_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
3743	f2c_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
3744	f2c_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,1)
3745	f2c_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2)-cfratio(3))/ cfratio(3)
3746
3747	CALL MPI_SEND( c2f_dims_cg (:,map_coord(1),map_coord(2)), 6, &
3748	MPI_INTEGER, target_idex, 100, comm_inter, ierr )
3749
3750	CALL MPI_SEND( f2c_dims_cg (:,map_coord(1),map_coord(2)), 6, &
3751	MPI_INTEGER, target_idex, 101, comm_inter, ierr )
3752
3753	end do
3754	end do
3755
3756	!-- A derived data type to pack 3 Z-levels of CG to set FG top BC
3757	MTV_X = ( nxr - nxl + 1 ) + 2*nbgp
3758	MTV_Y = ( nyn - nys + 1 ) + 2*nbgp
3759	MTV_Z = nzt+1 - nzb +1
3760
3761	MTV_RX = ( c2f_dims_cg (1,offset(1),offset(2)) - c2f_dims_cg (0,offset(1),offset(2)) ) +1+2
3762	MTV_RY = ( c2f_dims_cg (3,offset(1),offset(2)) - c2f_dims_cg (2,offset(1),offset(2)) ) +1+2
3763	MTV_RZ = ( c2f_dims_cg (5,offset(1),offset(2)) - c2f_dims_cg (4,offset(1),offset(2)) ) +1
3764
3765	CALL MPI_TYPE_EXTENT(MPI_REAL, SIZEOFREAL, IERR)
3766
3767	CALL MPI_TYPE_VECTOR ( MTV_RY, MTV_RZ, MTV_Z, MPI_REAL, TYPE_INT_YZ, IERR)
3768	CALL MPI_TYPE_HVECTOR( MTV_RX, 1, MTV_ZMTV_YSIZEOFREAL, &
3769	TYPE_INT_YZ, TYPE_VNEST_BC, IERR)
3770	CALL MPI_TYPE_FREE(TYPE_INT_YZ, IERR)
3771	CALL MPI_TYPE_COMMIT(TYPE_VNEST_BC, IERR)
3772
3773
3774	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
3775
3776	ALLOCATE( c2f_dims_fg (0:5) )
3777	ALLOCATE( f2c_dims_fg (0:5) )
3778
3779	offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
3780	offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
3781	map_coord(1) = offset(1)
3782	map_coord(2) = offset(2)
3783	target_idex = c_rnk_lst(map_coord(1),map_coord(2))
3784
3785	bdims (1,1) = nxl
3786	bdims (1,2) = nxr
3787	bdims (2,1) = nys
3788	bdims (2,2) = nyn
3789	bdims (3,1) = nzb
3790	bdims (3,2) = nzt
3791
3792	CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 10, &
3793	comm_inter, ierr )
3794
3795	!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
3796	!-- One FG can have only one CG partner
3797	CALL MPI_RECV( c2f_dims_fg, 6, MPI_INTEGER, target_idex, 100, &
3798	comm_inter,status, ierr )
3799
3800	CALL MPI_RECV( f2c_dims_fg, 6, MPI_INTEGER, target_idex, 101, &
3801	comm_inter,status, ierr )
3802
3803	!-- Store the CG dimensions that correspond to the FG partner; needed for anterpolation
3804
3805	n_cell_c = (f2c_dims_fg(1)-f2c_dims_fg(0)+1) * &
3806	(f2c_dims_fg(3)-f2c_dims_fg(2)+1) * &
3807	(f2c_dims_fg(5)-f2c_dims_fg(4)+0)
3808
3809	CALL MPI_TYPE_CONTIGUOUS(n_cell_c, MPI_REAL, TYPE_VNEST_ANTER, IERR)
3810	CALL MPI_TYPE_COMMIT(TYPE_VNEST_ANTER, ierr)
3811
3812	ENDIF
3813	#endif
3814	END SUBROUTINE vnest_init_pegrid_domain
3815
3816
3817	SUBROUTINE vnest_init_grid
3818
3819	#if defined( __parallel )
3820	USE arrays_3d, &
3821	ONLY: zu, zw
3822
3823	USE control_parameters, &
3824	ONLY: coupling_mode, message_string, number_stretch_level_start
3825
3826	USE indices, &
3827	ONLY: nzt
3828
3829	USE kinds
3830
3831	USE pegrid
3832
3833	IMPLICIT NONE
3834
3835	!
3836	!-- Allocate and Exchange zuc and zuf, zwc and zwf
3837	IF ( coupling_mode(1:8) == 'vnested_' ) THEN
3838
3839	ALLOCATE( zuc(0:nzc+1), zuf(0:nzf+1) )
3840	ALLOCATE( zwc(0:nzc+1), zwf(0:nzf+1) )
3841
3842	IF ( coupling_mode == 'vnested_crse' ) THEN
3843
3844	zuc = zu
3845	zwc = zw
3846
3847	IF ( myid == 0 ) THEN
3848
3849	CALL MPI_SEND( zuc, nzt+2, MPI_REAL, numprocs, 41, comm_inter, &
3850	ierr )
3851	CALL MPI_RECV( zuf, nzf+2, MPI_REAL, numprocs, 42, comm_inter, &
3852	status, ierr )
3853
3854	CALL MPI_SEND( zwc, nzt+2, MPI_REAL, numprocs, 43, comm_inter, &
3855	ierr )
3856	CALL MPI_RECV( zwf, nzf+2, MPI_REAL, numprocs, 44, comm_inter, &
3857	status, ierr )
3858
3859	ENDIF
3860
3861	CALL MPI_BCAST( zuf,nzf+2,MPI_REAL, 0, comm2d, ierr )
3862	CALL MPI_BCAST( zwf,nzf+2,MPI_REAL, 0, comm2d, ierr )
3863
3864	ELSEIF ( coupling_mode == 'vnested_fine' ) THEN
3865
3866	!
3867	!-- Check if stretching is used within the nested domain
3868	IF ( number_stretch_level_start > 0 ) THEN
3869	message_string = 'Stretching in the nested domain is not '//&
3870	'allowed'
3871	CALL message( 'vnest_init_grid', 'PA0498', 1, 2, 0, 6, 0 )
3872	ENDIF
3873
3874	zuf = zu
3875	zwf = zw
3876
3877	IF ( myid == 0 ) THEN
3878
3879	CALL MPI_RECV( zuc,nzc+2, MPI_REAL, 0, 41, comm_inter, status, &
3880	ierr )
3881	CALL MPI_SEND( zuf,nzt+2, MPI_REAL, 0, 42, comm_inter, ierr )
3882
3883	CALL MPI_RECV( zwc,nzc+2, MPI_REAL, 0, 43, comm_inter, status, &
3884	ierr )
3885	CALL MPI_SEND( zwf,nzt+2, MPI_REAL, 0, 44, comm_inter, ierr )
3886	ENDIF
3887
3888	CALL MPI_BCAST( zuc,nzc+2,MPI_REAL, 0, comm2d, ierr )
3889	CALL MPI_BCAST( zwc,nzc+2,MPI_REAL, 0, comm2d, ierr )
3890
3891	ENDIF
3892	ENDIF
3893
3894	#endif
3895	END SUBROUTINE vnest_init_grid
3896
3897
3898	SUBROUTINE vnest_check_parameters
3899	#if defined( __parallel )
3900
3901	USE pegrid, &
3902	ONLY: myid
3903
3904	IMPLICIT NONE
3905
3906	IF (myid==0) PRINT, ' vnest: check parameters not implemented yet *'
3907
3908	#endif
3909	END SUBROUTINE vnest_check_parameters
3910
3911
3912	SUBROUTINE vnest_timestep_sync
3913
3914	#if defined( __parallel )
3915	USE control_parameters, &
3916	ONLY: coupling_mode, dt_3d
3917
3918	USE interfaces
3919
3920	USE kinds
3921
3922	USE pegrid
3923
3924	IMPLICIT NONE
3925
3926	IF ( coupling_mode == 'vnested_crse') THEN
3927	dtc = dt_3d
3928	if (myid == 0) then
3929	CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id, &
3930	31, comm_inter, ierr )
3931	CALL MPI_RECV( dtf, 1, MPI_REAL, &
3932	target_id, 32, comm_inter, status, ierr )
3933
3934	endif
3935	CALL MPI_BCAST( dtf, 1, MPI_REAL, 0, comm2d, ierr )
3936	ELSE
3937	dtf = dt_3d
3938	if (myid == 0) then
3939	CALL MPI_RECV( dtc, 1, MPI_REAL, &
3940	target_id, 31, comm_inter, status, ierr )
3941	CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id, &
3942	32, comm_inter, ierr )
3943
3944	endif
3945	CALL MPI_BCAST( dtc, 1, MPI_REAL, 0, comm2d, ierr )
3946
3947	ENDIF
3948	!-- Identical timestep for coarse and fine grids
3949	dt_3d = MIN( dtc, dtf )
3950	#endif
3951	END SUBROUTINE vnest_timestep_sync
3952
3953	SUBROUTINE vnest_deallocate
3954	#if defined( __parallel )
3955	USE control_parameters, &
3956	ONLY: coupling_mode
3957
3958	IMPLICIT NONE
3959
3960	IF ( ALLOCATED(c_rnk_lst) ) DEALLOCATE (c_rnk_lst)
3961	IF ( ALLOCATED(f_rnk_lst) ) DEALLOCATE (f_rnk_lst)
3962
3963	IF ( coupling_mode == 'vnested_crse') THEN
3964	IF ( ALLOCATED (c2f_dims_cg) ) DEALLOCATE (c2f_dims_cg)
3965	IF ( ALLOCATED (f2c_dims_cg) ) DEALLOCATE (f2c_dims_cg)
3966	ELSEIF( coupling_mode == 'vnested_fine' ) THEN
3967	IF ( ALLOCATED (c2f_dims_fg) ) DEALLOCATE (c2f_dims_fg)
3968	IF ( ALLOCATED (f2c_dims_fg) ) DEALLOCATE (f2c_dims_fg)
3969	ENDIF
3970	#endif
3971	END SUBROUTINE vnest_deallocate
3972
3973	#endif
3974	END MODULE vertical_nesting_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |