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

Last change on this file since 4101 was 4101, checked in by gronemeier, 2 years ago

timestep.f90:

  • consider 2*Km within diffusion criterion as Km is considered twice within the diffusion of e,
  • in RANS mode, instead of considering each wind component individually use the wind speed of 3d wind vector in CFL criterion
  • do not limit the increase of dt based on its previous value in RANS mode

other:

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