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

Last change on this file since 2868 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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