Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |