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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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