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

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