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

Last change on this file since 3917 was 3802, checked in by raasch, 6 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
Line 
1!> @file vertical_nesting_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Karlsruhe Institute of Technology
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: vertical_nesting_mod.f90 3802 2019-03-17 13:33:42Z raasch $
28! unused subroutines commented out
29!
30! 3655 2019-01-07 16:51:22Z knoop
31! unused variables removed
32!
33! 3232 2018-09-07 12:21:44Z raasch
34! references to mrun replaced by palmrun, and updated
35!
36! 3083 2018-06-19 14:03:12Z gronemeier
37! Error messages revised
38!
39! 3065 2018-06-12 07:03:02Z Giersch
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
44! Error messages revised
45!
46! 3045 2018-05-28 07:55:41Z Giersch
47! Error message revised
48!
49! 2718 2018-01-02 08:49:38Z maronga
50! Corrected "Former revisions" section
51!
52! 2712 2017-12-20 17:32:50Z kanani
53! Formatting and clean-up (SadiqHuq)
54!
55! 2696 2017-12-14 17:12:51Z kanani
56! Change in file header (GPL part)
57! Renamed diffusivities to tcm_diffusivities (TG)
58!
59! 2516 2017-10-04 11:03:04Z suehring
60! Remove tabs
61!
62! 2514 2017-10-04 09:52:37Z suehring
63! Remove tabs
64!
65! 2514 2017-10-04 09:52:37Z suehring
66! Added todo list
67!
68! 2365 2017-08-21 14:59:59Z kanani
69! Initial revision (SadiqHuq)
70!
71!
72!
73!
74! Description:
75! ------------
76!> Module for vertical nesting.
77!>
78!> Definition of parameters and variables for vertical nesting
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
90!>
91!> @todo Replace dz(1) appropriatly to account for grid stretching
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
100!------------------------------------------------------------------------------!
101 MODULE vertical_nesting_mod
102
103    USE kinds
104
105    IMPLICIT NONE
106
107    LOGICAL                                   ::  vnested = .FALSE.            !> set to true if palmrun
108                                                                               !> provides specific information via stdin
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
113
114
115
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
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
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
130
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
134
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
143
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
149
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
158
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
165
166
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
173    SAVE
174
175!-- Public functions
176    PUBLIC vnest_init_fine, vnest_boundary_conds, vnest_anterpolate,          &
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
180
181!-- Public constants and variables
182    PUBLIC vnested, vnest_init, vnest_start_time
183
184    PRIVATE bdims, bdims_rem,                                                 &
185            work3d, work2dusws, work2dvsws, work2dts, work2dus,               &
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,                                     &
192            ngp_c, ngp_f, target_idex, n_cell_c,                              &
193            offset, map_coord, TYPE_VNEST_BC, TYPE_VNEST_ANTER               
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
244#if defined( __parallel )
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
260        USE turbulence_closure_mod,                                            &
261            ONLY :  tcm_diffusivities
262       
263   
264        IMPLICIT NONE
265   
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
272   
273   
274        if (myid ==0 ) print *, ' TIME TO INIT FINE from COARSE', simulated_time
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
285            CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,       &
286                target_id, 0,                                                  &
287                terminate_coupled_remote, 1, MPI_INTEGER,                      &
288                target_id, 0,                                                  &
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
295            WRITE( message_string, * ) 'remote model "',                       &
296                TRIM( coupling_mode_remote ),                                  &
297                '" terminated',                                                &
298                '&with terminate_coupled_remote = ',                           &
299                terminate_coupled_remote,                                      &
300                '&local model  "', TRIM( coupling_mode ),                      &
301                '" has',                                                       &
302                '&terminate_coupled = ',                                       &
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,  &
361                        101, comm_inter, ierr)
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,  &
367                        102, comm_inter, ierr)
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,  &
373                        103, comm_inter, ierr)
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,  &
379                        105, comm_inter, ierr)
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,  &
386                        116, comm_inter, ierr)
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,  &
393                        104, comm_inter, ierr)
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,  &
399                        106, comm_inter, ierr)
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,  &
405                        107, comm_inter, ierr)
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   
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!   
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
481            call interpolate_to_fine_u
482   
483            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 102, &
484                comm_inter,status, ierr )
485            interpol3d => v
486            call interpolate_to_fine_v
487   
488            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 103, &
489                comm_inter,status, ierr )
490            interpol3d => w
491            call interpolate_to_fine_w
492   
493            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 105, &
494                comm_inter,status, ierr )
495            interpol3d => pt
496            call interpolate_to_fine_s
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
502            call interpolate_to_fine_s
503            ENDIF
504
505            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 104, &
506                comm_inter,status, ierr )
507            interpol3d => e
508            call interpolate_to_fine_s
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 )
513            call interpolate_to_fine_kh
514   
515            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 107, &
516                comm_inter,status, ierr )
517            call interpolate_to_fine_km
518   
519            DEALLOCATE(   work3d       )
520            NULLIFY   (   interpol3d   )
521   
522!-- Recv Surface Fluxes   
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   
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 )
555   
556            DEALLOCATE(   work2dusws   )
557            DEALLOCATE(   work2dvsws   )
558            DEALLOCATE(   work2dts     )
559            DEALLOCATE(   work2dus     )
560          ENDIF
561   
562          IF ( .NOT. constant_diffusion ) THEN
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)
566   
567                         IF ( e(kkf,jjf,iif) < 0.0 ) THEN
568                              e(kkf,jjf,iif) = 1E-15_wp
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
598!
599!-- Temperature boundary conditions at the bottom boundary
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
618               CALL tcm_diffusivities( pt, pt_reference )
619           ELSE
620               CALL tcm_diffusivities( vpt, pt_reference )
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
642
643    CONTAINS
644   
645       SUBROUTINE interpolate_to_fine_w
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     
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
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     
697                       DO iif = bottomx, topx
698     
699                           eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
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                     
705                           wprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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     
722                   DO iif = nxl, nxr
723                       DO jjf = bottomy, topy
724                     
725                           eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
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     
731                           wprs(k,jjf,iif) = eminus * wprf(k,j-1,iif) &
732                               + edot  * wprf(k,j,iif)   &
733                               + eplus  * wprf(k,j+1,iif)
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     
749               DO jjf = nys, nyn
750                   DO iif = nxl, nxr
751                       DO kkf = bottomz, topz
752     
753                           w(kkf,jjf,iif) = wprs(k,jjf,iif) + ( zwf(kkf) - zwc(k) ) &
754                               * ( wprs(k+1,jjf,iif) - wprs(k,jjf,iif) ) / dzc
755     
756                       END DO
757                   END DO
758               END DO
759     
760           END DO
761     
762           DO jjf = nys, nyn
763               DO iif = nxl, nxr
764     
765                   w(nzt,jjf,iif) = wprs(nztop,jjf,iif)
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     
776       SUBROUTINE interpolate_to_fine_u
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     
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
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
831                       DO jjf = bottomy, topy
832     
833                           eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
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     
839                           uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
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     
857               DO jjf = nys, nyn
858                   DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
859                       DO kkf = bottomz, topz
860                     
861                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
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     
867                           uprs(kkf,jjf,i) = eminus * uprf(k-1,jjf,i) &
868                               + edot  * uprf(k,jjf,i)   &
869                               + eplus  * uprf(k+1,jjf,i)
870     
871                       END DO
872                   END DO
873               END DO
874     
875           END DO
876     
877           DO jjf = nys, nyn
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     
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)
889     
890               END DO
891           END DO
892     
893           !
894           !-- Interpolation in x-direction (linear)
895     
896           DO kkf = nzb+1, nzt+1
897               DO jjf = nys, nyn
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     
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
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     
921       SUBROUTINE interpolate_to_fine_v
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
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
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     
974                       DO iif = bottomx, topx
975     
976                           eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
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     
982                           vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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
1001                   DO iif = nxl, nxr
1002                       DO kkf = bottomz, topz
1003     
1004                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
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     
1010                           vprs(kkf,j,iif) = eminus * vprf(k-1,j,iif) &
1011                               + edot  * vprf(k,j,iif)   &
1012                               + eplus  * vprf(k+1,j,iif)
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
1021               DO iif = nxl, nxr
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     
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)
1032     
1033               END DO
1034           END DO
1035     
1036           !
1037           !-- Interpolation in y-direction (linear)
1038     
1039           DO kkf = nzb+1, nzt+1
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     
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
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     
1064       SUBROUTINE interpolate_to_fine_s
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
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
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     
1118                       DO iif = bottomx, topx
1119     
1120                           eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
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     
1126                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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     
1144                   DO iif = nxl, nxr
1145                       DO jjf = bottomy, topy
1146     
1147                           eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
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     
1153                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1154                               + edot  * ptprf(k,j,iif)   &
1155                               + eplus  * ptprf(k,j+1,iif)
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     
1171               DO jjf = nys, nyn
1172                   DO iif = nxl, nxr
1173                       DO kkf = bottomz, topz
1174                       
1175                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
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       
1181                           interpol3d(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1182                               + edot  * ptprs(k,jjf,iif)   &
1183                               + eplus  * ptprs(k+1,jjf,iif)
1184     
1185                       END DO
1186                   END DO
1187               END DO
1188     
1189           END DO
1190               
1191           DO jjf = nys, nyn
1192               DO iif = nxl, nxr
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     
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)
1203     
1204               END DO
1205           END DO
1206     
1207     
1208           DEALLOCATE( ptprf, ptprs ) 
1209     
1210       END SUBROUTINE interpolate_to_fine_s
1211     
1212     
1213       SUBROUTINE interpolate_to_fine_kh
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
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
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     
1270                       DO iif = bottomx, topx
1271     
1272                           eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
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     
1278                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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     
1296                   DO iif = nxl, nxr
1297                       DO jjf = bottomy, topy
1298     
1299                           eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
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     
1305                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1306                               + edot  * ptprf(k,j,iif)   &
1307                               + eplus  * ptprf(k,j+1,iif)
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     
1323               DO jjf = nys, nyn
1324                   DO iif = nxl, nxr
1325                       DO kkf = bottomz, topz
1326                       
1327                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
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       
1333                           kh(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1334                               + edot  * ptprs(k,jjf,iif)   &
1335                               + eplus  * ptprs(k+1,jjf,iif)
1336     
1337                       END DO
1338                   END DO
1339               END DO
1340     
1341           END DO
1342               
1343           DO jjf = nys, nyn
1344               DO iif = nxl, nxr
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     
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)
1355     
1356               END DO
1357           END DO
1358     
1359     
1360           DEALLOCATE( ptprf, ptprs ) 
1361     
1362       END SUBROUTINE interpolate_to_fine_kh
1363     
1364       SUBROUTINE interpolate_to_fine_km
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
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
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     
1421                       DO iif = bottomx, topx
1422     
1423                           eps    = ( iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
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     
1429                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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     
1447                   DO iif = nxl, nxr
1448                       DO jjf = bottomy, topy
1449     
1450                           eps    = ( jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
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     
1456                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
1457                               + edot  * ptprf(k,j,iif)   &
1458                               + eplus  * ptprf(k,j+1,iif)
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     
1474               DO jjf = nys, nyn
1475                   DO iif = nxl, nxr
1476                       DO kkf = bottomz, topz
1477                       
1478                           eps    = ( zuf(kkf) - zuc(k) ) / dzc
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       
1484                           km(kkf,jjf,iif) = eminus * ptprs(k-1,jjf,iif) &
1485                               + edot  * ptprs(k,jjf,iif)   &
1486                               + eplus  * ptprs(k+1,jjf,iif)
1487     
1488                       END DO
1489                   END DO
1490               END DO
1491     
1492           END DO
1493               
1494           DO jjf = nys, nyn
1495               DO iif = nxl, nxr
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     
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)
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     
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!
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)
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
1642   
1643   
1644#endif       
1645    END SUBROUTINE vnest_init_fine
1646   
1647    SUBROUTINE vnest_boundary_conds
1648#if defined( __parallel )
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   
1668        INTEGER(iwp)                          :: i
1669        INTEGER(iwp)                          :: j
1670        INTEGER(iwp)                          :: iif
1671        INTEGER(iwp)                          :: jjf
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
1774               CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 209, &
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
1783            DO jjf = nys, nyn
1784                DO iif = nxl, nxr
1785                   e(nzt+1,jjf,iif) = e(nzt,jjf,iif)
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
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
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       
1855                   DO iif = bottomx, topx
1856       
1857                       eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
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
1862                       wprf(j,iif) = eminus * work3d(bdims_rem(3,1),j,i-1) &
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       
1880               DO iif = nxl, nxr
1881       
1882                   DO jjf = bottomy, topy
1883       
1884                       eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
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       
1894                       w(nzt,jjf,iif) = eminus * wprf(j-1,iif) &
1895                           + edot  * wprf(j,iif)   &
1896                           + eplus  * wprf(j+1,iif)
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
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
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
1952                       DO jjf = bottomy, topy
1953       
1954                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
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       
1960                           uprf(k,jjf,i) = eminus * work3d(k,j-1,i) &
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       
1972           DO jjf = nys, nyn
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
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)
1982               END DO
1983           END DO
1984       
1985           !
1986           !-- Interpolation in x-direction
1987       
1988           DO jjf = nys, nyn
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       
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
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
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
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       
2057                       DO iif = bottomx, topx
2058       
2059                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
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
2064                           vprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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
2077               DO iif = nxl, nxr
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
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)
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)
2095               DO iif = nxl, nxr
2096       
2097                   bottomy = (nyf+1)/(nyc+1) * j
2098                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
2099       
2100                   DO jjf = bottomy, topy
2101       
2102                       v(nzt+1,jjf,iif) = vprs(j,iif) + ( jjf * dyf - j * dyc ) * ( vprs(j+1,iif) - vprs(j,iif) ) / dyc
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
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
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       
2166                       DO iif = bottomx, topx
2167       
2168                           eps = (iif * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
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       
2178                           ptprf(k,j,iif) = eminus * work3d(k,j,i-1) &
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       
2199                   DO iif = nxl, nxr
2200       
2201                       DO jjf = bottomy, topy
2202       
2203                           eps = (jjf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
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       
2213                           ptprs(k,jjf,iif) = eminus * ptprf(k,j-1,iif) &
2214                               + edot  * ptprf(k,j,iif)   &
2215                               + eplus  * ptprf(k,j+1,iif)
2216                       END DO
2217       
2218                   END DO
2219       
2220               END DO
2221       
2222           END DO
2223       
2224           !
2225           !-- Interpolation in z-direction
2226       
2227           DO jjf = nys, nyn
2228               DO iif = nxl, nxr
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       
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)
2243       
2244               END DO
2245           END DO         
2246           
2247           DEALLOCATE ( ptprf, ptprs )
2248       
2249       
2250       
2251       END SUBROUTINE vnest_set_topbc_s
2252#endif
2253    END SUBROUTINE vnest_boundary_conds
2254   
2255 
2256    SUBROUTINE vnest_boundary_conds_khkm
2257#if defined( __parallel )
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   
2278        INTEGER(iwp)                          :: i
2279        INTEGER(iwp)                          :: j
2280        INTEGER(iwp)                          :: iif
2281        INTEGER(iwp)                          :: jjf
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,    &
2319                        207, comm_inter, ierr)
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,    &
2325                        208, comm_inter, ierr)
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
2367        DO jjf = nys, nyn
2368            DO iif = nxl, nxr
2369               kh(nzt+1,jjf,iif) = kh(nzt,jjf,iif)
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
2377        DO jjf = nys, nyn
2378            DO iif = nxl, nxr
2379               km(nzt+1,jjf,iif) = km(nzt,jjf,iif)
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
2386        !!    DO jjf = nys-1, nyn+1
2387        !!       DO iif = nxl-1, nxr+1
2388        !!
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)
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   
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
2541       
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
2678   
2679 
2680#endif
2681    END SUBROUTINE vnest_boundary_conds_khkm
2682   
2683   
2684   
2685    SUBROUTINE vnest_anterpolate
2686
2687#if defined( __parallel )
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   
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
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
2736            WRITE( message_string, * ) 'remote model "',                       &
2737                TRIM( coupling_mode_remote ),                                  &
2738                '" terminated',                                                &
2739                '&with terminate_coupled_remote = ',                           &
2740                terminate_coupled_remote,                                      &
2741                '&local model  "', TRIM( coupling_mode ),                      &
2742                '" has',                                                       &
2743                '&terminate_coupled = ',                                       &
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   
2817                 IF ( humidity ) THEN
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 )
2824                 ENDIF
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
2910     
2911            CALL anterpolate_to_crse_u
2912            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
2913                101, comm_inter, ierr)
2914   
2915            anterpol3d => v
2916   
2917            CALL anterpolate_to_crse_v
2918            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
2919                102, comm_inter, ierr)
2920   
2921            anterpol3d => pt
2922   
2923            CALL anterpolate_to_crse_s
2924            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
2925                105, comm_inter, ierr)
2926   
2927   
2928          IF ( humidity ) THEN
2929   
2930            anterpol3d => q
2931   
2932            CALL anterpolate_to_crse_s
2933            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
2934                106, comm_inter, ierr)
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
2943            CALL anterpolate_to_crse_w
2944            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
2945                103, comm_inter, ierr)
2946   
2947            NULLIFY   (   anterpol3d                )
2948            DEALLOCATE(   work3d                    )
2949   
2950        ENDIF
2951   
2952   
2953
2954    CONTAINS
2955       SUBROUTINE anterpolate_to_crse_u
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
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
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     
2996                       iif = (nxf+1) / (nxc+1) * i
2997     
2998                       aweight   = 0.0
2999     
3000                       DO kkf = bottomz, topz
3001                           DO jjf = bottomy, topy
3002     
3003                               aweight   = aweight + anterpol3d(kkf,jjf,iif) *        &
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     
3022       SUBROUTINE anterpolate_to_crse_v
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
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
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     
3058                   jjf = (nyf+1) / (nyc+1) * j
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     
3067                       DO kkf = bottomz, topz
3068                           DO iif = bottomx, topx
3069     
3070                               aweight   = aweight + anterpol3d(kkf,jjf,iif) *        &
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     
3088       SUBROUTINE anterpolate_to_crse_w
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
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
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     
3119               kkf = cfratio(3) * k
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     
3133                       DO jjf = bottomy, topy
3134                           DO iif = bottomx, topx
3135     
3136                               aweight   = aweight + anterpol3d (kkf,jjf,iif) *        &
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     
3154       SUBROUTINE anterpolate_to_crse_s
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
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
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         
3201                       DO kkf = bottomz, topz
3202                           DO jjf = bottomy, topy
3203                               DO iif = bottomx, topx
3204     
3205                                   aweight   = aweight  + anterpol3d(kkf,jjf,iif) *                  &
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
3222#endif       
3223    END SUBROUTINE vnest_anterpolate
3224   
3225   
3226   
3227    SUBROUTINE vnest_anterpolate_e
3228#if defined( __parallel )
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   
3246        REAL(wp)                              :: time_since_reference_point_rem
3247        INTEGER(iwp)                          :: i
3248        INTEGER(iwp)                          :: j
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
3259            CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,       &
3260                target_id, 0,                                                  &
3261                terminate_coupled_remote, 1, MPI_INTEGER,                      &
3262                target_id, 0,                                                  &
3263                comm_inter, status, ierr )
3264        ENDIF
3265        CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d,   &
3266            ierr )
3267   
3268        IF ( terminate_coupled_remote > 0 )  THEN
3269            WRITE( message_string, * ) 'remote model "',                       &
3270                TRIM( coupling_mode_remote ),                                  &
3271                '" terminated',                                                &
3272                '&with terminate_coupled_remote = ',                           &
3273                terminate_coupled_remote,                                      &
3274                '&local model  "', TRIM( coupling_mode ),                      & 
3275                '" has',                                                       &
3276                '&terminate_coupled = ',                                       &
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 (3,2) = f2c_dims_cg (5,map_coord(1),map_coord(2)) 
3316
3317
3318                     n_cell_c = (bdims(1,2)-bdims(1,1)+1) * &
3319                         (bdims(2,2)-bdims(2,1)+1) * &
3320                         (bdims(3,2)-bdims(3,1)+0)
3321
3322
3323                    CALL MPI_RECV( e( bdims(3,1)+1:bdims(3,2), &
3324                        bdims(2,1)  :bdims(2,2), &
3325                        bdims(1,1)  :bdims(1,2)),&
3326                        n_cell_c, MPI_REAL, target_idex, 104, &
3327                        comm_inter,status, ierr )
3328                end do
3329            end do
3330   
3331   
3332            !
3333            !-- Boundary conditions
3334   
3335            IF ( .NOT. constant_diffusion ) THEN
3336                e(nzb,:,:)   = e(nzb+1,:,:)
3337            END IF
3338   
3339            IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
3340   
3341        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
3342                      ! Send data to coarse grid for anterpolation
3343   
3344            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
3345            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
3346            map_coord(1) = offset(1)
3347            map_coord(2) = offset(2)
3348            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
3349
3350            bdims_rem (1,1) =  f2c_dims_fg (0)
3351            bdims_rem (1,2) =  f2c_dims_fg (1)
3352            bdims_rem (2,1) =  f2c_dims_fg (2)
3353            bdims_rem (2,2) =  f2c_dims_fg (3)
3354            bdims_rem (3,1) =  f2c_dims_fg (4)
3355            bdims_rem (3,2) =  f2c_dims_fg (5)
3356
3357            ALLOCATE( work3d (                   &
3358                bdims_rem(3,1)+1:bdims_rem(3,2), &
3359                bdims_rem(2,1)  :bdims_rem(2,2), &
3360                bdims_rem(1,1)  :bdims_rem(1,2)))
3361   
3362            anterpol3d => e
3363   
3364            CALL anterpolate_to_crse_e
3365
3366            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
3367                           104, comm_inter, ierr)
3368
3369            NULLIFY   (   anterpol3d                )
3370            DEALLOCATE(   work3d                    )
3371        ENDIF
3372   
3373   
3374    CONTAINS
3375   
3376   
3377   
3378   
3379   
3380       SUBROUTINE anterpolate_to_crse_e
3381     
3382     
3383           USE arrays_3d
3384           USE control_parameters
3385           USE grid_variables
3386           USE indices
3387           USE pegrid
3388           
3389     
3390           IMPLICIT NONE
3391
3392           INTEGER(iwp)                       :: i
3393           INTEGER(iwp)                       :: j
3394           INTEGER(iwp)                       :: k
3395           INTEGER(iwp)                       :: iif
3396           INTEGER(iwp)                       :: jjf
3397           INTEGER(iwp)                       :: kkf
3398           INTEGER(iwp)                       :: bottomx
3399           INTEGER(iwp)                       :: bottomy
3400           INTEGER(iwp)                       :: bottomz
3401           INTEGER(iwp)                       :: topx
3402           INTEGER(iwp)                       :: topy
3403           INTEGER(iwp)                       :: topz
3404           REAL(wp)                           :: aweight_a
3405           REAL(wp)                           :: aweight_b
3406           REAL(wp)                           :: aweight_c
3407           REAL(wp)                           :: aweight_d
3408           REAL(wp)                           :: aweight_e
3409           REAL(wp)                           :: energ 
3410         
3411           DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
3412     
3413               bottomz = (dzc/dzf) * (k-1) + 1
3414               topz    = (dzc/dzf) * k
3415             
3416               DO j = bdims_rem(2,1), bdims_rem(2,2)
3417     
3418                   bottomy = (nyf+1) / (nyc+1) * j
3419                   topy    = (nyf+1) / (nyc+1) * (j+1) - 1
3420     
3421                   DO i = bdims_rem(1,1), bdims_rem(1,2)
3422     
3423                       bottomx = (nxf+1) / (nxc+1) * i
3424                       topx    = (nxf+1) / (nxc+1) * (i+1) - 1
3425     
3426                       aweight_a   = 0.0
3427                       aweight_b   = 0.0
3428                       aweight_c   = 0.0
3429                       aweight_d   = 0.0
3430                       aweight_e   = 0.0
3431     
3432                       DO kkf = bottomz, topz
3433                           DO jjf = bottomy, topy
3434                               DO iif = bottomx, topx
3435     
3436                                   aweight_a = aweight_a + anterpol3d(kkf,jjf,iif)  *                     &
3437                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3438                             
3439     
3440                                   energ = ( 0.5 * ( u(kkf,jjf,iif)   + u(kkf,jjf,iif+1) ) )**2.0 +          &
3441                                       ( 0.5 * ( v(kkf,jjf,iif)   + v(kkf,jjf+1,iif) ) )**2.0 +              &
3442                                       ( 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) )**2.0
3443     
3444                                   aweight_b   =   aweight_b + energ         *                         &
3445                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3446       
3447                                   aweight_c   =   aweight_c + 0.5 * ( u(kkf,jjf,iif) + u(kkf,jjf,iif+1) ) * &
3448                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3449     
3450                                   aweight_d   =   aweight_d + 0.5 * ( v(kkf,jjf,iif) + v(kkf,jjf+1,iif) ) * &
3451                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3452       
3453                                   aweight_e   =   aweight_e + 0.5 * ( w(kkf-1,jjf,iif) + w(kkf,jjf,iif) ) * &
3454                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
3455     
3456     
3457                               END DO
3458                           END DO
3459                       END DO
3460     
3461                       work3d(k,j,i)  = aweight_a + 0.5 * ( aweight_b -  &
3462                           aweight_c**2.0 -                              &
3463                           aweight_d**2.0 -                              &
3464                           aweight_e**2.0                 )
3465     
3466                   END DO
3467     
3468               END DO
3469     
3470           END DO
3471     
3472       
3473     
3474       END SUBROUTINE anterpolate_to_crse_e
3475#endif       
3476    END SUBROUTINE vnest_anterpolate_e
3477
3478    SUBROUTINE vnest_init_pegrid_rank
3479#if defined( __parallel )
3480! Domain decomposition and exchange of grid variables between coarse and fine
3481! Given processor coordinates as index f_rnk_lst(pcoord(1), pcoord(2))
3482!   returns the rank. A single coarse block will have to send data to multiple
3483!   fine blocks. In the coarse grid the pcoords of the remote block is first found and then using
3484!   f_rnk_lst the target_idex is identified.
3485! blk_dim stores the index limits of a given block. blk_dim_remote is received
3486! from the asscoiated nest partner.
3487! cf_ratio(1:3) is the ratio between fine and coarse grid: nxc/nxf, nyc/nyf and
3488! ceiling(dxc/dxf)
3489
3490
3491       USE control_parameters,                                                 &
3492           ONLY:  coupling_mode
3493
3494       USE kinds
3495         
3496       USE pegrid
3497 
3498
3499       IMPLICIT NONE
3500
3501       INTEGER(iwp)                           :: dest_rnk
3502       INTEGER(iwp)                           ::  i                        !<
3503
3504       IF (myid == 0) THEN
3505           IF ( coupling_mode == 'vnested_crse') THEN
3506               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 33, comm_inter, &
3507                   ierr )
3508               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, numprocs, 66,     &
3509                   comm_inter, status, ierr )
3510           ELSEIF ( coupling_mode == 'vnested_fine') THEN
3511               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 33,      &
3512                   comm_inter, status, ierr )
3513               CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 66, comm_inter, &
3514                   ierr )
3515           ENDIF
3516       ENDIF
3517
3518
3519       IF ( coupling_mode == 'vnested_crse') THEN
3520           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3521           ALLOCATE( c_rnk_lst( 0:(pdims(1)-1)    ,0:(pdims(2)-1) )   )
3522           ALLOCATE( f_rnk_lst( 0:(pdims_partner(1)-1)  ,0:(pdims_partner(2)-1) ) )
3523           do i=0,numprocs-1
3524               CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
3525               call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
3526               c_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
3527           end do
3528       ELSEIF ( coupling_mode == 'vnested_fine') THEN
3529           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3530           ALLOCATE( c_rnk_lst( 0:(pdims_partner(1)-1)  ,0:(pdims_partner(2)-1) ) )
3531           ALLOCATE( f_rnk_lst( 0:(pdims(1)-1)    ,0:(pdims(2)-1) )   )
3532
3533           do i=0,numprocs-1
3534               CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
3535               call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
3536               f_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
3537           enddo
3538       ENDIF
3539
3540
3541       IF ( coupling_mode == 'vnested_crse') THEN
3542           if (myid == 0) then
3543               CALL MPI_SEND( c_rnk_lst, pdims(1)*pdims(2),     MPI_INTEGER, numprocs, 0, comm_inter, ierr )
3544               CALL MPI_RECV( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, numprocs, 4, comm_inter,status, ierr )
3545           end if
3546           CALL MPI_BCAST( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
3547       ELSEIF ( coupling_mode == 'vnested_fine') THEN
3548           if (myid == 0) then
3549               CALL MPI_RECV( c_rnk_lst, pdims_partner(1)*pdims_partner(2),  MPI_INTEGER, 0, 0, comm_inter,status, ierr )
3550               CALL MPI_SEND( f_rnk_lst, pdims(1)*pdims(2),      MPI_INTEGER, 0, 4, comm_inter, ierr )
3551           end if
3552           CALL MPI_BCAST( c_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
3553       ENDIF
3554
3555       !-- Reason for MPI error unknown; solved if three lines duplicated
3556       CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
3557       CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
3558       CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
3559
3560
3561#endif
3562 
3563    END SUBROUTINE vnest_init_pegrid_rank
3564
3565
3566    SUBROUTINE vnest_init_pegrid_domain
3567#if defined( __parallel )
3568
3569       USE control_parameters,                                                 &
3570           ONLY:  coupling_mode, coupling_topology, dz,                        &
3571                  dz_stretch_level_start, message_string
3572
3573       USE grid_variables,                                                     &
3574           ONLY:  dx, dy
3575           
3576       USE indices,                                                            &
3577           ONLY: nbgp, nx, ny, nz, nxl, nxr, nys, nyn, nzb, nzt
3578
3579       USE kinds
3580         
3581       USE pegrid
3582       
3583       IMPLICIT NONE
3584
3585       INTEGER(iwp)                           :: i              !<
3586       INTEGER(iwp)                           :: j              !<
3587       INTEGER(iwp)                           :: tempx
3588       INTEGER(iwp)                           :: tempy
3589       INTEGER(iwp)                           :: TYPE_INT_YZ
3590       INTEGER(iwp)                           :: SIZEOFREAL
3591       INTEGER(iwp)                           :: MTV_X
3592       INTEGER(iwp)                           :: MTV_Y
3593       INTEGER(iwp)                           :: MTV_Z
3594       INTEGER(iwp)                           :: MTV_RX
3595       INTEGER(iwp)                           :: MTV_RY
3596       INTEGER(iwp)                           :: MTV_RZ
3597   
3598       !
3599       !--    Pass the number of grid points of the coarse model to
3600       !--    the nested model and vice versa
3601       IF ( coupling_mode ==  'vnested_crse' )  THEN
3602
3603           nxc = nx
3604           nyc = ny
3605           nzc = nz
3606           dxc = dx
3607           dyc = dy
3608           dzc = dz(1)
3609           cg_nprocs = numprocs
3610
3611           IF ( myid == 0 )  THEN
3612
3613               CALL MPI_SEND( nxc, 1, MPI_INTEGER  , numprocs, 1,  comm_inter, &
3614                   ierr )
3615               CALL MPI_SEND( nyc, 1, MPI_INTEGER  , numprocs, 2,  comm_inter, &
3616                   ierr )
3617               CALL MPI_SEND( nzc, 1, MPI_INTEGER  , numprocs, 3,  comm_inter, &
3618                   ierr )
3619               CALL MPI_SEND( dxc, 1, MPI_REAL     , numprocs, 4,  comm_inter, &
3620                   ierr )
3621               CALL MPI_SEND( dyc, 1, MPI_REAL     , numprocs, 5,  comm_inter, &
3622                   ierr )
3623               CALL MPI_SEND( dzc, 1, MPI_REAL     , numprocs, 6,  comm_inter, &
3624                   ierr )
3625               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7,  comm_inter, &
3626                   ierr )
3627               CALL MPI_SEND( cg_nprocs, 1, MPI_INTEGER, numprocs, 8,  comm_inter,  &
3628                   ierr )
3629               CALL MPI_RECV( nxf, 1, MPI_INTEGER,   numprocs, 21, comm_inter, &
3630                   status, ierr )
3631               CALL MPI_RECV( nyf, 1, MPI_INTEGER,   numprocs, 22, comm_inter, &
3632                   status, ierr )
3633               CALL MPI_RECV( nzf, 1, MPI_INTEGER,   numprocs, 23, comm_inter, &
3634                   status, ierr )
3635               CALL MPI_RECV( dxf, 1, MPI_REAL,      numprocs, 24, comm_inter, &
3636                   status, ierr )
3637               CALL MPI_RECV( dyf, 1, MPI_REAL,      numprocs, 25, comm_inter, &
3638                   status, ierr )
3639               CALL MPI_RECV( dzf, 1, MPI_REAL,      numprocs, 26, comm_inter, &
3640                   status, ierr )
3641               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER,                   &
3642                   numprocs, 27, comm_inter, status, ierr )
3643               CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER,                       &
3644                   numprocs, 28, comm_inter, status, ierr )
3645           ENDIF
3646
3647           CALL MPI_BCAST( nxf, 1,   MPI_INTEGER, 0, comm2d, ierr )
3648           CALL MPI_BCAST( nyf, 1,   MPI_INTEGER, 0, comm2d, ierr )
3649           CALL MPI_BCAST( nzf, 1,   MPI_INTEGER, 0, comm2d, ierr )
3650           CALL MPI_BCAST( dxf, 1,   MPI_REAL,    0, comm2d, ierr )
3651           CALL MPI_BCAST( dyf, 1,   MPI_REAL,    0, comm2d, ierr )
3652           CALL MPI_BCAST( dzf, 1,   MPI_REAL,    0, comm2d, ierr )
3653           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
3654           CALL MPI_BCAST( fg_nprocs,  1, MPI_INTEGER, 0, comm2d, ierr )
3655           
3656!
3657!--        Check if stretching is used within the nested domain. ABS(...) is
3658!--        necessary because of the default value of -9999999.9_wp (negative)
3659           IF ( ABS( dz_stretch_level_start(1) ) <= (nzf+1)*dzf )  THEN       
3660               message_string = 'Stretching in the parent domain is '//        &
3661                                'only allowed above the nested domain'
3662               CALL message( 'vnest_init_pegrid_domain', 'PA0497', 1, 2, 0, 6, 0 )
3663           ENDIF
3664
3665       ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
3666
3667           nxf = nx
3668           nyf = ny
3669           nzf = nz
3670           dxf = dx
3671           dyf = dy
3672           dzf = dz(1)
3673           fg_nprocs = numprocs
3674
3675           IF ( myid == 0 ) THEN
3676
3677               CALL MPI_RECV( nxc, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
3678                   ierr )
3679               CALL MPI_RECV( nyc, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
3680                   ierr )
3681               CALL MPI_RECV( nzc, 1, MPI_INTEGER, 0, 3, comm_inter, status, &
3682                   ierr )
3683               CALL MPI_RECV( dxc, 1, MPI_REAL,    0, 4, comm_inter, status, &
3684                   ierr )
3685               CALL MPI_RECV( dyc, 1, MPI_REAL,    0, 5, comm_inter, status, &
3686                   ierr )
3687               CALL MPI_RECV( dzc, 1, MPI_REAL,    0, 6, comm_inter, status, &
3688                   ierr )
3689               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 7, comm_inter, &
3690                   status, ierr )
3691               CALL MPI_RECV( cg_nprocs,  1, MPI_INTEGER, 0, 8, comm_inter, &
3692                   status, ierr )
3693               CALL MPI_SEND( nxf, 1, MPI_INTEGER, 0, 21,  comm_inter, ierr )
3694               CALL MPI_SEND( nyf, 1, MPI_INTEGER, 0, 22,  comm_inter, ierr )
3695               CALL MPI_SEND( nzf, 1, MPI_INTEGER, 0, 23, comm_inter, ierr )
3696               CALL MPI_SEND( dxf, 1, MPI_REAL,    0, 24, comm_inter, ierr )
3697               CALL MPI_SEND( dyf, 1, MPI_REAL,    0, 25, comm_inter, ierr )
3698               CALL MPI_SEND( dzf, 1, MPI_REAL,    0, 26, comm_inter, ierr )
3699               CALL MPI_SEND( pdims,2,MPI_INTEGER, 0, 27, comm_inter, ierr )
3700               CALL MPI_SEND( fg_nprocs,1,MPI_INTEGER, 0, 28, comm_inter, ierr )
3701           ENDIF
3702
3703           CALL MPI_BCAST( nxc, 1, MPI_INTEGER, 0, comm2d, ierr)
3704           CALL MPI_BCAST( nyc, 1, MPI_INTEGER, 0, comm2d, ierr)
3705           CALL MPI_BCAST( nzc, 1, MPI_INTEGER, 0, comm2d, ierr)
3706           CALL MPI_BCAST( dxc, 1, MPI_REAL,    0, comm2d, ierr)
3707           CALL MPI_BCAST( dyc, 1, MPI_REAL,    0, comm2d, ierr)
3708           CALL MPI_BCAST( dzc, 1, MPI_REAL,    0, comm2d, ierr)
3709           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr)
3710           CALL MPI_BCAST( cg_nprocs,  1, MPI_INTEGER, 0, comm2d, ierr)
3711
3712       ENDIF
3713       
3714       ngp_c = ( nxc+1 + 2 * nbgp ) * ( nyc+1 + 2 * nbgp )
3715       ngp_f = ( nxf+1 + 2 * nbgp ) * ( nyf+1 + 2 * nbgp )
3716
3717       IF ( coupling_mode(1:8) == 'vnested_')  coupling_topology = 1
3718
3719
3720!-- Nesting Ratio: For each coarse grid cell how many fine grid cells exist
3721       cfratio(1) = INT ( (nxf+1) / (nxc+1) )
3722       cfratio(2) = INT ( (nyf+1) / (nyc+1) )
3723       cfratio(3) = CEILING (  dzc    /  dzf  )
3724
3725!-- target_id is used only for exhange of information like simulated_time
3726!-- which are then MPI_BCAST to other processors in the group
3727          IF ( myid == 0 )  THEN
3728   
3729              IF ( TRIM( coupling_mode ) == 'vnested_crse' )  THEN
3730                  target_id = numprocs
3731              ELSE IF ( TRIM( coupling_mode ) == 'vnested_fine' )  THEN
3732                  target_id = 0
3733              ENDIF
3734   
3735          ENDIF
3736
3737!-- Store partner grid dimenstions and create MPI derived types
3738
3739         IF ( coupling_mode == 'vnested_crse' )  THEN
3740
3741            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
3742            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
3743     
3744            tempx =  ( pdims_partner(1) / pdims(1) ) - 1
3745            tempy =  ( pdims_partner(2) / pdims(2) ) - 1
3746            ALLOCATE( c2f_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
3747            ALLOCATE( f2c_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
3748
3749            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
3750                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
3751                   map_coord(1) = i+offset(1)
3752                   map_coord(2) = j+offset(2)
3753     
3754                   target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
3755
3756                   CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 10, &
3757                       comm_inter,status, ierr )
3758
3759!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
3760!-- One CG can have multiple FG partners. The 3D array is mapped by partner proc co-ord
3761                   c2f_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
3762                   c2f_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
3763                   c2f_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
3764                   c2f_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
3765                   c2f_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,2) / cfratio(3)
3766                   c2f_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2) / cfratio(3)) + 2
3767 
3768!-- Store the CG dimensions that  correspond to the FG partner; needed for anterpolation
3769                   f2c_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
3770                   f2c_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
3771                   f2c_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
3772                   f2c_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
3773                   f2c_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,1)
3774                   f2c_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2)-cfratio(3))/ cfratio(3)
3775   
3776                   CALL MPI_SEND( c2f_dims_cg (:,map_coord(1),map_coord(2)), 6, &
3777                      MPI_INTEGER, target_idex, 100, comm_inter, ierr )
3778     
3779                   CALL MPI_SEND( f2c_dims_cg (:,map_coord(1),map_coord(2)), 6, &
3780                      MPI_INTEGER, target_idex, 101, comm_inter, ierr )
3781     
3782                end do
3783            end do
3784     
3785!-- A derived data type to pack 3 Z-levels of CG to set FG top BC
3786            MTV_X = ( nxr - nxl + 1 ) + 2*nbgp
3787            MTV_Y = ( nyn - nys + 1 ) + 2*nbgp
3788            MTV_Z = nzt+1 - nzb +1
3789           
3790            MTV_RX = ( c2f_dims_cg (1,offset(1),offset(2)) - c2f_dims_cg (0,offset(1),offset(2)) ) +1+2
3791            MTV_RY = ( c2f_dims_cg (3,offset(1),offset(2)) - c2f_dims_cg (2,offset(1),offset(2)) ) +1+2
3792            MTV_RZ = ( c2f_dims_cg (5,offset(1),offset(2)) - c2f_dims_cg (4,offset(1),offset(2)) ) +1
3793           
3794            CALL MPI_TYPE_EXTENT(MPI_REAL, SIZEOFREAL, IERR)
3795
3796            CALL MPI_TYPE_VECTOR ( MTV_RY, MTV_RZ, MTV_Z, MPI_REAL, TYPE_INT_YZ, IERR)
3797            CALL MPI_TYPE_HVECTOR( MTV_RX, 1, MTV_Z*MTV_Y*SIZEOFREAL, &
3798                                              TYPE_INT_YZ, TYPE_VNEST_BC, IERR)
3799            CALL MPI_TYPE_FREE(TYPE_INT_YZ, IERR)
3800            CALL MPI_TYPE_COMMIT(TYPE_VNEST_BC, IERR)   
3801
3802
3803         ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
3804
3805           ALLOCATE( c2f_dims_fg (0:5) )
3806           ALLOCATE( f2c_dims_fg (0:5) )
3807   
3808           offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
3809           offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
3810           map_coord(1) = offset(1)
3811           map_coord(2) = offset(2)
3812           target_idex = c_rnk_lst(map_coord(1),map_coord(2))
3813   
3814           bdims (1,1) = nxl
3815           bdims (1,2) = nxr
3816           bdims (2,1) = nys
3817           bdims (2,2) = nyn
3818           bdims (3,1) = nzb
3819           bdims (3,2) = nzt
3820
3821           CALL MPI_SEND( bdims,    6, MPI_INTEGER, target_idex, 10, &
3822               comm_inter, ierr )
3823           
3824!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
3825!-- One FG can have only one CG partner
3826           CALL MPI_RECV( c2f_dims_fg,   6, MPI_INTEGER, target_idex, 100, &
3827               comm_inter,status, ierr )
3828           
3829           CALL MPI_RECV( f2c_dims_fg,   6, MPI_INTEGER, target_idex, 101, &
3830               comm_inter,status, ierr )
3831 
3832!-- Store the CG dimensions that  correspond to the FG partner; needed for anterpolation
3833
3834           n_cell_c = (f2c_dims_fg(1)-f2c_dims_fg(0)+1) * &
3835                      (f2c_dims_fg(3)-f2c_dims_fg(2)+1) * &
3836                      (f2c_dims_fg(5)-f2c_dims_fg(4)+0)
3837
3838           CALL MPI_TYPE_CONTIGUOUS(n_cell_c, MPI_REAL, TYPE_VNEST_ANTER, IERR)
3839           CALL MPI_TYPE_COMMIT(TYPE_VNEST_ANTER, ierr) 
3840
3841        ENDIF
3842#endif   
3843       END SUBROUTINE vnest_init_pegrid_domain
3844
3845
3846       SUBROUTINE vnest_init_grid
3847 
3848#if defined( __parallel )
3849          USE arrays_3d,                                                       &
3850              ONLY:  zu, zw
3851             
3852          USE control_parameters,                                              &
3853              ONLY:  coupling_mode, message_string, number_stretch_level_start
3854             
3855          USE indices,                                                         &
3856              ONLY:  nzt
3857         
3858          USE kinds
3859         
3860          USE pegrid
3861
3862          IMPLICIT NONE
3863     
3864!
3865!--       Allocate and Exchange zuc and zuf, zwc and zwf
3866          IF ( coupling_mode(1:8)  == 'vnested_' )  THEN
3867         
3868             ALLOCATE( zuc(0:nzc+1), zuf(0:nzf+1) )
3869             ALLOCATE( zwc(0:nzc+1), zwf(0:nzf+1) )
3870         
3871             IF ( coupling_mode ==  'vnested_crse' )  THEN
3872               
3873                zuc = zu
3874                zwc = zw
3875               
3876                IF ( myid == 0 )  THEN
3877         
3878                   CALL MPI_SEND( zuc, nzt+2, MPI_REAL, numprocs, 41, comm_inter,  &
3879                                  ierr )
3880                   CALL MPI_RECV( zuf, nzf+2, MPI_REAL, numprocs, 42, comm_inter,  &
3881                                  status, ierr )
3882         
3883                   CALL MPI_SEND( zwc, nzt+2, MPI_REAL, numprocs, 43, comm_inter,  &
3884                                  ierr )
3885                   CALL MPI_RECV( zwf, nzf+2, MPI_REAL, numprocs, 44, comm_inter,  &
3886                                  status, ierr )
3887         
3888                ENDIF
3889         
3890                CALL MPI_BCAST( zuf,nzf+2,MPI_REAL,    0, comm2d, ierr ) 
3891                CALL MPI_BCAST( zwf,nzf+2,MPI_REAL,    0, comm2d, ierr ) 
3892         
3893             ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
3894         
3895!
3896!--             Check if stretching is used within the nested domain
3897                IF ( number_stretch_level_start > 0 )  THEN
3898                   message_string = 'Stretching in the nested domain is not '//&
3899                           'allowed'
3900                   CALL message( 'vnest_init_grid', 'PA0498', 1, 2, 0, 6, 0 )
3901                ENDIF
3902               
3903                zuf = zu
3904                zwf = zw
3905               
3906                IF ( myid == 0 )  THEN
3907         
3908                   CALL MPI_RECV( zuc,nzc+2, MPI_REAL, 0, 41, comm_inter, status, &
3909                                  ierr )
3910                   CALL MPI_SEND( zuf,nzt+2, MPI_REAL, 0, 42, comm_inter, ierr )
3911         
3912                   CALL MPI_RECV( zwc,nzc+2, MPI_REAL, 0, 43, comm_inter, status, &
3913                                  ierr )
3914                   CALL MPI_SEND( zwf,nzt+2, MPI_REAL, 0, 44, comm_inter, ierr )
3915                ENDIF
3916         
3917                CALL MPI_BCAST( zuc,nzc+2,MPI_REAL,  0, comm2d, ierr ) 
3918                CALL MPI_BCAST( zwc,nzc+2,MPI_REAL,  0, comm2d, ierr ) 
3919         
3920             ENDIF
3921          ENDIF
3922
3923#endif
3924       END SUBROUTINE vnest_init_grid
3925
3926
3927       SUBROUTINE vnest_check_parameters
3928#if defined( __parallel )
3929
3930          USE pegrid,                                                                &
3931              ONLY:  myid
3932
3933          IMPLICIT NONE
3934
3935          IF (myid==0) PRINT*, '*** vnest: check parameters not implemented yet ***'
3936     
3937#endif
3938       END SUBROUTINE vnest_check_parameters
3939
3940
3941       SUBROUTINE vnest_timestep_sync
3942
3943#if defined( __parallel )
3944         USE control_parameters,                                                    &
3945             ONLY:  coupling_mode, dt_3d, old_dt
3946     
3947         USE interfaces
3948     
3949         USE kinds
3950
3951         USE pegrid
3952
3953         IMPLICIT NONE
3954
3955         IF ( coupling_mode == 'vnested_crse')  THEN
3956            dtc = dt_3d
3957               if (myid == 0) then
3958                   CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id,                     &
3959                        31, comm_inter, ierr )
3960                   CALL MPI_RECV( dtf, 1, MPI_REAL,                                  &
3961                        target_id, 32, comm_inter, status, ierr )
3962
3963                endif
3964                CALL MPI_BCAST( dtf, 1,   MPI_REAL,    0, comm2d, ierr ) 
3965         ELSE
3966                dtf = dt_3d
3967                if (myid == 0) then
3968                   CALL MPI_RECV( dtc, 1, MPI_REAL,                                  &
3969                        target_id, 31, comm_inter, status, ierr )
3970                   CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id,                     &
3971                        32, comm_inter, ierr )
3972
3973                endif
3974                CALL MPI_BCAST( dtc, 1,   MPI_REAL,    0, comm2d, ierr ) 
3975         
3976         ENDIF
3977!-- Identical timestep for coarse and fine grids
3978          dt_3d = MIN( dtc, dtf )
3979          !> @fixme setting old_dt might be obsolete at this point
3980          !>   Due to changes in timestep routine, setting of old_dt might be
3981          !>   not necessary any more at this point. However, could not be
3982          !>   tested so far.
3983          !>   2018-05-18, gronemeier
3984          old_dt = dt_3d
3985#endif
3986       END SUBROUTINE vnest_timestep_sync
3987       
3988       SUBROUTINE vnest_deallocate
3989#if defined( __parallel )
3990          USE control_parameters,                                                    &
3991              ONLY:  coupling_mode
3992             
3993          IMPLICIT NONE
3994         
3995          IF ( ALLOCATED(c_rnk_lst) ) DEALLOCATE (c_rnk_lst)
3996          IF ( ALLOCATED(f_rnk_lst) ) DEALLOCATE (f_rnk_lst)
3997         
3998          IF ( coupling_mode == 'vnested_crse')  THEN
3999             IF ( ALLOCATED (c2f_dims_cg) ) DEALLOCATE (c2f_dims_cg)
4000             IF ( ALLOCATED (f2c_dims_cg) ) DEALLOCATE (f2c_dims_cg)
4001          ELSEIF( coupling_mode == 'vnested_fine' )  THEN
4002             IF ( ALLOCATED (c2f_dims_fg) ) DEALLOCATE (c2f_dims_fg)
4003             IF ( ALLOCATED (f2c_dims_fg) ) DEALLOCATE (f2c_dims_fg)
4004          ENDIF
4005#endif
4006       END SUBROUTINE vnest_deallocate
4007
4008 END MODULE vertical_nesting_mod
Note: See TracBrowser for help on using the repository browser.