Home

Context Navigation

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

Last change on this file since 3760 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 150.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |