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

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

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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