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

Last change on this file since 4121 was 4102, checked in by suehring, 5 years ago

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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