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

Last change on this file since 3065 was 3065, checked in by Giersch, 3 years ago

New vertical stretching procedure has been introduced

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