Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |