Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |