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

Last change on this file since 3186 was 3083, checked in by gronemeier, 6 years ago

merge with branch rans

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