Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |