Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |