Home

Context Navigation

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

Last change on this file since 2941 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 153.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |