Home

Context Navigation

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

Last change on this file since 4143 was 4102, checked in by suehring, 6 years ago
Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.
Property svn:keywords set to `Id`
File size: 150.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |