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

Last change on this file since 4360 was 4360, checked in by suehring, 19 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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