Home

Context Navigation

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

Last change on this file since 2712 was 2712, checked in by kanani, 6 years ago
Formatting and clean-up of vertical_nesting_mod
Property svn:keywords set to `Id`
File size: 153.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |