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

Last change on this file since 2516 was 2516, checked in by suehring, 4 years ago

document changes

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