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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

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