Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |