Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |