Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |