Home

Context Navigation

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

Last change on this file since 4555 was 4481, checked in by maronga, 5 years ago
Bugfix for copyright updates in document_changes; copyright update applied to all files
Property svn:keywords set to `Id`
File size: 150.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |