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

Last change on this file since 2712 was 2712, checked in by kanani, 4 years ago

Formatting and clean-up of vertical_nesting_mod

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