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

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

Merge of branch palm4u into trunk

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