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

Last change on this file since 4444 was 4444, checked in by raasch, 14 months ago

bugfix: cpp-directives for serial mode added

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