Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |