Home

Context Navigation

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

Last change on this file since 2604 was 2516, checked in by suehring, 7 years ago
document changes
Property svn:keywords set to `Id`
File size: 139.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |