Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |