Home

Context Navigation

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

Last change on this file since 4039 was 3802, checked in by raasch, 6 years ago
unused variables removed, unused subroutines commented out, type conversion added to avoid compiler warning about constant integer division truncation, script document_changes made bash compatible
Property svn:keywords set to `Id`
File size: 150.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |