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

Last change on this file since 4470 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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