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

Last change on this file since 3232 was 3232, checked in by raasch, 3 years ago

references to mrun replaced by palmrun, and updated

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