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

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

Vertical nesting implemented (SadiqHuq?)

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