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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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