Home

Context Navigation

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

Last change on this file since 2501 was 2374, checked in by kanani, 7 years ago
Added todo list for vertical_nesting_mod
Property svn:keywords set to `Id`
File size: 138.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |