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

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

Remove tabs from code, causing problems during merging

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