Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |