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

Last change on this file since 2712 was 2712, checked in by kanani, 4 years ago

Formatting and clean-up of vertical_nesting_mod

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