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

Last change on this file since 3550 was 3467, checked in by suehring, 6 years ago

Branch salsa @3446 re-integrated into trunk

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