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

Last change on this file since 3802 was 3802, checked in by raasch, 2 years ago

unused variables removed, unused subroutines commented out, type conversion added to avoid compiler warning about constant integer division truncation, script document_changes made bash compatible

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