Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |