Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |