Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |