Home

Context Navigation

vertical_nesting_mod.f90 @ 4101

Last change on this file since 4101 was 4101, checked in by gronemeier, 5 years ago

timestep.f90:

consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
in RANS mode, instead of considering each wind component individually use the wind speed of 3d wind vector in CFL criterion
do not limit the increase of dt based on its previous value in RANS mode

other:

remove dt_old

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: