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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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