source: palm/trunk/SOURCE/pmc_interface.f90 @ 1766

Last change on this file since 1766 was 1766, checked in by raasch, 6 years ago

pmc now runs with pointer version too

  • Property svn:keywords set to Id
File size: 146.5 KB
Line 
1MODULE pmc_interface
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! modifications to allow for using PALM's pointer version,
23! +new routine pmci_set_swaplevel
24!
25! Former revisions:
26! -----------------
27! $Id: pmc_interface.f90 1766 2016-02-29 08:37:15Z raasch $
28!
29! 1764 2016-02-28 12:45:19Z raasch
30! +cpl_parent_id,
31! cpp-statements for nesting replaced by __parallel statements,
32! errors output with message-subroutine,
33! index bugfixes in pmci_interp_tril_all,
34! some adjustments to PALM style
35!
36! 1762 2016-02-25 12:31:13Z hellstea
37! Initial revision by A. Hellsten
38!
39! Description:
40! ------------
41! Domain nesting interface routines. The low-level inter-domain communication   
42! is conducted by the PMC-library routines.
43!------------------------------------------------------------------------------!
44
45#if defined( __nopointer )
46!-- TO_DO: is it really required to use arrays te_m, tu_m, etc. and to set them
47!--        to zero withi buildings?
48    USE arrays_3d,                                                             &
49        ONLY:  dzu, dzw, e, e_p, pt, pt_p, q, q_p, te_m, tu_m, tv_m, tw_m, u,  &
50               u_p, v, v_p, w, w_p, zu, zw, z0
51#else
52   USE arrays_3d,                                                              &
53        ONLY:  dzu, dzw, e, e_p, e_1, e_2, pt, pt_p, pt_1, pt_2, q, q_p, q_1, q_2, te_m, tu_m, &
54               tv_m, tw_m, u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw, z0
55#endif
56
57    USE control_parameters,                                                    &
58        ONLY:  dt_3d, dz, humidity, message_string, nest_bound_l,              &
59               nest_bound_r, nest_bound_s, nest_bound_n, passive_scalar,       &
60               simulated_time, topography, volume_flow
61
62    USE cpulog,                                                                &
63        ONLY:  cpu_log, log_point_s
64
65    USE grid_variables,                                                        &
66        ONLY:  dx, dy
67
68    USE indices,                                                               &
69        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
70               nysv, nz, nzb, nzb_s_inner, nzb_u_inner, nzb_u_outer,           &
71               nzb_v_inner, nzb_v_outer, nzb_w_inner, nzb_w_outer, nzt
72
73    USE kinds
74
75#if defined( __parallel )
76#if defined( __lc )
77    USE MPI
78#else
79    INCLUDE "mpif.h"
80#endif
81
82    USE pegrid,                                                                &
83        ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
84               numprocs
85
86    USE pmc_client,                                                            &
87        ONLY:  pmc_clientinit, pmc_c_getnextarray, pmc_c_get_2d_index_list,    &
88               pmc_c_getbuffer, pmc_c_putbuffer, pmc_c_setind_and_allocmem,    &
89               pmc_c_set_dataarray, pmc_set_dataarray_name
90
91    USE pmc_general,                                                           &
92        ONLY:  da_namelen, pmc_max_modell, pmc_status_ok
93
94    USE pmc_handle_communicator,                                               &
95        ONLY:  pmc_get_local_model_info, pmc_init_model, pmc_is_rootmodel,     &
96               pmc_no_namelist_found, pmc_server_for_client
97
98    USE pmc_mpi_wrapper,                                                       &
99        ONLY:  pmc_bcast, pmc_recv_from_client, pmc_recv_from_server,          &
100               pmc_send_to_client, pmc_send_to_server
101
102    USE pmc_server,                                                            &
103        ONLY:  pmc_serverinit, pmc_s_fillbuffer, pmc_s_getdata_from_buffer,    &
104               pmc_s_getnextarray, pmc_s_setind_and_allocmem,                  &
105               pmc_s_set_active_data_array, pmc_s_set_dataarray,               &
106               pmc_s_set_2d_index_list
107
108#endif
109
110    IMPLICIT NONE
111
112!-- TO_DO: a lot of lines (including comments) in this file exceed the 80 char
113!--        limit. Try to reduce as much as possible
114
115!-- TO_DO: shouldn't we use public as default here? Only a minority of the
116!-- variables is private.
117    PRIVATE    !:  Note that the default publicity is here set to private.
118
119!
120!-- Constants
121    INTEGER(iwp), PARAMETER, PUBLIC ::  client_to_server = 2   !:
122    INTEGER(iwp), PARAMETER, PUBLIC ::  server_to_client = 1   !:
123
124!
125!-- Coupler setup
126    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_id  = 1            !:
127    CHARACTER(LEN=32), PUBLIC, SAVE ::  cpl_name               !:
128    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_npex               !:
129    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_npey               !:
130    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_parent_id          !:
131
132!
133!-- Control parameters, will be made input parameters later
134    CHARACTER(LEN=7), PUBLIC, SAVE ::  nesting_mode = 'two-way'  !: steering parameter for one- or two-way nesting
135
136    LOGICAL, PUBLIC, SAVE ::  nested_run = .FALSE.  !: general switch if nested run or not
137
138    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_l = -1.0_wp   !:
139    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_r = -1.0_wp   !:
140    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_s = -1.0_wp   !:
141    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_n = -1.0_wp   !:
142    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_t = -1.0_wp   !:
143
144!
145!-- Geometry
146    REAL(wp), PUBLIC, SAVE                    ::  area_t               !:
147    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_x              !:
148    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_y              !:
149    REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_x   !:
150    REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_y   !:
151
152!
153!-- Client coarse data arrays
154    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ec   !:
155    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ptc  !:
156    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  uc   !:
157    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  vc   !:
158    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  wc   !:
159    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  qc   !:
160
161    INTEGER(iwp), DIMENSION(5)                          ::  coarse_bound   !:
162    REAL(wp), PUBLIC, SAVE                              ::  xexl           !:
163    REAL(wp), PUBLIC, SAVE                              ::  xexr           !:
164    REAL(wp), PUBLIC, SAVE                              ::  yexs           !:
165    REAL(wp), PUBLIC, SAVE                              ::  yexn           !:
166    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !:
167    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !:
168    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !:
169    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !:
170    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !:
171
172!
173!-- Client interpolation coefficients and client-array indices to be precomputed and stored.
174    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !:
175    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !:
176    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !:
177    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !:
178    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !:
179    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !:
180    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !:
181    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !:
182    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !:
183    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !:
184    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !:
185    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !:
186    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !:
187    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !:
188    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !:
189    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !:
190    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
191    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
192
193!
194!-- Client index arrays and log-ratio arrays for the log-law near-wall corrections.
195!-- These are not truly 3-D arrays but multiply 2-D arrays.
196    INTEGER(iwp), PUBLIC, SAVE :: ncorr  !: ncorr is the 4th dimension of the log_ratio-arrays
197    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !:
198    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !:
199    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !:
200    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !:
201    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !:
202    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !:
203    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !:
204    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !:
205    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !:
206    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !:
207    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !:
208    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !:
209    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !:
210    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !:
211    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !:
212    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !:
213    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !:
214    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !:
215    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !:
216    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !:
217    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !:
218    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !:
219    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
220    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
221
222!
223!-- Upper bounds for k in anterpolation.
224    INTEGER(iwp), PUBLIC, SAVE ::  kceu   !:
225    INTEGER(iwp), PUBLIC, SAVE ::  kcew   !:
226
227!
228!-- Upper bound for k in log-law correction in interpolation.
229    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_l   !:
230    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_n   !:
231    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_r   !:
232    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_s   !:
233
234!
235!-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
236    INTEGER(iwp), PUBLIC, SAVE ::  nhll   !:
237    INTEGER(iwp), PUBLIC, SAVE ::  nhlr   !:
238    INTEGER(iwp), PUBLIC, SAVE ::  nhls   !:
239    INTEGER(iwp), PUBLIC, SAVE ::  nhln   !:
240
241!
242!-- Spatial under-relaxation coefficients for anterpolation.
243    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !:
244    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
245    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
246
247!
248!-- Client-array indices to be precomputed and stored for anterpolation.
249    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !:
250    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !:
251    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !:
252    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !:
253    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !:
254    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !:
255    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !:
256    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !:
257    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !:
258    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !:
259    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
260    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
261
262!
263!-- Module private variables.
264    INTEGER(iwp), DIMENSION(3)          ::  define_coarse_grid_int    !:
265    REAL(wp), DIMENSION(9)              ::  define_coarse_grid_real   !:
266
267    TYPE coarsegrid_def
268       INTEGER(iwp)                        ::  nx
269       INTEGER(iwp)                        ::  ny
270       INTEGER(iwp)                        ::  nz
271       REAL(wp)                            ::  dx
272       REAL(wp)                            ::  dy
273       REAL(wp)                            ::  dz
274       REAL(wp)                            ::  lower_left_coord_x
275       REAL(wp)                            ::  lower_left_coord_y
276       REAL(wp)                            ::  xend
277       REAL(wp)                            ::  yend
278       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_x
279       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_y
280       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu       
281       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw       
282       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu       
283       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw       
284    END TYPE coarsegrid_def
285                                         
286    TYPE(coarsegrid_def), SAVE             ::  cg   !:
287
288
289    INTERFACE pmci_client_datatrans
290       MODULE PROCEDURE pmci_client_datatrans
291    END INTERFACE
292
293    INTERFACE pmci_client_initialize
294       MODULE PROCEDURE pmci_client_initialize
295    END INTERFACE
296
297    INTERFACE pmci_client_synchronize
298       MODULE PROCEDURE pmci_client_synchronize
299    END INTERFACE
300
301    INTERFACE pmci_ensure_nest_mass_conservation
302       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
303    END INTERFACE
304
305    INTERFACE pmci_init
306       MODULE PROCEDURE pmci_init
307    END INTERFACE
308
309    INTERFACE pmci_modelconfiguration
310       MODULE PROCEDURE pmci_modelconfiguration
311    END INTERFACE
312
313    INTERFACE pmci_server_initialize
314       MODULE PROCEDURE pmci_server_initialize
315    END INTERFACE
316
317    INTERFACE pmci_server_synchronize
318       MODULE PROCEDURE pmci_server_synchronize
319    END INTERFACE
320
321    INTERFACE pmci_set_swaplevel
322       MODULE PROCEDURE pmci_set_swaplevel
323    END INTERFACE pmci_set_swaplevel
324
325    INTERFACE pmci_update_new
326       MODULE PROCEDURE pmci_update_new
327    END INTERFACE
328
329    PUBLIC pmci_client_datatrans
330    PUBLIC pmci_client_initialize
331    PUBLIC pmci_client_synchronize
332    PUBLIC pmci_ensure_nest_mass_conservation
333    PUBLIC pmci_init
334    PUBLIC pmci_modelconfiguration
335    PUBLIC pmci_server_datatrans
336    PUBLIC pmci_server_initialize
337    PUBLIC pmci_server_synchronize
338    PUBLIC pmci_set_swaplevel
339    PUBLIC pmci_update_new
340
341
342 CONTAINS
343
344
345 SUBROUTINE pmci_init( world_comm )
346
347    IMPLICIT NONE
348
349    INTEGER, INTENT(OUT) ::  world_comm   !:
350
351#if defined( __parallel )
352
353    INTEGER(iwp)         ::  ierr         !:
354    INTEGER(iwp)         ::  istat        !:
355    INTEGER(iwp)         ::  pmc_status   !:
356
357
358    CALL pmc_init_model( world_comm, nesting_mode, pmc_status )
359
360    IF ( pmc_status == pmc_no_namelist_found )  THEN
361!
362!--    This is not a nested run
363!
364!--    TO_DO: this wouldn't be required any more?
365       world_comm = MPI_COMM_WORLD
366       cpl_id     = 1
367       cpl_name   = ""
368       cpl_npex   = 2
369       cpl_npey   = 2
370       lower_left_coord_x = 0.0_wp
371       lower_left_coord_y = 0.0_wp
372       RETURN
373    ELSE
374!
375!--    Set the general steering switch which tells PALM that its a nested run
376       nested_run = .TRUE.
377    ENDIF
378
379    CALL pmc_get_local_model_info( my_cpl_id = cpl_id,                         &
380                                   my_cpl_parent_id = cpl_parent_id,           &
381                                   cpl_name = cpl_name,                        &
382                                   npe_x = cpl_npex, npe_y = cpl_npey,         &
383                                   lower_left_x = lower_left_coord_x,          &
384                                   lower_left_y = lower_left_coord_y )
385!
386!-- Message that communicators for nesting are initialized.
387!-- Attention: myid has been set at the end of pmc_init_model in order to
388!-- guarantee that only PE0 of the root domain does the output.
389    CALL location_message( 'finished', .TRUE. )
390!
391!-- Reset myid to its default value
392    myid = 0
393#else
394!
395!-- Nesting cannot be used in serial mode. cpl_id is set to root domain (1)
396!-- because no location messages would be generated otherwise.
397!-- world_comm is given a dummy value to avoid compiler warnings (INTENT(OUT)
398!-- should get an explicit value)
399    cpl_id     = 1
400    nested_run = .FALSE.
401    world_comm = 1
402#endif
403
404 END SUBROUTINE pmci_init
405
406
407
408 SUBROUTINE pmci_modelconfiguration
409
410    IMPLICIT NONE
411
412    CALL location_message( 'setup the nested model configuration', .FALSE. )
413    CALL pmci_setup_coordinates   !:  Compute absolute coordinates valid for all models
414    CALL pmci_setup_client        !:  Initialize PMC Client (Must be called before pmc_setup_server)
415    CALL pmci_setup_server        !:  Initialize PMC Server
416    CALL location_message( 'finished', .TRUE. )
417
418 END SUBROUTINE pmci_modelconfiguration
419
420
421
422 SUBROUTINE pmci_setup_server
423
424#if defined( __parallel )
425    IMPLICIT NONE
426
427    INTEGER(iwp)               ::  client_id        !:
428    INTEGER(iwp)               ::  ierr             !:
429    INTEGER(iwp)               ::  i                !:
430    INTEGER(iwp)               ::  j                !:
431    INTEGER(iwp)               ::  k                !:
432    INTEGER(iwp)               ::  m                !:
433    INTEGER(iwp)               ::  nomatch          !:
434    INTEGER(iwp)               ::  nx_cl            !:
435    INTEGER(iwp)               ::  ny_cl            !:
436    INTEGER(iwp)               ::  nz_cl            !:
437    INTEGER(iwp), DIMENSION(5) ::  val              !:
438    REAL(wp), DIMENSION(1)     ::  fval             !:
439    REAL(wp)                   ::  dx_cl            !:
440    REAL(wp)                   ::  dy_cl            !:
441    REAL(wp)                   ::  xez              !:
442    REAL(wp)                   ::  yez              !:
443    CHARACTER(len=32)          ::  mychannel
444    CHARACTER(len=32)          ::  myname
445    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !:
446    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !:
447   
448
449!
450!   Initialize the PMC server
451    CALL pmc_serverinit
452
453!
454!-- Get coordinates from all clients
455    DO  m = 1, SIZE( pmc_server_for_client ) - 1
456       client_id = pmc_server_for_client(m)
457       IF ( myid == 0 )  THEN       
458
459          CALL pmc_recv_from_client( client_id, val,  size(val),  0, 123, ierr )
460          CALL pmc_recv_from_client( client_id, fval, size(fval), 0, 124, ierr )
461         
462          nx_cl = val(1)
463          ny_cl = val(2)
464          dx_cl = val(4)
465          dy_cl = val(5)
466
467          nz_cl = nz
468
469!
470!--       Find the highest client level in the coarse grid for the reduced z
471!--       transfer
472          DO  k = 1, nz                 
473             IF ( zw(k) > fval(1) )  THEN
474                nz_cl = k
475                EXIT
476             ENDIF
477          ENDDO
478
479!   
480!--       Get absolute coordinates from the client
481          ALLOCATE( cl_coord_x(-nbgp:nx_cl+nbgp) )
482          ALLOCATE( cl_coord_y(-nbgp:ny_cl+nbgp) )
483         
484          CALL pmc_recv_from_client( client_id, cl_coord_x, SIZE( cl_coord_x ),&
485                                     0, 11, ierr )
486          CALL pmc_recv_from_client( client_id, cl_coord_y, SIZE( cl_coord_y ),&
487                                     0, 12, ierr )
488          WRITE ( 0, * )  'receive from pmc Client ', client_id, nx_cl, ny_cl
489         
490          define_coarse_grid_real(1) = lower_left_coord_x
491          define_coarse_grid_real(2) = lower_left_coord_y
492!--       TO_DO: remove this?
493          define_coarse_grid_real(3) = 0             !  KK currently not used.
494          define_coarse_grid_real(4) = 0
495          define_coarse_grid_real(5) = dx
496          define_coarse_grid_real(6) = dy
497          define_coarse_grid_real(7) = lower_left_coord_x + ( nx + 1 ) * dx
498          define_coarse_grid_real(8) = lower_left_coord_y + ( ny + 1 ) * dy
499          define_coarse_grid_real(9) = dz
500
501          define_coarse_grid_int(1)  = nx
502          define_coarse_grid_int(2)  = ny
503          define_coarse_grid_int(3)  = nz_cl
504
505!
506!--       Check that the client domain is completely inside the server domain.
507          nomatch = 0
508          xez = ( nbgp + 1 ) * dx
509          yez = ( nbgp + 1 ) * dy
510          IF ( cl_coord_x(0) < define_coarse_grid_real(1) + xez )          nomatch = 1 
511          IF ( cl_coord_x(nx_cl + 1) > define_coarse_grid_real(7) - xez )  nomatch = 1 
512          IF ( cl_coord_y(0) < define_coarse_grid_real(2) + yez )          nomatch = 1 
513          IF ( cl_coord_y(ny_cl + 1) > define_coarse_grid_real(8) - yez )  nomatch = 1
514
515          DEALLOCATE( cl_coord_x )
516          DEALLOCATE( cl_coord_y )
517
518!
519!--       Send coarse grid information to client
520          CALL pmc_send_to_client( client_id, Define_coarse_grid_real, 9, 0,   &
521                                   21, ierr )
522          CALL pmc_send_to_client( client_id, Define_coarse_grid_int,  3, 0,   &
523                                   22, ierr )
524
525!
526!--       Send local grid to client
527          CALL pmc_send_to_client( client_id, coord_x, nx+1+2*nbgp, 0, 24, ierr )
528          CALL pmc_send_to_client( client_id, coord_y, ny+1+2*nbgp, 0, 25, ierr )
529
530!
531!--       Also send the dzu-, dzw-, zu- and zw-arrays here
532          CALL pmc_send_to_client( client_id, dzu, nz_cl + 1, 0, 26, ierr )
533          CALL pmc_send_to_client( client_id, dzw, nz_cl + 1, 0, 27, ierr )
534          CALL pmc_send_to_client( client_id, zu,  nz_cl + 2, 0, 28, ierr )
535          CALL pmc_send_to_client( client_id, zw,  nz_cl + 2, 0, 29, ierr )
536
537       ENDIF
538
539       CALL MPI_BCAST( nomatch, 1, MPI_INTEGER, 0, comm2d, ierr )
540       IF ( nomatch /= 0 ) THEN
541          WRITE ( message_string, * )  'Error: nested client domain does ',    &
542                                       'not fit into its server domain'
543          CALL message( 'pmc_palm_setup_server', 'PA0XYZ', 1, 2, 0, 6, 0 )
544       ENDIF
545     
546       CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
547       
548       CALL pmci_create_index_list
549
550!
551!--    Include couple arrays into server content
552       DO  WHILE ( pmc_s_getnextarray( client_id, myname ) )
553          CALL pmci_set_array_pointer( myname, client_id = client_id,          &
554                                       nz_cl = nz_cl )
555       ENDDO
556       CALL pmc_s_setind_and_allocmem( client_id )
557    ENDDO
558
559 CONTAINS
560
561
562   SUBROUTINE pmci_create_index_list
563
564       IMPLICIT NONE
565
566       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  coarse_bound_all   !:
567       INTEGER(iwp)                               ::  i                  !:
568       INTEGER(iwp)                               ::  ic                 !:
569       INTEGER(iwp)                               ::  ierr               !:
570       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list         !:
571       INTEGER(iwp)                               ::  j                  !:
572       INTEGER(iwp)                               ::  k                  !:
573       INTEGER(iwp)                               ::  m                  !:
574       INTEGER(iwp)                               ::  n                  !:
575       INTEGER(iwp)                               ::  npx                !:
576       INTEGER(iwp)                               ::  npy                !:
577       INTEGER(iwp)                               ::  nrx                !:
578       INTEGER(iwp)                               ::  nry                !:
579       INTEGER(iwp)                               ::  px                 !:
580       INTEGER(iwp)                               ::  py                 !:
581       INTEGER(iwp), DIMENSION(2)                 ::  scoord             !:
582       INTEGER(iwp)                               ::  server_pe          !:
583       INTEGER(iwp), DIMENSION(2)                 ::  size_of_array      !:
584
585
586       IF ( myid == 0 )  THEN
587          CALL pmc_recv_from_client( client_id, size_of_array, 2, 0, 40, ierr )
588          ALLOCATE( coarse_bound_all(size_of_array(1),size_of_array(2)) )
589          CALL pmc_recv_from_client( client_id, coarse_bound_all,              &
590                                     SIZE( coarse_bound_all ), 0, 41, ierr )
591
592!
593!--       Compute size of index_list.
594          ic = 0
595          DO  k = 1, size_of_array(2)
596             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
597                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
598                   ic = ic + 1
599                ENDDO
600             ENDDO
601          ENDDO
602
603          ALLOCATE( index_list(6,ic) )
604
605          CALL MPI_COMM_SIZE( comm1dx, npx, ierr )
606          CALL MPI_COMM_SIZE( comm1dy, npy, ierr )
607
608          nrx = nxr - nxl + 1   !  +1 in index because FORTRAN indexing starts with 1, palm with 0
609          nry = nyn - nys + 1
610          ic  = 0
611          DO  k = 1, size_of_array(2)                                   !  loop over all client PEs
612             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)       !  area in y required by actual client PE
613                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)    !  area in x required by actual client PE
614                   px = i / nrx
615                   py = j / nry
616                   scoord(1) = px
617                   scoord(2) = py
618                   CALL MPI_CART_RANK( comm2d, scoord, server_pe, ierr )
619                 
620                   ic = ic + 1
621                   index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp       !  First index in Server Array
622                   index_list(2,ic) = j - ( py * nry ) + 1 + nbgp       !  Second index in Server Array
623                   index_list(3,ic) = i - coarse_bound_all(1,k) + 1     !  x Index client coarse grid
624                   index_list(4,ic) = j - coarse_bound_all(3,k) + 1     !  y Index client coarse grid
625                   index_list(5,ic) = k - 1                             !  PE Number client
626                   index_list(6,ic) = server_pe                         !  PE Number server
627                ENDDO
628             ENDDO
629          ENDDO
630          CALL pmc_s_set_2d_index_list( client_id, index_list(:,1:ic) )
631       ELSE
632          ALLOCATE( index_list(6,1) )    !  Dummy allocate
633          CALL pmc_s_set_2d_index_list( client_id, index_list )
634       ENDIF
635
636       DEALLOCATE(index_list)
637
638     END SUBROUTINE pmci_create_index_list
639
640#endif
641 END SUBROUTINE pmci_setup_server
642
643
644
645 SUBROUTINE pmci_setup_client
646
647#if defined( __parallel )
648    IMPLICIT NONE
649
650    CHARACTER(LEN=DA_Namelen)  ::  myname     !:
651
652    INTEGER(iwp)               ::  i          !:
653    INTEGER(iwp)               ::  ierr       !:
654    INTEGER(iwp)               ::  icl        !:
655    INTEGER(iwp)               ::  icr        !:
656    INTEGER(iwp)               ::  j          !:
657    INTEGER(iwp)               ::  jcn        !:
658    INTEGER(iwp)               ::  jcs        !:
659    INTEGER(iwp), DIMENSION(5) ::  val        !:
660   
661    REAL(wp)                   ::  area_t_l   !:
662    REAL(wp), DIMENSION(1)     ::  fval       !:
663    REAL(wp)                   ::  xcs        !:
664    REAL(wp)                   ::  xce        !:
665    REAL(wp)                   ::  ycs        !:
666    REAL(wp)                   ::  yce        !:
667    REAL(wp)                   ::  zct        !:
668    REAL(wp), DIMENSION(4)     ::  ztt        !:
669                                             
670
671!-- TO_DO: describe what is happening in this if-clause
672!-- Root Model does not have Server and is not a client
673    IF ( .NOT. pmc_is_rootmodel() )  THEN
674       CALL pmc_clientinit
675       
676       CALL pmc_set_dataarray_name( 'coarse', 'u'  ,'fine', 'u',  ierr )
677       CALL pmc_set_dataarray_name( 'coarse', 'v'  ,'fine', 'v',  ierr )
678       CALL pmc_set_dataarray_name( 'coarse', 'w'  ,'fine', 'w',  ierr )
679       CALL pmc_set_dataarray_name( 'coarse', 'e'  ,'fine', 'e',  ierr )
680       CALL pmc_set_dataarray_name( 'coarse', 'pt' ,'fine', 'pt', ierr )
681       IF ( humidity  .OR.  passive_scalar )  THEN
682          CALL pmc_set_dataarray_name( 'coarse', 'q'  ,'fine', 'q',  ierr )
683       ENDIF
684
685!
686!--    Update this list appropritely and also in create_client_arrays and in
687!--    pmci_set_array_pointer.
688!--    If a variable is removed, it only has to be removed from here.
689       CALL pmc_set_dataarray_name( lastentry = .TRUE. )
690
691!
692!--    Send grid to server
693       val(1)  = nx
694       val(2)  = ny
695       val(3)  = nz
696       val(4)  = dx
697       val(5)  = dy
698       fval(1) = zw(nzt+1)
699
700       IF ( myid == 0 )  THEN
701          CALL pmc_send_to_server( val, SIZE( val ), 0, 123, ierr )
702          CALL pmc_send_to_server( fval, SIZE( fval ), 0, 124, ierr )
703          CALL pmc_send_to_server( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
704          CALL pmc_send_to_server( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
705
706!
707!--       Receive Coarse grid information.
708          CALL pmc_recv_from_server( define_coarse_grid_real, 9, 0, 21, ierr )
709          CALL pmc_recv_from_server( define_coarse_grid_int,  3, 0, 22, ierr )
710
711!
712!--       Receive also the dz-,zu- and zw-arrays here.
713!--       TO_DO: what is the meaning of above comment + remove write statements
714!--              and give this informations in header
715          WRITE(0,*) 'Coarse grid from Server '
716          WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
717          WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
718          WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(7)
719          WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(8)
720          WRITE(0,*) 'dx            = ',define_coarse_grid_real(5)
721          WRITE(0,*) 'dy            = ',define_coarse_grid_real(6)
722          WRITE(0,*) 'dz            = ',define_coarse_grid_real(9)
723          WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
724          WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
725          WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)       
726       ENDIF
727
728       CALL MPI_BCAST( define_coarse_grid_real, 9, MPI_REAL, 0, comm2d, ierr )
729       CALL MPI_BCAST( define_coarse_grid_int, 3, MPI_INTEGER, 0, comm2d, ierr )
730
731       cg%dx = define_coarse_grid_real(5)
732       cg%dy = define_coarse_grid_real(6)
733       cg%dz = define_coarse_grid_real(9)
734       cg%nx = define_coarse_grid_int(1)
735       cg%ny = define_coarse_grid_int(2)
736       cg%nz = define_coarse_grid_int(3)
737
738!
739!--    Get Server coordinates on coarse grid
740       ALLOCATE( cg%coord_x(-nbgp:cg%nx+nbgp) )
741       ALLOCATE( cg%coord_y(-nbgp:cg%ny+nbgp) )
742     
743       ALLOCATE( cg%dzu(1:cg%nz+1) )
744       ALLOCATE( cg%dzw(1:cg%nz+1) )
745       ALLOCATE( cg%zu(0:cg%nz+1) )
746       ALLOCATE( cg%zw(0:cg%nz+1) )
747
748!
749!--    Get coarse grid coordinates and vales of the z-direction from server
750       IF ( myid == 0) THEN
751          CALL pmc_recv_from_server( cg%coord_x, cg%nx + 1 + 2 * nbgp, 0, 24,  &
752                                     ierr )
753          CALL pmc_recv_from_server( cg%coord_y, cg%ny + 1 + 2 * nbgp, 0, 25,  &
754                                     ierr )
755          CALL pmc_recv_from_server( cg%dzu, cg%nz + 1, 0, 26, ierr )
756          CALL pmc_recv_from_server( cg%dzw, cg%nz + 1, 0, 27, ierr )
757          CALL pmc_recv_from_server( cg%zu, cg%nz + 2, 0, 28, ierr )
758          CALL pmc_recv_from_server( cg%zw, cg%nz + 2, 0, 29, ierr )
759       ENDIF
760
761!
762!--    and broadcast this information
763       CALL MPI_BCAST( cg%coord_x, cg%nx + 1 + 2 * nbgp, MPI_REAL, 0, comm2d,  &
764                       ierr )
765       CALL MPI_BCAST( cg%coord_y, cg%ny + 1 + 2 * nbgp, MPI_REAL, 0, comm2d,  &
766                       ierr )
767       CALL MPI_BCAST( cg%dzu, cg%nz + 1, MPI_REAL, 0, comm2d, ierr )
768       CALL MPI_BCAST( cg%dzw, cg%nz + 1, MPI_REAL, 0, comm2d, ierr )
769       CALL MPI_BCAST( cg%zu, cg%nz + 2,  MPI_REAL, 0, comm2d, ierr )
770       CALL MPI_BCAST( cg%zw, cg%nz + 2,  MPI_REAL, 0, comm2d, ierr )
771       
772       CALL pmci_map_fine_to_coarse_grid
773       CALL pmc_c_get_2d_index_list
774
775!
776!--    Include couple arrays into client content.
777       DO  WHILE ( pmc_c_getnextarray( myname ) )
778!--       TO_DO: Klaus, why the c-arrays are still up to cg%nz??
779          CALL pmci_create_client_arrays ( myname, icl, icr, jcs, jcn, cg%nz )
780       ENDDO
781       CALL pmc_c_setind_and_allocmem
782
783!
784!--    Precompute interpolation coefficients and client-array indices
785       CALL pmci_init_interp_tril
786
787!
788!--    Precompute the log-law correction index- and ratio-arrays
789       CALL pmci_init_loglaw_correction
790
791!
792!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio
793       CALL pmci_init_tkefactor
794
795!
796!--    Two-way coupling     
797       IF ( nesting_mode == 'two-way' )  THEN
798          CALL pmci_init_anterp_tophat
799       ENDIF
800
801!
802!--    Finally, compute the total area of the top-boundary face of the domain.
803!--    This is needed in the pmc_ensure_nest_mass_conservation     
804       area_t_l = ( nxr + 1 - nxl ) * ( nyn + 1 - nys ) * dx * dy
805       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
806       CALL MPI_ALLREDUCE( area_t_l, area_t, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
807
808!
809!--    TO_DO: Why not just simply? test this!
810       !area_t_l = ( nx + 1 ) * (ny + 1 ) * dx * dy           
811
812    ENDIF
813
814 CONTAINS
815
816
817    SUBROUTINE pmci_map_fine_to_coarse_grid
818
819        IMPLICIT NONE
820
821        INTEGER(iwp), DIMENSION(5,numprocs)  ::  coarse_bound_all   !:
822        INTEGER(iwp), DIMENSION(2)           ::  size_of_array      !:
823                                             
824        REAL(wp)                             ::  coarse_dx   !:
825        REAL(wp)                             ::  coarse_dy   !:
826        REAL(wp)                             ::  loffset     !:
827        REAL(wp)                             ::  noffset     !:
828        REAL(wp)                             ::  roffset     !:
829        REAL(wp)                             ::  soffset     !:
830
831!
832!--     Determine indices of interpolation/anterpolation area in coarse grid.
833        coarse_dx = cg%dx
834        coarse_dy = cg%dy
835       
836        loffset = MOD( coord_x(nxl), coarse_dx )       !  If the fine- and coarse grid nodes do not match.
837        xexl    = coarse_dx + loffset
838        nhll    = CEILING( xexl / coarse_dx )          !  This is needed in the anterpolation phase.
839        xcs     = coord_x(nxl) - xexl
840        DO  i = 0, cg%nx
841           IF ( cg%coord_x(i) > xcs )  THEN
842              icl = MAX( -1, i-1 ) 
843              EXIT
844           ENDIF
845        ENDDO
846
847        roffset = MOD( coord_x(nxr + 1), coarse_dx )   !  If the fine- and coarse grid nodes do not match.
848        xexr    = coarse_dx + roffset
849        nhlr    = CEILING( xexr / coarse_dx )          !  This is needed in the anterpolation phase.
850        xce     = coord_x(nxr) + xexr
851        DO  i = cg%nx, 0 , -1
852           IF ( cg%coord_x(i) < xce )  THEN
853              icr = MIN( cg%nx + 1, i + 1 )
854              EXIT
855           ENDIF
856        ENDDO
857
858        soffset = MOD( coord_y(nys), coarse_dy )       !  If the fine- and coarse grid nodes do not match
859        yexs    = coarse_dy + soffset
860        nhls    = CEILING( yexs / coarse_dy )          !  This is needed in the anterpolation phase.
861        ycs     = coord_y(nys) - yexs
862        DO  j = 0, cg%ny
863           IF ( cg%coord_y(j) > ycs )  THEN
864              jcs = MAX( -nbgp, j - 1 ) 
865              EXIT
866           ENDIF
867        ENDDO
868
869        noffset = MOD( coord_y(nyn + 1), coarse_dy )   !  If the fine- and coarse grid nodes do not match
870        yexn    = coarse_dy + noffset
871        nhln    = CEILING( yexn / coarse_dy )          !  This is needed in the anterpolation phase.
872        yce     = coord_y(nyn) + yexn
873        DO  j = cg%ny, 0, -1
874           IF ( cg%coord_y(j) < yce )  THEN
875              jcn = MIN( cg%ny + nbgp, j + 1 ) 
876              EXIT
877           ENDIF
878        ENDDO
879
880        WRITE( 0, * )  'Coarse area ', myid, icl, icr, jcs, jcn
881
882        coarse_bound(1) = icl
883        coarse_bound(2) = icr
884        coarse_bound(3) = jcs
885        coarse_bound(4) = jcn
886        coarse_bound(5) = myid
887!
888!--     Note that MPI_Gather receives data from all processes in the rank order.
889        CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5, &
890                         MPI_INTEGER, 0, comm2d, ierr )
891
892        IF ( myid == 0 )  THEN
893           size_of_array(1) = SIZE( coarse_bound_all, 1 )
894           size_of_array(2) = SIZE( coarse_bound_all, 2 )
895           CALL pmc_send_to_server( size_of_array, 2, 0, 40, ierr )
896           CALL pmc_send_to_server( coarse_bound_all, SIZE( coarse_bound_all ), 0, 41, ierr )
897        ENDIF
898
899   END SUBROUTINE pmci_map_fine_to_coarse_grid
900
901
902
903   SUBROUTINE pmci_init_interp_tril
904
905!
906!--   Precomputation of the interpolation coefficients and client-array indices
907!--   to be used by the interpolation routines interp_tril_lr, interp_tril_ns and
908!--   interp_tril_t. Constant dz is still assumed.
909!
910!--                                  Antti Hellsten 3.3.2015.
911!
912!--   Modified for variable dz, but not yet tested.
913!--                                  Antti Hellsten 27.3.2015.
914!
915      IMPLICIT NONE
916
917      INTEGER(iwp) ::  i       !:
918      INTEGER(iwp) ::  i1      !:
919      INTEGER(iwp) ::  j       !:
920      INTEGER(iwp) ::  j1      !:
921      INTEGER(iwp) ::  k       !:
922      INTEGER(iwp) ::  kc      !:
923
924      REAL(wp) ::  coarse_dx   !:
925      REAL(wp) ::  coarse_dy   !:
926      REAL(wp) ::  coarse_dz   !:
927      REAL(wp) ::  xb          !:
928      REAL(wp) ::  xcsu        !:
929      REAL(wp) ::  xfso        !:
930      REAL(wp) ::  xcso        !:
931      REAL(wp) ::  xfsu        !:
932      REAL(wp) ::  yb          !:
933      REAL(wp) ::  ycso        !:
934      REAL(wp) ::  ycsv        !:
935      REAL(wp) ::  yfso        !:
936      REAL(wp) ::  yfsv        !:
937      REAL(wp) ::  zcso        !:
938      REAL(wp) ::  zcsw        !:
939      REAL(wp) ::  zfso        !:
940      REAL(wp) ::  zfsw        !:
941     
942
943      coarse_dx = cg%dx
944      coarse_dy = cg%dy
945      coarse_dz = cg%dz
946      xb        = nxl * dx
947      yb        = nys * dy
948     
949      ALLOCATE( icu(nxlg:nxrg) )
950      ALLOCATE( ico(nxlg:nxrg) )
951      ALLOCATE( jcv(nysg:nyng) )
952      ALLOCATE( jco(nysg:nyng) )
953      ALLOCATE( kcw(nzb:nzt+1) )
954      ALLOCATE( kco(nzb:nzt+1) )
955      ALLOCATE( r1xu(nxlg:nxrg) )
956      ALLOCATE( r2xu(nxlg:nxrg) )
957      ALLOCATE( r1xo(nxlg:nxrg) )
958      ALLOCATE( r2xo(nxlg:nxrg) )
959      ALLOCATE( r1yv(nysg:nyng) )
960      ALLOCATE( r2yv(nysg:nyng) )
961      ALLOCATE( r1yo(nysg:nyng) )
962      ALLOCATE( r2yo(nysg:nyng) )
963      ALLOCATE( r1zw(nzb:nzt+1) )
964      ALLOCATE( r2zw(nzb:nzt+1) )
965      ALLOCATE( r1zo(nzb:nzt+1) )
966      ALLOCATE( r2zo(nzb:nzt+1) )
967
968!
969!--   Note that the node coordinates xfs... and xcs... are relative
970!--   to the lower-left-bottom corner of the fc-array, not the actual
971!--   client domain corner.
972      DO  i = nxlg, nxrg
973         xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
974         xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
975         icu(i)  = icl + FLOOR( xfsu / coarse_dx )
976         ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * coarse_dx ) / coarse_dx )
977         xcsu    = ( icu(i) - icl ) * coarse_dx
978         xcso    = ( ico(i) - icl ) * coarse_dx + 0.5_wp * coarse_dx
979         r2xu(i) = ( xfsu - xcsu ) / coarse_dx
980         r2xo(i) = ( xfso - xcso ) / coarse_dx
981         r1xu(i) = 1.0_wp - r2xu(i)
982         r1xo(i) = 1.0_wp - r2xo(i)
983      ENDDO
984
985      DO  j = nysg, nyng
986         yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
987         yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
988         jcv(j)  = jcs + FLOOR( yfsv/coarse_dy )           
989         jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * coarse_dy ) / coarse_dy )
990         ycsv    = ( jcv(j) - jcs ) * coarse_dy
991         ycso    = ( jco(j) - jcs ) * coarse_dy + 0.5_wp * coarse_dy         
992         r2yv(j) = ( yfsv - ycsv ) / coarse_dy           
993         r2yo(j) = ( yfso - ycso ) / coarse_dy 
994         r1yv(j) = 1.0_wp - r2yv(j)
995         r1yo(j) = 1.0_wp - r2yo(j)       
996      ENDDO
997
998      DO  k = nzb, nzt + 1
999         zfsw = zw(k)
1000         zfso = zu(k)
1001
1002         kc = 0
1003         DO  WHILE ( cg%zw(kc) <= zfsw ) 
1004            kc = kc + 1
1005         ENDDO
1006         kcw(k) = kc - 1
1007         
1008         kc = 0
1009         DO  WHILE ( cg%zu(kc) <= zfso ) 
1010            kc = kc + 1
1011         ENDDO
1012         kco(k) = kc - 1
1013
1014         zcsw    = cg%zw(kcw(k))
1015         zcso    = cg%zu(kco(k))
1016         r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k) + 1 )
1017         r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k) + 1 )
1018         r1zw(k) = 1.0_wp - r2zw(k)
1019         r1zo(k) = 1.0_wp - r2zo(k)
1020      ENDDO
1021     
1022   END SUBROUTINE pmci_init_interp_tril
1023
1024
1025
1026   SUBROUTINE pmci_init_loglaw_correction
1027
1028!
1029!--   Precomputation of the index and log-ratio arrays for the log-law corrections
1030!--   for near-wall nodes after the nest-BC interpolation.
1031!--   These are used by the interpolation routines interp_tril_lr and interp_tril_ns.
1032!
1033!--                                  Antti Hellsten 30.12.2015.
1034!
1035      IMPLICIT NONE
1036      INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i. 
1037      INTEGER(iwp) ::  i            !:
1038      INTEGER(iwp) ::  icorr        !:
1039      INTEGER(iwp) ::  inc          !:  Wall outward-normal index increment -1 or 1, for direction=1, inc=1 always.
1040      INTEGER(iwp) ::  iw           !:
1041      INTEGER(iwp) ::  j            !:
1042      INTEGER(iwp) ::  jcorr        !:
1043      INTEGER(iwp) ::  jw           !:
1044      INTEGER(iwp) ::  k            !:
1045      INTEGER(iwp) ::  kb           !:
1046      INTEGER(iwp) ::  kcorr        !:
1047      INTEGER(iwp) ::  lc           !:
1048      INTEGER(iwp) ::  ni           !:
1049      INTEGER(iwp) ::  nj           !:
1050      INTEGER(iwp) ::  nk           !:
1051      INTEGER(iwp) ::  nzt_topo_max !:
1052      INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
1053
1054      REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
1055      LOGICAL  ::  more             !:
1056
1057!
1058!--   First determine the maximum k-index needed for the near-wall corrections.
1059!--   This maximum is individual for each boundary to minimize the storage requirements
1060!--   and to minimize the corresponding loop k-range in the interpolation routines.
1061      nzt_topo_nestbc_l = nzb
1062      IF ( nest_bound_l )  THEN
1063         DO  i = nxl - 1, nxl
1064            DO  j = nys, nyn
1065               nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),   &
1066                                        nzb_v_inner(j,i), nzb_w_inner(j,i) )
1067            ENDDO
1068         ENDDO
1069         nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
1070      ENDIF
1071     
1072      nzt_topo_nestbc_r = nzb
1073      IF ( nest_bound_r )  THEN
1074         i = nxr + 1
1075         DO  j = nys, nyn
1076            nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),      &
1077                                     nzb_v_inner(j,i), nzb_w_inner(j,i) )
1078         ENDDO
1079         nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1       
1080      ENDIF
1081
1082      nzt_topo_nestbc_s = nzb
1083      IF ( nest_bound_s )  THEN
1084         DO  j = nys - 1, nys
1085            DO  i = nxl, nxr
1086               nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),   &
1087                                        nzb_v_inner(j,i), nzb_w_inner(j,i) )
1088            ENDDO
1089         ENDDO
1090         nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1     
1091      ENDIF
1092
1093      nzt_topo_nestbc_n = nzb
1094      IF ( nest_bound_n )  THEN
1095         j = nyn + 1
1096         DO  i = nxl, nxr
1097            nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),      &
1098                                     nzb_v_inner(j,i), nzb_w_inner(j,i) )
1099         ENDDO
1100         nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
1101      ENDIF
1102
1103!
1104!--   Then determine the maximum number of near-wall nodes per wall point based
1105!--   on the grid-spacing ratios.
1106      nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,                &
1107                          nzt_topo_nestbc_s, nzt_topo_nestbc_n )
1108
1109!
1110!--   Note that the outer division must be integer division.
1111      ni = CEILING( cg%dx / dx ) / 2
1112      nj = CEILING( cg%dy / dy ) / 2
1113      nk = 1
1114      DO  k = 1, nzt_topo_max
1115         nk = MAX( nk, CEILING( cg%dzu(k) / dzu(k) ) )
1116      ENDDO
1117      nk = nk / 2   !  Note that this must be integer division.
1118      ncorr =  MAX( ni, nj, nk )
1119
1120      ALLOCATE( lcr(0:ncorr - 1) )
1121      lcr = 1.0_wp
1122
1123!
1124!--   First horizontal walls
1125!--   Left boundary
1126      IF ( nest_bound_l )   THEN
1127         ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
1128         ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
1129         ALLOCATE( logc_ratio_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
1130         ALLOCATE( logc_ratio_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
1131         logc_u_l       = 0
1132         logc_v_l       = 0
1133         logc_ratio_u_l = 1.0_wp
1134         logc_ratio_v_l = 1.0_wp
1135         direction      = 1
1136         inc            = 1
1137
1138         DO  j = nys, nyn     
1139!
1140!--         Left boundary for u
1141            i   = 0
1142            kb  = nzb_u_inner(j,i)
1143            k   = kb + 1
1144            wall_index = kb
1145            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1146                                     wall_index, z0(j,i), kb, direction, ncorr )
1147            logc_u_l(k,j,1) = lc
1148            logc_ratio_u_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
1149            lcr(0:ncorr-1) = 1.0_wp
1150
1151!
1152!--         Left boundary for v
1153            i   = -1
1154            kb  =  nzb_v_inner(j,i)
1155            k   =  kb + 1
1156            wall_index = kb
1157            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1158                                     wall_index, z0(j,i), kb, direction, ncorr )
1159            logc_v_l(k,j,1) = lc
1160            logc_ratio_v_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
1161            lcr(0:ncorr-1) = 1.0_wp
1162
1163         ENDDO
1164      ENDIF
1165
1166!
1167!--   Right boundary
1168      IF ( nest_bound_r )  THEN
1169         ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
1170         ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
1171         ALLOCATE( logc_ratio_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
1172         ALLOCATE( logc_ratio_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
1173         logc_u_r       = 0
1174         logc_v_r       = 0
1175         logc_ratio_u_r = 1.0_wp
1176         logc_ratio_v_r = 1.0_wp
1177         direction      = 1
1178         inc            = 1 
1179         DO  j = nys, nyn
1180!
1181!--         Right boundary for u.
1182            i   = nxr + 1
1183            kb  = nzb_u_inner(j,i)
1184            k   = kb + 1
1185            wall_index = kb
1186            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1187                                     wall_index, z0(j,i), kb, direction, ncorr )
1188            logc_u_r(k,j,1) = lc
1189            logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
1190            lcr(0:ncorr-1) = 1.0_wp
1191
1192!
1193!--         Right boundary for v.
1194            i   = nxr + 1
1195            kb  = nzb_v_inner(j,i)
1196            k   = kb + 1
1197            wall_index = kb
1198            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1199                                     wall_index, z0(j,i), kb, direction, ncorr )
1200            logc_v_r(k,j,1) = lc
1201            logc_ratio_v_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
1202            lcr(0:ncorr-1) = 1.0_wp
1203
1204         ENDDO
1205      ENDIF
1206
1207!
1208!--   South boundary
1209      IF ( nest_bound_s )  THEN
1210         ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
1211         ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
1212         ALLOCATE( logc_ratio_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
1213         ALLOCATE( logc_ratio_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
1214         logc_u_s       = 0
1215         logc_v_s       = 0
1216         logc_ratio_u_s = 1.0_wp
1217         logc_ratio_v_s = 1.0_wp
1218         direction      = 1
1219         inc            = 1 
1220         DO  i = nxl, nxr
1221!
1222!--         South boundary for u.
1223            j   = -1
1224            kb  =  nzb_u_inner(j,i)
1225            k   =  kb + 1
1226            wall_index = kb
1227            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1228                                     wall_index, z0(j,i), kb, direction, ncorr )
1229            logc_u_s(k,i,1) = lc
1230            logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
1231            lcr(0:ncorr-1) = 1.0_wp
1232
1233!
1234!--         South boundary for v
1235            j   = 0
1236            kb  = nzb_v_inner(j,i)
1237            k   = kb + 1
1238            wall_index = kb
1239            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1240                                     wall_index, z0(j,i), kb, direction, ncorr )
1241            logc_v_s(k,i,1) = lc
1242            logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
1243            lcr(0:ncorr-1) = 1.0_wp
1244         ENDDO
1245      ENDIF
1246
1247!
1248!--   North boundary
1249      IF ( nest_bound_n )  THEN
1250         ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
1251         ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
1252         ALLOCATE( logc_ratio_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
1253         ALLOCATE( logc_ratio_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
1254         logc_u_n       = 0
1255         logc_v_n       = 0
1256         logc_ratio_u_n = 1.0_wp
1257         logc_ratio_v_n = 1.0_wp 
1258         direction      = 1
1259         inc            = 1 
1260         DO  i = nxl, nxr
1261!
1262!--         North boundary for u.
1263            j   = nyn + 1
1264            kb  = nzb_u_inner(j,i)
1265            k   = kb + 1
1266            wall_index = kb
1267            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1268                                     wall_index, z0(j,i), kb, direction, ncorr )
1269            logc_u_n(k,i,1) = lc
1270            logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
1271            lcr(0:ncorr-1) = 1.0_wp
1272
1273!
1274!--         North boundary for v.
1275            j   = nyn + 1
1276            kb  = nzb_v_inner(j,i)
1277            k   = kb + 1
1278            wall_index = kb
1279            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
1280                                     wall_index, z0(j,i), kb, direction, ncorr )
1281            logc_v_n(k,i,1) = lc
1282            logc_ratio_v_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
1283            lcr(0:ncorr-1) = 1.0_wp
1284
1285         ENDDO
1286      ENDIF
1287
1288!       
1289!--   Then vertical walls and corners if necessary.
1290      IF ( topography /= 'flat' ) THEN       
1291         kb = 0       ! kb is not used when direction > 1
1292!       
1293!--      Left boundary
1294         IF ( nest_bound_l )  THEN
1295            ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
1296            ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
1297            logc_w_l       = 0                   
1298            logc_ratio_w_l = 1.0_wp
1299            direction      = 2
1300            DO  j = nys, nyn
1301               DO  k = nzb, nzt_topo_nestbc_l
1302
1303!
1304!--               Wall for u on the south side, but not on the north side.
1305                  i  = 0
1306                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.     &
1307                       ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )  THEN
1308                     inc        =  1
1309                     wall_index =  j
1310                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1311
1312!
1313!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1314                     logc_u_l(k,j,2) = inc * lc
1315                     logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1316                     lcr(0:ncorr-1) = 1.0_wp
1317                  ENDIF
1318
1319!
1320!--               Wall for u on the north side, but not on the south side.
1321                  i  = 0
1322                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND.     &
1323                       ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
1324                     inc        = -1                 
1325                     wall_index =  j + 1
1326                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1327
1328!
1329!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1330                     logc_u_l(k,j,2) = inc * lc
1331                     logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1332                     lcr(0:ncorr-1) = 1.0_wp
1333                  ENDIF
1334
1335!
1336!--               Wall for w on the south side, but not on the north side.
1337                  i  = -1
1338                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.    &
1339                       ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) )  THEN
1340                     inc        =  1
1341                     wall_index =  j
1342                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1343
1344!
1345!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1346                     logc_w_l(k,j,2) = inc * lc
1347                     logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1348                     lcr(0:ncorr-1) = 1.0_wp
1349                  ENDIF
1350
1351!
1352!--               Wall for w on the north side, but not on the south side.
1353                  i  = -1
1354                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND.     &
1355                       ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
1356                     inc        = -1
1357                     wall_index =  j+1 
1358                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1359
1360!
1361!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1362                     logc_w_l(k,j,2) = inc * lc
1363                     logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1364                     lcr(0:ncorr-1) = 1.0_wp
1365                  ENDIF
1366               ENDDO
1367            ENDDO
1368         ENDIF   !  IF ( nest_bound_l )
1369
1370!       
1371!--      Right boundary.
1372         IF ( nest_bound_r )  THEN
1373            ALLOCATE( logc_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
1374            ALLOCATE( logc_ratio_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
1375            logc_w_r       = 0
1376            logc_ratio_w_r = 1.0_wp
1377            direction      = 2
1378            i  = nxr + 1       
1379            DO  j = nys, nyn
1380               DO  k = nzb, nzt_topo_nestbc_r
1381
1382!
1383!--               Wall for u on the south side, but not on the north side.
1384                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) )  .AND.    &
1385                       ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )  THEN
1386                     inc        = 1
1387                     wall_index = j
1388                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1389
1390!
1391!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1392                     logc_u_r(k,j,2) = inc * lc
1393                     logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1394                     lcr(0:ncorr-1) = 1.0_wp
1395                  ENDIF
1396
1397!
1398!--               Wall for u on the north side, but not on the south side.
1399                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) )  .AND.    &
1400                       ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
1401                     inc        = -1                 
1402                     wall_index =  j+1
1403                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1404
1405!
1406!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1407                     logc_u_r(k,j,2) = inc * lc
1408                     logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1409                     lcr(0:ncorr-1) = 1.0_wp
1410                  ENDIF
1411
1412!
1413!--               Wall for w on the south side, but not on the north side.
1414                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.    &
1415                       ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) )  THEN
1416                     inc        =  1
1417                     wall_index =  j
1418                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1419
1420!
1421!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1422                     logc_w_r(k,j,2) = inc * lc
1423                     logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1424                     lcr(0:ncorr-1) = 1.0_wp
1425                  ENDIF
1426!
1427!--               Wall for w on the north side, but not on the south side.
1428                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND.     &
1429                       ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
1430                     inc        = -1
1431                     wall_index =  j+1 
1432                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
1433
1434!
1435!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1436                     logc_w_r(k,j,2) = inc * lc
1437                     logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
1438                     lcr(0:ncorr-1) = 1.0_wp
1439                  ENDIF
1440               ENDDO
1441            ENDDO
1442         ENDIF   !  IF ( nest_bound_r )
1443
1444!       
1445!--      South boundary.
1446         IF ( nest_bound_s )  THEN
1447            ALLOCATE( logc_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2) )
1448            ALLOCATE( logc_ratio_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2, 0:ncorr-1) )
1449            logc_w_s       = 0
1450            logc_ratio_w_s = 1.0_wp 
1451            direction      = 3
1452            DO  i = nxl, nxr
1453               DO  k = nzb, nzt_topo_nestbc_s
1454
1455!
1456!--               Wall for v on the left side, but not on the right side.
1457                  j  = 0
1458                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.    &
1459                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
1460                     inc        =  1
1461                     wall_index =  i
1462                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1463
1464!
1465!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1466                     logc_v_s(k,i,2) = inc * lc
1467                     logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1468                     lcr(0:ncorr-1) = 1.0_wp
1469                  ENDIF
1470!
1471!--               Wall for v on the right side, but not on the left side.
1472                  j  = 0
1473                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.    &
1474                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
1475                     inc        = -1
1476                     wall_index =  i+1
1477                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1478
1479!
1480!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1481                     logc_v_s(k,i,2) = inc * lc
1482                     logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1483                     lcr(0:ncorr-1) = 1.0_wp
1484                  ENDIF
1485
1486!
1487!--               Wall for w on the left side, but not on the right side.
1488                  j  = -1
1489                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.    &
1490                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
1491                     inc        =  1
1492                     wall_index =  i
1493                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1494
1495!
1496!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1497                     logc_w_s(k,i,2) = inc * lc
1498                     logc_ratio_w_s(k,i,2,0:ncorr - 1) = lcr(0:ncorr-1)
1499                     lcr(0:ncorr-1) = 1.0_wp
1500                  ENDIF
1501
1502!
1503!--               Wall for w on the right side, but not on the left side.
1504                  j  = -1
1505                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.    &
1506                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
1507                     inc        = -1
1508                     wall_index =  i+1 
1509                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1510
1511!
1512!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1513                     logc_w_s(k,i,2) = inc * lc
1514                     logc_ratio_w_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1515                     lcr(0:ncorr-1) = 1.0_wp
1516                  ENDIF
1517               ENDDO
1518            ENDDO
1519         ENDIF   !  IF (nest_bound_s )
1520
1521!       
1522!--      North boundary.
1523         IF ( nest_bound_n )  THEN
1524            ALLOCATE( logc_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2) )
1525            ALLOCATE( logc_ratio_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2, 0:ncorr-1) )
1526            logc_w_n       = 0
1527            logc_ratio_w_n = 1.0_wp
1528            direction      = 3
1529            j  = nyn + 1
1530            DO  i = nxl, nxr
1531               DO  k = nzb, nzt_topo_nestbc_n
1532
1533!
1534!--               Wall for v on the left side, but not on the right side.
1535                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.    &
1536                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
1537                     inc        = 1
1538                     wall_index = i
1539                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1540
1541!
1542!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1543                     logc_v_n(k,i,2) = inc * lc
1544                     logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1545                     lcr(0:ncorr-1) = 1.0_wp
1546                  ENDIF
1547
1548!
1549!--               Wall for v on the right side, but not on the left side.
1550                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.    &
1551                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
1552                     inc        = -1                 
1553                     wall_index =  i + 1
1554                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1555
1556!
1557!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1558                     logc_v_n(k,i,2) = inc * lc
1559                     logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1560                     lcr(0:ncorr-1) = 1.0_wp
1561                  ENDIF
1562
1563!
1564!--               Wall for w on the left side, but not on the right side.
1565                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.    &
1566                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
1567                     inc        = 1
1568                     wall_index = i
1569                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1570
1571!
1572!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1573                     logc_w_n(k,i,2) = inc * lc
1574                     logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1575                     lcr(0:ncorr-1) = 1.0_wp
1576                  ENDIF
1577
1578!
1579!--               Wall for w on the right side, but not on the left side.
1580                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND.     &
1581                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
1582                     inc        = -1
1583                     wall_index =  i+1 
1584                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
1585
1586!
1587!--                  The direction of the wall-normal index is stored as the sign of the logc-element.
1588                     logc_w_n(k,i,2) = inc * lc
1589                     logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
1590                     lcr(0:ncorr-1) = 1.0_wp
1591                  ENDIF
1592               ENDDO
1593            ENDDO
1594         ENDIF   !  IF ( nest_bound_n )
1595      ENDIF   !  IF ( topography /= 'flat' )
1596
1597   END SUBROUTINE pmci_init_loglaw_correction
1598
1599
1600
1601   SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc, wall_index, z0_l, kb, direction, ncorr )
1602      IMPLICIT NONE
1603      INTEGER(iwp), INTENT(IN)  ::  direction                 !:
1604      INTEGER(iwp), INTENT(IN)  ::  ij                        !:
1605      INTEGER(iwp), INTENT(IN)  ::  inc                       !:
1606      INTEGER(iwp), INTENT(IN)  ::  k                         !:
1607      INTEGER(iwp), INTENT(IN)  ::  kb                        !:
1608      INTEGER(iwp), INTENT(OUT) ::  lc                        !:
1609      INTEGER(iwp), INTENT(IN)  ::  ncorr                     !:
1610      INTEGER(iwp), INTENT(IN)  ::  wall_index                !:
1611      REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !:
1612      REAL(wp), INTENT(IN)      ::  z0_l                      !:
1613     
1614      INTEGER(iwp) ::  alcorr       !:
1615      INTEGER(iwp) ::  corr_index   !:
1616      INTEGER(iwp) ::  lcorr        !:
1617      REAL(wp)     ::  logvelc1     !:
1618      LOGICAL      ::  more         !:
1619     
1620
1621      SELECT CASE ( direction )
1622
1623         CASE (1)   !  k
1624            more = .TRUE.
1625            lcorr = 0
1626            DO  WHILE ( more .AND. lcorr <= ncorr - 1 )
1627               corr_index = k + lcorr
1628               IF ( lcorr == 0 )  THEN
1629                  CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
1630               ENDIF
1631               
1632               IF ( corr_index < lc )  THEN
1633                  lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
1634                  more = .TRUE.
1635               ELSE
1636                  lcr(lcorr) = 1.0
1637                  more = .FALSE.
1638               ENDIF
1639               lcorr = lcorr + 1             
1640            ENDDO
1641
1642         CASE (2)   !  j
1643            more = .TRUE.
1644            lcorr  = 0
1645            alcorr = 0
1646            DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
1647               corr_index = ij + lcorr   !  In this case (direction = 2) ij is j
1648               IF ( lcorr == 0 )  THEN
1649                  CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index, z0_l, inc )
1650               ENDIF
1651
1652!
1653!--            The role of inc here is to make the comparison operation "<" valid in both directions.
1654               IF ( inc * corr_index < inc * lc )  THEN
1655                  lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy  &
1656                              - coord_y(wall_index) ) / z0_l ) / logvelc1
1657                  more = .TRUE.
1658               ELSE
1659                  lcr(alcorr) = 1.0_wp
1660                  more = .FALSE.
1661               ENDIF
1662               lcorr  = lcorr + inc
1663               alcorr = ABS( lcorr )
1664            ENDDO
1665
1666         CASE (3)   !  i         
1667            more = .TRUE.
1668            lcorr  = 0
1669            alcorr = 0
1670            DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
1671               corr_index = ij + lcorr   !  In this case (direction = 3) ij is i
1672               IF ( lcorr == 0 )  THEN
1673                  CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index, z0_l, inc )
1674               ENDIF
1675
1676!
1677!--            The role of inc here is to make the comparison operation "<" valid in both directions.
1678               IF ( inc * corr_index < inc * lc )  THEN
1679                  lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx &
1680                       - coord_x(wall_index) ) / z0_l ) / logvelc1
1681                  more = .TRUE.
1682               ELSE
1683                  lcr(alcorr) = 1.0_wp
1684                  more = .FALSE.
1685               ENDIF
1686               lcorr  = lcorr + inc
1687               alcorr = ABS( lcorr )
1688            ENDDO
1689 
1690      END SELECT
1691
1692   END SUBROUTINE pmci_define_loglaw_correction_parameters
1693
1694
1695
1696   SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
1697
1698!
1699!--   Finds the pivot node and te log-law factor for near-wall nodes for
1700!--   which the wall-parallel velocity components will be log-law corrected
1701!--   after interpolation. This subroutine is only for horizontal walls.
1702!
1703!--                               Antti Hellsten 30.12.2015
1704      IMPLICIT NONE
1705      REAL(wp),INTENT(OUT)  ::  logzc1     !:
1706      REAL(wp), INTENT(IN)  ::  z0_l       !:
1707      INTEGER(iwp), INTENT(IN)    ::  kb   !:
1708      INTEGER(iwp), INTENT(OUT)   ::  lc   !:
1709     
1710      REAL(wp)     ::  zuc1   !:
1711      INTEGER(iwp) ::  kbc    !:
1712      INTEGER(iwp) ::  k1     !:
1713
1714
1715      kbc = nzb + 1
1716      DO  WHILE ( cg%zu(kbc) < zu(kb) )   !  kbc is the first coarse-grid point above the surface.
1717         kbc = kbc + 1
1718      ENDDO
1719      zuc1  = cg%zu(kbc)
1720      k1    = kb + 1 
1721      DO  WHILE ( zu(k1) < zuc1 )         !  Important: must be <, not <=
1722         k1 = k1 + 1
1723      ENDDO
1724      logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
1725      lc = k1
1726
1727   END SUBROUTINE pmci_find_logc_pivot_k
1728
1729
1730
1731   SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
1732
1733!
1734!--   Finds the pivot node and te log-law factor for near-wall nodes for
1735!--   which the wall-parallel velocity components will be log-law corrected
1736!--   after interpolation. This subroutine is only for vertical walls on
1737!--   south/north sides of the node.
1738!
1739!--                               Antti Hellsten 5.1.2016
1740      IMPLICIT NONE
1741      REAL(wp), INTENT(IN)      ::  z0_l   !:
1742      INTEGER(iwp), INTENT(IN)  ::  inc    !:  increment must be 1 or -1.
1743      INTEGER(iwp), INTENT(IN)  ::  j      !:
1744      INTEGER(iwp), INTENT(IN)  ::  jw     !:
1745      INTEGER(iwp), INTENT(OUT) ::  lc     !:
1746     
1747      REAL(wp)     ::  logyc1   !:
1748      REAL(wp)     ::  yc1      !:
1749      INTEGER(iwp) ::  j1       !:
1750
1751!
1752!--   yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from the wall.
1753      yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy 
1754
1755!
1756!--   j1 is the first fine-grid index futher away from the wall than yc1.
1757      j1 = j   
1758      DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )   !  Important: must be <, not <=
1759         j1 = j1 + inc
1760      ENDDO
1761
1762      logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
1763      lc = j1
1764
1765   END SUBROUTINE pmci_find_logc_pivot_j
1766
1767
1768
1769   SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
1770
1771!
1772!--   Finds the pivot node and the log-law factor for near-wall nodes for
1773!--   which the wall-parallel velocity components will be log-law corrected
1774!--   after interpolation. This subroutine is only for vertical walls on
1775!--   south/north sides of the node.
1776!
1777!--                               Antti Hellsten 8.1.2016
1778
1779      IMPLICIT NONE
1780      REAL(wp), INTENT(IN)      ::  z0_l   !:
1781      INTEGER(iwp), INTENT(IN)  ::  i      !:
1782      INTEGER(iwp), INTENT(IN)  ::  inc    !: increment must be 1 or -1.
1783      INTEGER(iwp), INTENT(IN)  ::  iw     !:
1784      INTEGER(iwp), INTENT(OUT) ::  lc     !:
1785
1786      REAL(wp)     ::  logxc1   !:
1787      REAL(wp)     ::  xc1      !:
1788      INTEGER(iwp) ::  i1       !:
1789
1790!
1791!--   xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from the wall.
1792      xc1  = coord_x(iw) + 0.5_wp *inc * cg%dx 
1793
1794!
1795!--   i1 is the first fine-grid index futher away from the wall than xc1.
1796      i1 = i   
1797      DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc *xc1 )   !  Important: must be <, not <=
1798         i1 = i1 + inc
1799      ENDDO
1800     
1801      logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
1802      lc = i1
1803
1804   END SUBROUTINE pmci_find_logc_pivot_i
1805
1806
1807
1808   SUBROUTINE pmci_init_anterp_tophat
1809!
1810!--   Precomputation of the client-array indices for
1811!--   corresponding coarse-grid array index and the
1812!--   Under-relaxation coefficients to be used by anterp_tophat.
1813!
1814!--                               Antti Hellsten 9.10.2015.
1815      IMPLICIT NONE
1816      INTEGER(iwp) ::  i        !:
1817      INTEGER(iwp) ::  istart   !:
1818      INTEGER(iwp) ::  j        !:
1819      INTEGER(iwp) ::  jstart   !:
1820      INTEGER(iwp) ::  k        !:
1821      INTEGER(iwp) ::  kstart   !:
1822      INTEGER(iwp) ::  l        !:
1823      INTEGER(iwp) ::  m        !:
1824      INTEGER(iwp) ::  n        !:
1825      REAL(wp)     ::  xi       !:
1826      REAL(wp)     ::  eta      !:
1827      REAL(wp)     ::  zeta     !:
1828     
1829!
1830!--   Default values:
1831      IF ( anterp_relax_length_l < 0.0_wp )  THEN
1832         anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx
1833      ENDIF
1834      IF ( anterp_relax_length_r < 0.0_wp )  THEN
1835         anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx
1836      ENDIF
1837      IF ( anterp_relax_length_s < 0.0_wp )  THEN
1838         anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy
1839      ENDIF
1840      IF ( anterp_relax_length_n < 0.0_wp )  THEN
1841         anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy
1842      ENDIF
1843      IF ( anterp_relax_length_t < 0.0_wp )  THEN
1844         anterp_relax_length_t = 0.1_wp * zu(nzt)
1845      ENDIF
1846
1847!
1848!--   First determine kceu and kcew that are the coarse-grid upper bounds for index k.
1849      n = 0
1850      DO  WHILE ( cg%zu(n) < zu(nzt) )
1851         n = n + 1
1852      ENDDO
1853      kceu = n - 1 
1854
1855      n = 0
1856      DO  WHILE ( cg%zw(n) < zw(nzt-1) )
1857         n = n + 1
1858      ENDDO
1859      kcew = n - 1 
1860
1861      ALLOCATE( iflu(icl:icr) )
1862      ALLOCATE( iflo(icl:icr) )
1863      ALLOCATE( ifuu(icl:icr) )
1864      ALLOCATE( ifuo(icl:icr) )
1865      ALLOCATE( jflv(jcs:jcn) )
1866      ALLOCATE( jflo(jcs:jcn) )
1867      ALLOCATE( jfuv(jcs:jcn) )
1868      ALLOCATE( jfuo(jcs:jcn) )
1869      ALLOCATE( kflw(0:kcew) )
1870      ALLOCATE( kflo(0:kceu) )
1871      ALLOCATE( kfuw(0:kcew) )
1872      ALLOCATE( kfuo(0:kceu) )
1873
1874!
1875!--   i-indices of u for each l-index value.   
1876      istart = nxlg
1877      DO  l = icl, icr 
1878         i = istart
1879         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
1880            i = i + 1
1881         ENDDO
1882         iflu(l) = MIN( MAX( i, nxlg ), nxrg )
1883         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
1884            i = i + 1
1885         ENDDO
1886         ifuu(l) = MIN( MAX( i, nxlg ), nxrg )
1887         istart = iflu(l)
1888      ENDDO
1889
1890!
1891!--   i-indices of others for each l-index value.
1892      istart = nxlg
1893      DO  l = icl, icr   
1894         i = istart
1895         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) )  .AND.  ( i < nxrg ) ) 
1896            i = i + 1
1897         ENDDO
1898         iflo(l) = MIN( MAX( i, nxlg ), nxrg )
1899         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) + cg%dx )  .AND.  ( i < nxrg ) ) 
1900            i = i + 1
1901         ENDDO
1902         ifuo(l) = MIN(MAX(i,nxlg),nxrg)
1903         istart = iflo(l)
1904      ENDDO
1905
1906!
1907!--   j-indices of v for each m-index value.
1908      jstart = nysg
1909      DO  m = jcs, jcn
1910         j = jstart
1911         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
1912            j = j + 1
1913         ENDDO
1914         jflv(m) = MIN( MAX( j, nysg ), nyng ) 
1915         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
1916            j = j + 1
1917         ENDDO
1918         jfuv(m) = MIN( MAX( j, nysg ), nyng )
1919         jstart = jflv(m)
1920      ENDDO
1921
1922!
1923!--   j-indices of others for each m-index value.
1924      jstart = nysg
1925      DO  m = jcs, jcn
1926         j = jstart
1927         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) )  .AND.  ( j < nyng ) ) 
1928            j = j + 1
1929         ENDDO
1930         jflo(m) = MIN( MAX( j, nysg ), nyng )
1931         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) + cg%dy )  .AND.  ( j < nyng ) ) 
1932            j = j + 1
1933         ENDDO
1934         jfuo(m) = MIN( MAX( j, nysg ), nyng )
1935         jstart = jflv(m)
1936      ENDDO
1937
1938!
1939!--   k-indices of w for each n-index value.
1940      kstart  = 0
1941      kflw(0) = 0
1942      kfuw(0) = 0
1943      DO  n = 1, kcew
1944         k = kstart
1945         DO  WHILE ( ( zw(k) < cg%zw(n) - 0.5_wp * cg%dzw(n) )  .AND.  ( k < nzt ) ) 
1946            k = k + 1
1947         ENDDO
1948         kflw(n) = MIN( MAX( k, 1 ), nzt + 1 )
1949         DO  WHILE ( ( zw(k) < cg%zw(n) + 0.5_wp * cg%dzw(n+1) )  .AND.  ( k < nzt ) ) 
1950            k = k + 1
1951         ENDDO
1952         kfuw(n) = MIN( MAX( k, 1 ), nzt + 1 )
1953         kstart = kflw(n)
1954      ENDDO
1955
1956!
1957!--   k-indices of others for each n-index value.
1958      kstart  = 0
1959      kflo(0) = 0
1960      kfuo(0) = 0
1961      DO  n = 1, kceu
1962         k = kstart
1963         DO  WHILE ( ( zu(k) < cg%zu(n) - 0.5_wp * cg%dzu(n) )  .AND.  ( k < nzt ) ) 
1964            k = k + 1
1965         ENDDO
1966         kflo(n) = MIN( MAX( k, 1 ), nzt + 1 )
1967         DO  WHILE ( ( zu(k) < cg%zu(n) + 0.5_wp * cg%dzu(n+1) )  .AND.  ( k < nzt ) ) 
1968            k = k + 1
1969         ENDDO
1970         kfuo(n) = MIN( MAX( k-1, 1 ), nzt + 1 )
1971         kstart = kflo(n)
1972      ENDDO
1973     
1974!
1975!--   Spatial under-relaxation coefficients
1976      ALLOCATE( frax(icl:icr) )
1977      DO  l = icl, icr
1978         IF ( nest_bound_l ) THEN
1979            xi   = ( ( cg%coord_x(l) - lower_left_coord_x ) / anterp_relax_length_l )**4
1980         ELSEIF ( nest_bound_r ) THEN
1981            xi   = ( ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(l) ) / anterp_relax_length_r )**4
1982         ELSE
1983            xi   = 999999.9_wp
1984         ENDIF
1985         frax(l)  = xi / ( 1.0_wp + xi )
1986      ENDDO
1987
1988      ALLOCATE( fray(jcs:jcn) )
1989      DO  m = jcs, jcn
1990         IF ( nest_bound_s ) THEN           
1991            eta  = ( ( cg%coord_y(m) - lower_left_coord_y ) / anterp_relax_length_s )**4
1992         ELSEIF ( nest_bound_n ) THEN
1993            eta  = ( (lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(m)) / anterp_relax_length_n )**4
1994         ELSE
1995            eta  = 999999.9_wp
1996         ENDIF
1997         fray(m)  = eta / ( 1.0_wp + eta )
1998      ENDDO
1999     
2000      ALLOCATE( fraz(0:kceu) )
2001      DO  n = 0, kceu       
2002         zeta = ( ( zu(nzt) - cg%zu(n) ) / anterp_relax_length_t )**4       
2003         fraz(n) = zeta / ( 1.0_wp + zeta )
2004      ENDDO
2005
2006   END SUBROUTINE pmci_init_anterp_tophat
2007
2008
2009
2010   SUBROUTINE pmci_init_tkefactor
2011
2012!
2013!--   Computes the scaling factor for the SGS TKE from coarse grid to be used
2014!--   as BC for the fine grid. Based on the Kolmogorov energy spectrum
2015!--   for the inertial subrange and assumption of sharp cut-off of the resolved
2016!--   energy spectrum. Near the surface, the reduction of TKE is made
2017!--   smaller than further away from the surface.
2018!
2019!                      Antti Hellsten 4.3.2015
2020!
2021!--   Extended for non-flat topography and variable dz.
2022!
2023!                      Antti Hellsten 26.3.2015
2024!
2025!--   The current near-wall adaptation can be replaced by a new one which
2026!--   uses a step function [0,1] based on the logc-arrays. AH 30.12.2015
2027      IMPLICIT NONE
2028      REAL(wp), PARAMETER  ::  cfw = 0.2_wp          !:
2029      REAL(wp), PARAMETER  ::  c_tkef = 0.6_wp       !:
2030      REAL(wp)             ::  fw                    !:
2031      REAL(wp), PARAMETER  ::  fw0 = 0.9_wp          !:
2032      REAL(wp)             ::  glsf                  !:
2033      REAL(wp)             ::  glsc                  !:
2034      REAL(wp)             ::  height                !:
2035      REAL(wp), PARAMETER  ::  p13 = 1.0_wp/3.0_wp   !:
2036      REAL(wp), PARAMETER  ::  p23 = 2.0_wp/3.0_wp   !:
2037      INTEGER(iwp)         ::  k                     !:
2038      INTEGER(iwp)         ::  kc                    !:
2039       
2040
2041      IF ( nest_bound_l )  THEN
2042         ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )   
2043         tkefactor_l = 0.0_wp
2044         i = nxl - 1 
2045         DO  j = nysg, nyng
2046            DO  k = nzb_s_inner(j,i) + 1, nzt
2047               kc     = kco(k+1)
2048               glsf   = ( dx * dy * dzu(k) )**p13
2049               glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
2050               height = zu(k) - zu(nzb_s_inner(j,i))
2051               fw     = EXP( -cfw * height / glsf )
2052               tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
2053            ENDDO
2054            tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0
2055         ENDDO
2056      ENDIF
2057
2058      IF ( nest_bound_r )  THEN
2059         ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )   
2060         tkefactor_r = 0.0_wp
2061         i = nxr + 1
2062         DO  j = nysg, nyng
2063            DO  k = nzb_s_inner(j,i) + 1, nzt
2064               kc     = kco(k+1)
2065               glsf   = ( dx * dy * dzu(k) )**p13
2066               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
2067               height = zu(k) - zu(nzb_s_inner(j,i))
2068               fw     = EXP( -cfw * height / glsf )
2069               tkefactor_r(k,j) = c_tkef * (fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
2070            ENDDO
2071            tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0
2072         ENDDO
2073      ENDIF
2074
2075      IF ( nest_bound_s )  THEN
2076         ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
2077         tkefactor_s = 0.0_wp
2078         j = nys - 1
2079         DO  i = nxlg, nxrg
2080            DO  k = nzb_s_inner(j,i) + 1, nzt
2081               kc     = kco(k+1)
2082               glsf   = ( dx * dy * dzu(k) )**p13
2083               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
2084               height = zu(k) - zu(nzb_s_inner(j,i))
2085               fw     = EXP( -cfw*height / glsf )
2086               tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
2087            ENDDO
2088            tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0
2089         ENDDO
2090      ENDIF
2091
2092      IF ( nest_bound_n )  THEN
2093         ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
2094         tkefactor_n = 0.0_wp   
2095         j = nyn + 1
2096         DO  i = nxlg, nxrg
2097            DO  k = nzb_s_inner(j,i)+1, nzt
2098               kc     = kco(k+1)
2099               glsf   = ( dx * dy * dzu(k) )**p13
2100               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
2101               height = zu(k) - zu(nzb_s_inner(j,i))
2102               fw     = EXP( -cfw * height / glsf )
2103               tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
2104            ENDDO
2105            tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0
2106         ENDDO
2107      ENDIF
2108
2109      ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
2110      k = nzt
2111      DO  i = nxlg, nxrg
2112         DO  j = nysg, nyng
2113            kc     = kco(k+1)
2114            glsf   = ( dx * dy * dzu(k) )**p13
2115            glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
2116            height = zu(k) - zu(nzb_s_inner(j,i))
2117            fw     = EXP( -cfw * height / glsf )
2118            tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
2119         ENDDO
2120      ENDDO
2121     
2122    END SUBROUTINE pmci_init_tkefactor
2123
2124#endif
2125 END SUBROUTINE pmci_setup_client
2126
2127
2128
2129 SUBROUTINE pmci_setup_coordinates
2130
2131#if defined( __parallel )
2132    IMPLICIT NONE
2133
2134    INTEGER(iwp) ::  i   !:
2135    INTEGER(iwp) ::  j   !:
2136
2137!
2138!-- Create coordinate arrays.
2139    ALLOCATE( coord_x(-nbgp:nx+nbgp) )
2140    ALLOCATE( coord_y(-nbgp:ny+nbgp) )
2141     
2142    DO  i = -nbgp, nx + nbgp
2143       coord_x(i) = lower_left_coord_x + i * dx
2144    ENDDO
2145     
2146    DO  j = -nbgp, ny + nbgp
2147       coord_y(j) = lower_left_coord_y + j * dy
2148    ENDDO
2149
2150#endif
2151 END SUBROUTINE pmci_setup_coordinates
2152
2153
2154
2155 SUBROUTINE pmci_server_synchronize
2156
2157#if defined( __parallel )
2158!
2159!-- Unify the time steps for each model and synchronize.
2160!-- This is based on the assumption that the native time step
2161!-- (original dt_3d) of any server is always larger than the smallest
2162!-- native time step of it s clients.
2163    IMPLICIT NONE
2164    INTEGER(iwp)           ::  client_id   !:
2165    REAL(wp), DIMENSION(1) ::  dtc         !:
2166    REAL(wp), DIMENSION(1) ::  dtl         !:
2167    INTEGER(iwp)           ::  ierr        !:
2168    INTEGER(iwp)           ::  m           !:   
2169
2170!
2171!-- First find the smallest native time step of all the clients of the current server.
2172    dtl(1) = 999999.9_wp
2173    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
2174       client_id = PMC_Server_for_Client(m)
2175       IF ( myid == 0 )  THEN
2176          CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )
2177          dtl(1) = MIN( dtl(1), dtc(1) )
2178          dt_3d   = dtl(1)
2179       ENDIF
2180    ENDDO
2181
2182!
2183!-- Broadcast the unified time step to all server processes.
2184    CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
2185
2186!
2187!-- Send the new time step to all the clients of the current server.
2188    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
2189       client_id = PMC_Server_for_Client(m)
2190       IF ( myid == 0 )  THEN
2191          CALL pmc_send_to_client( client_id, dtl, SIZE( dtl ), 0, 102, ierr )
2192       ENDIF
2193    ENDDO
2194
2195#endif
2196 END SUBROUTINE pmci_server_synchronize
2197
2198
2199
2200 SUBROUTINE pmci_client_synchronize
2201
2202#if defined( __parallel )
2203!
2204!-- Unify the time steps for each model and synchronize.
2205!-- This is based on the assumption that the native time step
2206!-- (original dt_3d) of any server is always larger than the smallest
2207!-- native time step of it s clients.
2208
2209    IMPLICIT NONE
2210    REAL(wp), DIMENSION(1) ::  dtl    !:
2211    REAL(wp), DIMENSION(1) ::  dts    !:
2212    INTEGER(iwp)           ::  ierr   !:
2213   
2214
2215    dtl(1) = dt_3d
2216    IF ( cpl_id > 1 )  THEN   !  Root id is never a client
2217       IF ( myid==0 )  THEN
2218          CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
2219          CALL pmc_recv_from_server( dts, SIZE( dts ), 0, 102, ierr )
2220          dt_3d = dts(1)
2221       ENDIF
2222
2223!
2224!--    Broadcast the unified time step to all server processes.
2225       CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
2226    ENDIF
2227
2228#endif
2229 END SUBROUTINE pmci_client_synchronize
2230               
2231
2232
2233 SUBROUTINE pmci_set_swaplevel( swaplevel )
2234
2235    IMPLICIT NONE
2236
2237    INTEGER(iwp),INTENT(IN) ::  swaplevel  !: swaplevel (1 or 2) of PALM's timestep
2238
2239    INTEGER(iwp)            ::  client_id  !:
2240    INTEGER(iwp)            ::  m          !:
2241
2242!
2243!-- After each timestep, alternately set buffer one or buffer two active
2244    DO  m = 1, SIZE( pmc_server_for_client )-1
2245       client_id = pmc_server_for_client(m)
2246       CALL pmc_s_set_active_data_array( client_id, swaplevel )
2247    ENDDO
2248
2249 END SUBROUTINE pmci_set_swaplevel
2250
2251
2252
2253 SUBROUTINE pmci_server_datatrans( direction )
2254
2255    IMPLICIT NONE
2256
2257    INTEGER(iwp),INTENT(IN) ::  direction   !:
2258
2259#if defined( __parallel )
2260    INTEGER(iwp)            ::  client_id   !:
2261    INTEGER(iwp)            ::  i           !:
2262    INTEGER(iwp)            ::  j           !:
2263    INTEGER(iwp)            ::  ierr        !:
2264    INTEGER(iwp)            ::  m           !:
2265    REAL(wp)                ::  waittime    !:
2266    REAL(wp), DIMENSION(1)  ::  dtc         !:
2267    REAL(wp), DIMENSION(1)  ::  dtl         !:
2268
2269!
2270!-- First find the smallest native time step of all the clients of the current server.
2271    dtl(1) = 999999.9_wp
2272    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
2273       client_id = PMC_Server_for_Client(m)
2274       IF ( myid==0 )  THEN
2275          CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )
2276          dtl(1) = MIN( dtl(1), dtc(1) )
2277          dt_3d  = dtl(1)
2278       ENDIF
2279    ENDDO
2280
2281!
2282!-- Broadcast the unified time step to all server processes.
2283    CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
2284
2285    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
2286       client_id = PMC_Server_for_Client(m)
2287       CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
2288
2289!
2290!--    Send the new time step to all the clients of the current server.
2291       IF ( myid == 0 ) THEN
2292          CALL pmc_send_to_client( client_id, dtl, SIZE( dtl ), 0, 102, ierr )
2293       ENDIF
2294       CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
2295       
2296       IF ( direction == server_to_client )  THEN
2297          CALL cpu_log( log_point_s(71), 'PMC Server Send', 'start' )
2298          CALL pmc_s_fillbuffer( client_id, waittime=waittime )
2299          CALL cpu_log( log_point_s(71), 'PMC Server Send', 'stop' )
2300       ELSE   !  Communication from client to server.
2301          CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'start' )
2302          client_id = pmc_server_for_client(m)
2303          CALL pmc_s_getdata_from_buffer( client_id )
2304          CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'stop' )       
2305
2306!
2307!--       The anterpolated data is now available in u etc.
2308          IF ( topography /= 'flat' )  THEN
2309
2310!
2311!--          Inside buildings/topography reset velocities and TKE back to zero.
2312!--          Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
2313!--          maybe revise later.
2314             DO   i = nxlg, nxrg
2315                DO   j = nysg, nyng
2316                   u(nzb:nzb_u_inner(j,i),j,i)     = 0.0_wp
2317                   v(nzb:nzb_v_inner(j,i),j,i)     = 0.0_wp
2318                   w(nzb:nzb_w_inner(j,i),j,i)     = 0.0_wp
2319                   e(nzb:nzb_w_inner(j,i),j,i)     = 0.0_wp
2320                   tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
2321                   tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
2322                   tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
2323                   te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
2324                ENDDO
2325             ENDDO
2326          ENDIF
2327       ENDIF
2328    ENDDO
2329
2330#endif
2331 END SUBROUTINE pmci_server_datatrans
2332
2333
2334
2335 SUBROUTINE pmci_client_datatrans( direction )
2336
2337    IMPLICIT NONE
2338
2339    INTEGER(iwp), INTENT(IN) ::  direction   !:
2340
2341#if defined( __parallel )
2342    INTEGER(iwp)             ::  ierr        !:
2343    INTEGER(iwp)             ::  icl         !:
2344    INTEGER(iwp)             ::  icr         !:
2345    INTEGER(iwp)             ::  jcs         !:
2346    INTEGER(iwp)             ::  jcn         !:
2347   
2348    REAL(wp), DIMENSION(1)   ::  dtl         !:
2349    REAL(wp), DIMENSION(1)   ::  dts         !:
2350    REAL(wp)                 ::  waittime    !:
2351
2352
2353    dtl = dt_3d
2354    IF ( cpl_id > 1 )  THEN   !  Root id is never a client     
2355       CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
2356       IF ( myid==0 )  THEN
2357          CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
2358          CALL pmc_recv_from_server( dts, SIZE( dts ), 0, 102, ierr )
2359          dt_3d = dts(1)
2360       ENDIF
2361
2362!
2363!--    Broadcast the unified time step to all server processes.
2364       CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
2365       CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
2366
2367!
2368!--    Client domain boundaries in the server indice space.
2369       icl = coarse_bound(1)
2370       icr = coarse_bound(2)
2371       jcs = coarse_bound(3)
2372       jcn = coarse_bound(4)
2373
2374       IF ( direction == server_to_client )  THEN
2375          CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'start' )
2376          CALL pmc_c_getbuffer( WaitTime = WaitTime )
2377          CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'stop' )
2378
2379!
2380!--       The interpolation. Add IF-condition here: IF not vertical nesting
2381          IF ( nest_bound_l )  THEN   !  Left border pe
2382             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
2383                                       nzb_u_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'u' )
2384             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
2385                                       nzb_v_inner, logc_v_l, logc_ratio_v_l, nzt_topo_nestbc_l, 'l', 'v' )
2386             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
2387                                       nzb_w_inner, logc_w_l, logc_ratio_w_l, nzt_topo_nestbc_l, 'l', 'w' )
2388             CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2389                                       nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'e' )
2390             CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2391                                       nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
2392             IF ( humidity  .OR.  passive_scalar )  THEN
2393                CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
2394                                          nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
2395             ENDIF
2396             IF ( nesting_mode == 'one-way' )  THEN
2397                CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
2398                CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
2399                CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
2400                CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
2401                CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
2402                IF ( humidity  .OR.  passive_scalar )  THEN
2403                   CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
2404                ENDIF
2405             ENDIF
2406          ENDIF
2407          IF ( nest_bound_r )  THEN   !  Right border pe
2408             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
2409                                       nzb_u_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'u' )
2410             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
2411                                       nzb_v_inner, logc_v_r, logc_ratio_v_r, nzt_topo_nestbc_r, 'r', 'v' )
2412             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
2413                                       nzb_w_inner, logc_w_r, logc_ratio_w_r, nzt_topo_nestbc_r, 'r', 'w' )
2414             CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2415                                       nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'e' )
2416             CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2417                                       nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
2418             IF ( humidity  .OR.  passive_scalar )  THEN
2419                CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
2420                                          nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
2421             ENDIF
2422             IF ( nesting_mode == 'one-way' )  THEN
2423                CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
2424                CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'r', 'v' )
2425                CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'r', 'w' )
2426                CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'r', 'e' )
2427                CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'r', 's' )
2428                IF ( humidity  .OR.  passive_scalar )  THEN
2429                   CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'r', 's' )
2430                ENDIF
2431             ENDIF
2432          ENDIF
2433          IF ( nest_bound_s )  THEN   !  South border pe
2434             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
2435                                       nzb_u_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'u' )
2436             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
2437                                       nzb_v_inner, logc_v_s, logc_ratio_v_s, nzt_topo_nestbc_s, 's', 'v' )
2438             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
2439                                       nzb_w_inner, logc_w_s, logc_ratio_w_s, nzt_topo_nestbc_s, 's', 'w' )
2440             CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2441                                       nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'e' )
2442             CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2443                                       nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
2444             IF ( humidity  .OR.  passive_scalar )  THEN
2445                CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
2446                                          nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
2447             ENDIF
2448             IF ( nesting_mode == 'one-way' )  THEN
2449                CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 's', 'u' )
2450                CALL pmci_extrap_ifoutflow_sn( v, nzb_v_inner, 's', 'v' )
2451                CALL pmci_extrap_ifoutflow_sn( w, nzb_w_inner, 's', 'w' )
2452                CALL pmci_extrap_ifoutflow_sn( e, nzb_s_inner, 's', 'e' )
2453                CALL pmci_extrap_ifoutflow_sn( pt,nzb_s_inner, 's', 's' )
2454                IF ( humidity  .OR.  passive_scalar )  THEN
2455                   CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 's', 's' )
2456                ENDIF
2457             ENDIF
2458          ENDIF
2459          IF ( nest_bound_n )  THEN   !  North border pe
2460             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
2461                                       nzb_u_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'u' )
2462             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
2463                                       nzb_v_inner, logc_v_n, logc_ratio_v_n, nzt_topo_nestbc_n, 'n', 'v' )
2464             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
2465                                       nzb_w_inner, logc_w_n, logc_ratio_w_n, nzt_topo_nestbc_n, 'n', 'w' )
2466             CALL pmci_interp_tril_sn( e,  ec,  ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                  &
2467                                       nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'e' )
2468             CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
2469                                       nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
2470             IF ( humidity  .OR.  passive_scalar )  THEN
2471                CALL pmci_interp_tril_sn( q, qc, ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                 &
2472                                          nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
2473             ENDIF
2474             IF ( nesting_mode == 'one-way' )  THEN             
2475                CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 'n', 'u' )
2476                CALL pmci_extrap_ifoutflow_sn( v, nzb_v_inner, 'n', 'v' )
2477                CALL pmci_extrap_ifoutflow_sn( w, nzb_w_inner, 'n', 'w' )
2478                CALL pmci_extrap_ifoutflow_sn( e, nzb_s_inner, 'n', 'e' )
2479                CALL pmci_extrap_ifoutflow_sn( pt,nzb_s_inner, 'n', 's' )
2480                IF ( humidity  .OR.  passive_scalar )  THEN
2481                   CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 'n', 's' )
2482                ENDIF
2483             ENDIF
2484          ENDIF
2485
2486!
2487!--       All PEs are top-border PEs
2488          CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, 'u' )
2489          CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, 'v' )
2490          CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, 'w' )
2491          CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 'e' )
2492          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
2493          IF ( humidity  .OR.  passive_scalar )  THEN
2494             CALL pmci_interp_tril_t( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
2495          ENDIF
2496          IF ( nesting_mode == 'one-way' )  THEN
2497             CALL pmci_extrap_ifoutflow_t( u,  'u' )
2498             CALL pmci_extrap_ifoutflow_t( v,  'v' )
2499             CALL pmci_extrap_ifoutflow_t( w,  'w' )
2500             CALL pmci_extrap_ifoutflow_t( e,  'e' )
2501             CALL pmci_extrap_ifoutflow_t( pt, 's' )
2502             IF ( humidity  .OR.  passive_scalar )  THEN
2503                CALL pmci_extrap_ifoutflow_t( q, 's' )
2504             ENDIF
2505          ENDIF
2506       ELSE    !  IF ( direction == server_to_client )
2507          CALL pmci_anterp_tophat( u,  uc,  kceu, iflu, ifuu, jflo, jfuo, kflo, kfuo, nzb_u_inner, 'u' )
2508          CALL pmci_anterp_tophat( v,  vc,  kceu, iflo, ifuo, jflv, jfuv, kflo, kfuo, nzb_v_inner, 'v' )
2509          CALL pmci_anterp_tophat( w,  wc,  kcew, iflo, ifuo, jflo, jfuo, kflw, kfuw, nzb_w_inner, 'w' )
2510          CALL pmci_anterp_tophat( pt, ptc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
2511          IF ( humidity  .OR.  passive_scalar )  THEN
2512             CALL pmci_anterp_tophat( q, qc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
2513          ENDIF
2514          CALL cpu_log( log_point_s(74), 'PMC Client Send', 'start' )
2515          CALL pmc_c_putbuffer( WaitTime = WaitTime )
2516          CALL cpu_log( log_point_s(74), 'PMC Client Send', 'stop' )
2517       ENDIF
2518    ENDIF
2519   
2520 CONTAINS
2521
2522
2523
2524   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, &
2525                                  nzt_topo_nestbc, edge, var )
2526
2527!
2528!--   Interpolation of ghost-node values used as the client-domain boundary
2529!--   conditions. This subroutine handles the left and right boundaries.
2530!--   This subroutine is based on trilinear interpolation.
2531!--   Constant dz is still assumed.
2532!
2533!--                                   Antti Hellsten 22.2.2015.
2534!
2535!--   Rewritten so that all the coefficients and client-array indices are
2536!--   precomputed in the initialization phase by pmci_init_interp_tril.
2537!
2538!                                     Antti Hellsten 3.3.2015.
2539!
2540!--   Constant dz no more assumed.
2541!                                     Antti Hellsten 23.3.2015.
2542!
2543!--   Adapted for non-flat topography. However, the near-wall velocities
2544!--   are log-corrected only over horizontal surfaces, not yet near vertical
2545!--   walls.
2546!--                                   Antti Hellsten 26.3.2015.
2547!
2548!--   Indexing in the principal direction (i) is changed. Now, the nest-boundary
2549!--   values are interpolated only into the first ghost-node layers on each later
2550!--   boundary. These values are then simply copied to the second ghost-node layer.
2551!
2552!--                                   Antti Hellsten 6.10.2015.
2553      IMPLICIT NONE
2554      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT)          ::  f            !:
2555      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)                 ::  fc           !:
2556      REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2,0:ncorr-1), INTENT(IN) ::  logc_ratio   !:
2557      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r1x          !:
2558      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r2x          !:
2559      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r1y          !:
2560      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r2y          !:
2561      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r1z          !:
2562      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r2z          !:
2563     
2564      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                             ::  ic           !:
2565      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                             ::  jc           !:
2566      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)                   ::  kb           !:
2567      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                             ::  kc           !:
2568      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2), INTENT(IN)       ::  logc         !:
2569      INTEGER(iwp)                                                               ::  nzt_topo_nestbc   !:
2570
2571      CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
2572      CHARACTER(LEN=1),INTENT(IN) ::  var    !:
2573
2574      INTEGER(iwp) ::  i       !:
2575      INTEGER(iwp) ::  ib      !:
2576      INTEGER(iwp) ::  iw      !:
2577      INTEGER(iwp) ::  j       !:
2578      INTEGER(iwp) ::  jco     !:
2579      INTEGER(iwp) ::  jcorr   !:
2580      INTEGER(iwp) ::  jinc    !:
2581      INTEGER(iwp) ::  jw      !:
2582      INTEGER(iwp) ::  j1      !:
2583      INTEGER(iwp) ::  k       !:
2584      INTEGER(iwp) ::  kco     !:
2585      INTEGER(iwp) ::  kcorr   !:
2586      INTEGER(iwp) ::  k1      !:
2587      INTEGER(iwp) ::  l       !:
2588      INTEGER(iwp) ::  m       !:
2589      INTEGER(iwp) ::  n       !:
2590      INTEGER(iwp) ::  kbc     !:
2591     
2592      REAL(wp) ::  coarse_dx   !:
2593      REAL(wp) ::  coarse_dy   !:
2594      REAL(wp) ::  coarse_dz   !:
2595      REAL(wp) ::  fkj         !:
2596      REAL(wp) ::  fkjp        !:
2597      REAL(wp) ::  fkpj        !:
2598      REAL(wp) ::  fkpjp       !:
2599      REAL(wp) ::  fk          !:
2600      REAL(wp) ::  fkp         !:
2601     
2602!
2603!--   Check which edge is to be handled: left or right. Note the assumption that the same PE never
2604!--   holds both left and right nest boundaries. Should this be changed?
2605      IF ( edge == 'l' )  THEN
2606         IF ( var == 'u' )  THEN   !  For u, nxl is a ghost node, but not for the other variables.
2607            i  = nxl
2608            ib = nxl - 1 
2609         ELSE
2610            = nxl - 1
2611            ib = nxl - 2
2612         ENDIF
2613      ELSEIF ( edge == 'r' )  THEN
2614         = nxr + 1
2615         ib = nxr + 2
2616      ENDIF
2617     
2618      DO  j = nys, nyn + 1
2619         DO  k = kb(j,i), nzt + 1
2620            l = ic(i)
2621            m = jc(j)
2622            n = kc(k)
2623            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
2624            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
2625            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
2626            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
2627            fk       = r1y(j) * fkj  + r2y(j) * fkjp
2628            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
2629            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
2630         ENDDO
2631      ENDDO
2632
2633!
2634!--   Generalized log-law-correction algorithm.
2635!--   Doubly two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays
2636!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction.
2637!
2638!--   Solid surface below the node
2639      IF ( var == 'u' .OR. var == 'v' )  THEN           
2640         DO  j = nys, nyn
2641            k = kb(j,i) + 1
2642            IF ( ( logc(k,j,1) /= 0 )  .AND.  ( logc(k,j,2) == 0 ) )  THEN
2643               k1 = logc(k,j,1)
2644               DO  kcorr=0,ncorr - 1
2645                  kco = k + kcorr
2646                  f(kco,j,i) = logc_ratio(k,j,1,kcorr) * f(k1,j,i)
2647               ENDDO
2648            ENDIF
2649         ENDDO
2650      ENDIF
2651
2652!
2653!--   In case of non-flat topography, also vertical walls and corners need to be treated.
2654!--   Only single and double wall nodes are corrected. Triple and higher-multiple wall nodes
2655!--   are not corrected as the log law would not be valid anyway in such locations.
2656      IF ( topography /= 'flat' )  THEN
2657         IF ( var == 'u' .OR. var == 'w' )  THEN                 
2658
2659!
2660!--         Solid surface only on south/north side of the node                   
2661            DO  j = nys, nyn
2662               DO  k = kb(j,i) + 1, nzt_topo_nestbc
2663                  IF ( ( logc(k,j,2) /= 0 )  .AND.  ( logc(k,j,1) == 0 ) )  THEN
2664
2665!
2666!--                  Direction of the wall-normal index is carried in as the sign of logc.           
2667                     jinc = SIGN( 1, logc(k,j,2) )
2668                     j1   = ABS( logc(k,j,2) )
2669                     DO  jcorr=0, ncorr - 1
2670                        jco = j + jinc * jcorr
2671                        f(k,jco,i) = logc_ratio(k,j,2,jcorr) * f(k,j1,i)
2672                     ENDDO
2673                  ENDIF
2674               ENDDO
2675            ENDDO
2676         ENDIF
2677
2678!
2679!--      Solid surface on both below and on south/north side of the node           
2680         IF ( var == 'u' )  THEN
2681            DO  j = nys, nyn
2682               k = kb(j,i) + 1
2683               IF ( ( logc(k,j,2) /= 0 )  .AND.  ( logc(k,j,1) /= 0 ) )  THEN
2684                  k1   = logc(k,j,1)                 
2685                  jinc = SIGN( 1, logc(k,j,2) )
2686                  j1   = ABS( logc(k,j,2) )                 
2687                  DO  jcorr = 0, ncorr - 1
2688                     jco = j + jinc * jcorr
2689                     DO  kcorr = 0, ncorr - 1
2690                        kco = k + kcorr
2691                        f(kco,jco,i) = 0.5_wp * ( logc_ratio(k,j,1,kcorr) * f(k1,j,i)  &
2692                                     + logc_ratio(k,j,2,jcorr) * f(k,j1,i) )
2693                     ENDDO
2694                  ENDDO
2695               ENDIF
2696            ENDDO
2697         ENDIF
2698
2699      ENDIF  ! ( topography /= 'flat' )
2700
2701!
2702!--   Rescale if f is the TKE.
2703      IF ( var == 'e')  THEN
2704         IF ( edge == 'l' )  THEN
2705            DO  j = nys, nyn + 1
2706               DO  k = kb(j,i), nzt + 1
2707                  f(k,j,i) = tkefactor_l(k,j) * f(k,j,i)
2708               ENDDO
2709            ENDDO
2710         ELSEIF ( edge == 'r' )  THEN           
2711            DO  j = nys, nyn + 1
2712               DO  k = kb(j,i), nzt + 1
2713                  f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
2714               ENDDO
2715            ENDDO
2716         ENDIF
2717      ENDIF
2718
2719!
2720!--   Store the boundary values also into the second ghost node layer.
2721      f(0:nzt+1,nys:nyn+1,ib) = f(0:nzt+1,nys:nyn+1,i)
2722
2723   END SUBROUTINE pmci_interp_tril_lr
2724
2725
2726
2727   SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, &
2728                                  nzt_topo_nestbc, edge, var )
2729
2730!
2731!--   Interpolation of ghost-node values used as the client-domain boundary
2732!--   conditions. This subroutine handles the south and north boundaries.
2733!--   This subroutine is based on trilinear interpolation.
2734!--   Constant dz is still assumed.
2735!
2736!--                                      Antti Hellsten 22.2.2015.
2737!
2738!--   Rewritten so that all the coefficients and client-array indices are
2739!--   precomputed in the initialization phase by pmci_init_interp_tril.
2740!
2741!--                                      Antti Hellsten 3.3.2015.
2742!
2743!--   Constant dz no more assumed.
2744!--                                      Antti Hellsten 23.3.2015.
2745!
2746!--   Adapted for non-flat topography. However, the near-wall velocities
2747!--   are log-corrected only over horifontal surfaces, not yet near vertical
2748!--   walls.
2749!--                                      Antti Hellsten 26.3.2015.
2750!
2751!--   Indexing in the principal direction (j) is changed. Now, the nest-boundary
2752!--   values are interpolated only into the first ghost-node layers on each later
2753!--   boundary. These values are then simply copied to the second ghost-node layer.
2754!
2755!--                                      Antti Hellsten 6.10.2015.
2756      IMPLICIT NONE
2757      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT)          ::  f                 !:
2758      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)                 ::  fc                !:
2759      REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2,0:ncorr-1), INTENT(IN) ::  logc_ratio        !:
2760      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r1x               !:
2761      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r2x               !:
2762      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r1y               !:
2763      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r2y               !:
2764      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r1z               !:
2765      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r2z               !:
2766     
2767      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                             ::  ic                !:
2768      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                             ::  jc                !:
2769      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)                   ::  kb                !:
2770      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                             ::  kc                !:
2771      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2), INTENT(IN)       ::  logc              !:
2772      INTEGER(iwp)                                                               ::  nzt_topo_nestbc   !:
2773
2774      CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
2775      CHARACTER(LEN=1), INTENT(IN) ::  var    !:
2776     
2777      INTEGER(iwp) ::  i       !:
2778      INTEGER(iwp) ::  iinc    !:
2779      INTEGER(iwp) ::  icorr   !:
2780      INTEGER(iwp) ::  ico     !:
2781      INTEGER(iwp) ::  i1      !:
2782      INTEGER(iwp) ::  j       !:
2783      INTEGER(iwp) ::  jb      !:
2784      INTEGER(iwp) ::  k       !:
2785      INTEGER(iwp) ::  kcorr   !:
2786      INTEGER(iwp) ::  kco     !:
2787      INTEGER(iwp) ::  k1      !:
2788      INTEGER(iwp) ::  l       !:
2789      INTEGER(iwp) ::  m       !:
2790      INTEGER(iwp) ::  n       !:
2791                           
2792      REAL(wp) ::  coarse_dx   !:
2793      REAL(wp) ::  coarse_dy   !:
2794      REAL(wp) ::  coarse_dz   !:
2795      REAL(wp) ::  fk          !:
2796      REAL(wp) ::  fkj         !:
2797      REAL(wp) ::  fkjp        !:
2798      REAL(wp) ::  fkpj        !:
2799      REAL(wp) ::  fkpjp       !:
2800      REAL(wp) ::  fkp         !:
2801     
2802!
2803!--   Check which edge is to be handled: south or north. Note the assumption that the same PE never
2804!--   holds both south and north nest boundaries. Should this be changed?
2805      IF ( edge == 's' )  THEN
2806         IF ( var == 'v' )  THEN   !  For v, nys is a ghost node, but not for the other variables.
2807            j  = nys
2808            jb = nys - 1 
2809         ELSE
2810            = nys - 1
2811            jb = nys - 2
2812         ENDIF
2813      ELSEIF ( edge == 'n' )  THEN
2814         = nyn + 1
2815         jb = nyn + 2
2816      ENDIF
2817
2818      DO  i = nxl, nxr + 1
2819         DO  k = kb(j,i), nzt + 1
2820            l = ic(i)
2821            m = jc(j)
2822            n = kc(k)             
2823            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
2824            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
2825            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
2826            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
2827            fk       = r1y(j) * fkj  + r2y(j) * fkjp
2828            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
2829            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
2830         ENDDO
2831      ENDDO
2832
2833!
2834!--   Generalized log-law-correction algorithm.
2835!--   Multiply two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays
2836!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction.
2837!
2838!--   Solid surface below the node
2839      IF ( var == 'u'  .OR.  var == 'v' )  THEN           
2840         DO  i = nxl, nxr
2841            k = kb(j,i) + 1
2842            IF ( ( logc(k,i,1) /= 0 )  .AND.  ( logc(k,i,2) == 0 ) )  THEN
2843               k1 = logc(k,i,1)
2844               DO  kcorr = 0, ncorr - 1
2845                  kco = k + kcorr
2846                  f(kco,j,i) = logc_ratio(k,i,1,kcorr) * f(k1,j,i)
2847               ENDDO
2848            ENDIF
2849         ENDDO
2850      ENDIF
2851
2852!
2853!--   In case of non-flat topography, also vertical walls and corners need to be treated.
2854!--   Only single and double wall nodes are corrected.
2855!--   Triple and higher-multiple wall nodes are not corrected as it would be extremely complicated
2856!--   and the log law would not be valid anyway in such locations.
2857      IF ( topography /= 'flat' )  THEN
2858         IF ( var == 'v' .OR. var == 'w' )  THEN
2859            DO  i = nxl, nxr
2860               DO  k = kb(j,i), nzt_topo_nestbc
2861
2862!
2863!--               Solid surface only on left/right side of the node           
2864                  IF ( ( logc(k,i,2) /= 0 )  .AND.  ( logc(k,i,1) == 0 ) )  THEN
2865
2866!
2867!--                  Direction of the wall-normal index is carried in as the sign of logc.           
2868                     iinc = SIGN( 1, logc(k,i,2) )
2869                     i1  = ABS( logc(k,i,2) )
2870                     DO  icorr = 0, ncorr - 1
2871                        ico = i + iinc * icorr
2872                        f(k,j,ico) = logc_ratio(k,i,2,icorr) * f(k,j,i1)
2873                     ENDDO
2874                  ENDIF
2875               ENDDO
2876            ENDDO
2877         ENDIF
2878
2879!
2880!--      Solid surface on both below and on left/right side of the node           
2881         IF ( var == 'v' )  THEN
2882            DO  i = nxl, nxr
2883               k = kb(j,i) + 1
2884               IF ( ( logc(k,i,2) /= 0 )  .AND.  ( logc(k,i,1) /= 0 ) )  THEN
2885                  k1   = logc(k,i,1)         
2886                  iinc = SIGN( 1, logc(k,i,2) )
2887                  i1   = ABS( logc(k,i,2) )
2888                  DO  icorr = 0, ncorr - 1
2889                     ico = i + iinc * icorr
2890                     DO  kcorr = 0, ncorr - 1
2891                        kco = k + kcorr
2892                        f(kco,i,ico) = 0.5_wp * ( logc_ratio(k,i,1,kcorr) * f(k1,j,i)  &
2893                                     + logc_ratio(k,i,2,icorr) * f(k,j,i1) )
2894                     ENDDO
2895                  ENDDO
2896               ENDIF
2897            ENDDO
2898         ENDIF
2899         
2900      ENDIF  ! ( topography /= 'flat' )
2901
2902!
2903!--   Rescale if f is the TKE.
2904      IF ( var == 'e')  THEN
2905         IF ( edge == 's' )  THEN
2906            DO  i = nxl, nxr + 1
2907               DO  k = kb(j,i), nzt + 1
2908                  f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
2909               ENDDO
2910            ENDDO
2911         ELSEIF ( edge == 'n' )  THEN
2912            DO  i = nxl, nxr + 1
2913               DO  k = kb(j,i), nzt + 1
2914                  f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
2915               ENDDO
2916            ENDDO
2917         ENDIF
2918      ENDIF
2919
2920!
2921!--   Store the boundary values also into the second ghost node layer.
2922      f(0:nzt+1,jb,nxl:nxr+1) = f(0:nzt+1,j,nxl:nxr+1)
2923
2924   END SUBROUTINE pmci_interp_tril_sn
2925
2926 
2927
2928   SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, var )
2929
2930!
2931!--   Interpolation of ghost-node values used as the client-domain boundary
2932!--   conditions. This subroutine handles the top boundary.
2933!--   This subroutine is based on trilinear interpolation.
2934!--   Constant dz is still assumed.
2935!
2936!--                                      Antti Hellsten 23.2.2015.
2937!
2938!
2939!--   Rewritten so that all the coefficients and client-array indices are
2940!--   precomputed in the initialization phase by pmci_init_interp_tril.
2941!
2942!--                                      Antti Hellsten 3.3.2015.
2943!
2944!--   Constant dz no more assumed.
2945!--                                      Antti Hellsten 23.3.2015.
2946!
2947!--   Indexing in the principal direction (k) is changed. Now, the nest-boundary
2948!--   values are interpolated only into the first ghost-node layer. Actually there is
2949!--   only one ghost-node layer in the k-direction. 
2950!
2951!--                                      Antti Hellsten 6.10.2015.
2952!
2953      IMPLICIT NONE
2954      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
2955      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)        ::  fc    !:
2956      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r1x   !:
2957      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r2x   !:
2958      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r1y   !:
2959      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r2y   !:
2960      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r1z   !:
2961      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r2z   !:
2962     
2963      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                    ::  ic    !:
2964      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                    ::  jc    !:
2965      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                    ::  kc    !:
2966     
2967      CHARACTER(LEN=1), INTENT(IN) :: var   !:
2968
2969      INTEGER(iwp) ::  i   !:
2970      INTEGER(iwp) ::  j   !:
2971      INTEGER(iwp) ::  k   !:
2972      INTEGER(iwp) ::  l   !:
2973      INTEGER(iwp) ::  m   !:
2974      INTEGER(iwp) ::  n   !:
2975     
2976      REAL(wp) ::  coarse_dx   !:
2977      REAL(wp) ::  coarse_dy   !:
2978      REAL(wp) ::  coarse_dz   !:
2979      REAL(wp) ::  fk          !:
2980      REAL(wp) ::  fkj         !:
2981      REAL(wp) ::  fkjp        !:
2982      REAL(wp) ::  fkpj        !:
2983      REAL(wp) ::  fkpjp       !:
2984      REAL(wp) ::  fkp         !:
2985
2986     
2987      IF ( var == 'w' )  THEN
2988         = nzt
2989      ELSE
2990         = nzt + 1
2991      ENDIF
2992     
2993      DO  i = nxl - 1, nxr + 1
2994         DO  j = nys - 1, nyn + 1
2995            l = ic(i)
2996            m = jc(j)
2997            n = kc(k)             
2998            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
2999            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
3000            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
3001            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
3002            fk       = r1y(j) * fkj  + r2y(j) * fkjp
3003            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
3004            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
3005         ENDDO
3006      ENDDO
3007
3008!
3009!--   Just fill up the second ghost-node layer for w.
3010      IF ( var == 'w' )  THEN
3011         f(nzt+1,:,:) = f(nzt,:,:)
3012      ENDIF
3013
3014!
3015!--   Rescale if f is the TKE.
3016!--   It is assumed that the bottom surface never reaches the top
3017!---  boundary of a nest domain.
3018      IF ( var == 'e')  THEN
3019         DO  i = nxl, nxr
3020            DO  j = nys, nyn
3021               f(k,j,i) = tkefactor_t(j,i) * f(k,j,i)
3022            ENDDO
3023         ENDDO
3024      ENDIF
3025
3026   END SUBROUTINE pmci_interp_tril_t
3027
3028
3029
3030   SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
3031
3032!
3033!--   After the interpolation of ghost-node values for the client-domain boundary
3034!--   conditions, this subroutine checks if there is a local outflow through the
3035!--   boundary. In that case this subroutine overwrites the interpolated values
3036!--   by values extrapolated from the domain. This subroutine handles the left and
3037!--   right boundaries.
3038!--   However, this operation is only needed in case of one-way coupling.
3039!
3040!--                                       Antti Hellsten 9.3.2015.
3041!
3042!--   Indexing in the principal direction (i) is changed. Now, the nest-boundary
3043!--   values are interpolated only into the first ghost-node layers on each later
3044!--   boundary. These values are then simply copied to the second ghost-node layer.
3045!
3046!--                                       Antti Hellsten 6.10.2015.
3047!
3048      IMPLICIT NONE
3049      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
3050      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
3051
3052      CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
3053      CHARACTER(LEN=1),INTENT(IN) ::  var    !:
3054
3055      INTEGER(iwp) ::  i     !:
3056      INTEGER(iwp) ::  ib    !:
3057      INTEGER(iwp) ::  ied   !:
3058      INTEGER(iwp) ::  j     !:
3059      INTEGER(iwp) ::  k     !:
3060     
3061      REAL(wp) ::  outnor    !:
3062      REAL(wp) ::  vdotnor   !:
3063
3064!
3065!--   Check which edge is to be handled: left or right.
3066      IF ( edge == 'l' )  THEN
3067         IF ( var == 'u' )  THEN
3068            i   = nxl
3069            ib  = nxl - 1
3070            ied = nxl + 1
3071         ELSE
3072            i   = nxl - 1
3073            ib  = nxl - 2
3074            ied = nxl
3075         ENDIF
3076         outnor = -1.0_wp
3077      ELSEIF ( edge == 'r' )  THEN
3078         i      = nxr + 1
3079         ib     = nxr + 2
3080         ied    = nxr
3081         outnor = 1.0_wp
3082      ENDIF
3083
3084      DO  j = nys, nyn + 1
3085         DO  k = kb(j,i), nzt +1
3086            vdotnor = outnor * u(k,j,ied)
3087            IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
3088               f(k,j,i) = f(k,j,ied)
3089            ENDIF
3090         ENDDO
3091         IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
3092            f(kb(j,i),j,i) = 0.0_wp
3093         ENDIF
3094      ENDDO
3095
3096!
3097!--   Store the updated boundary values also into the second ghost node layer.
3098      f(0:nzt,nys:nyn+1,ib) = f(0:nzt,nys:nyn+1,i)
3099
3100   END SUBROUTINE pmci_extrap_ifoutflow_lr
3101
3102
3103
3104   SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
3105!
3106!--   After  the interpolation of ghost-node values for the client-domain boundary
3107!--   conditions, this subroutine checks if there is a local outflow through the
3108!--   boundary. In that case this subroutine overwrites the interpolated values
3109!--   by values extrapolated from the domain. This subroutine handles the south and
3110!--   north boundaries.
3111!
3112!--                                       Antti Hellsten 9.3.2015.
3113!
3114!--   Indexing in the principal direction (j) is changed. Now, the nest-boundary
3115!--   values are interpolated only into the first ghost-node layers on each later
3116!--   boundary. These values are then simply copied to the second ghost-node layer.
3117!
3118!--                                       Antti Hellsten 6.10.2015.
3119
3120      IMPLICIT NONE
3121      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
3122      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
3123      CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
3124      CHARACTER(LEN=1), INTENT(IN) ::  var    !:
3125     
3126      INTEGER(iwp) ::  i         !:
3127      INTEGER(iwp) ::  j         !:
3128      INTEGER(iwp) ::  jb        !:
3129      INTEGER(iwp) ::  jed       !:
3130      INTEGER(iwp) ::  k         !:
3131      REAL(wp)     ::  outnor    !:
3132      REAL(wp)     ::  vdotnor   !:
3133
3134!
3135!--   Check which edge is to be handled: left or right.
3136      IF ( edge == 's' )  THEN
3137         IF ( var == 'v' )  THEN
3138            j   = nys
3139            jb  = nys - 1
3140            jed = nys + 1
3141         ELSE
3142            j   = nys - 1
3143            jb  = nys - 2
3144            jed = nys
3145         ENDIF
3146         outnor = -1.0_wp
3147      ELSEIF ( edge == 'n' )  THEN
3148         j      = nyn + 1
3149         jb     = nyn + 2
3150         jed    = nyn
3151         outnor = 1.0_wp
3152      ENDIF
3153
3154      DO  i = nxl, nxr + 1
3155         DO  k = kb(j,i), nzt + 1
3156            vdotnor = outnor * v(k,jed,i)
3157            IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
3158               f(k,j,i) = f(k,jed,i)                 
3159            ENDIF
3160         ENDDO
3161         IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
3162            f(kb(j,i),j,i) = 0.0_wp
3163         ENDIF
3164      ENDDO
3165
3166!
3167!--   Store the updated boundary values also into the second ghost node layer.
3168      f(0:nzt,jb,nxl:nxr+1) = f(0:nzt,j,nxl:nxr+1)
3169
3170   END SUBROUTINE pmci_extrap_ifoutflow_sn
3171
3172 
3173
3174   SUBROUTINE pmci_extrap_ifoutflow_t( f, var )
3175
3176!
3177!--   Interpolation of ghost-node values used as the client-domain boundary
3178!--   conditions. This subroutine handles the top boundary.
3179!--   This subroutine is based on trilinear interpolation.
3180!     
3181!--                                       Antti Hellsten 23.2.2015.
3182!     
3183!     
3184!--   Rewritten so that all the coefficients and client-array indices are
3185!--   precomputed in the initialization phase by init_interp_tril.
3186!     
3187!--                                       Antti Hellsten 3.3.2015.
3188!     
3189!--   Indexing in the principal direction (k) is changed. Now, the nest-boundary
3190!--   values are extrapolated only into the first ghost-node layer. Actually there is
3191!--   only one ghost-node layer in the k-direction. 
3192!     
3193!--                                       Antti Hellsten 6.10.2015.
3194      IMPLICIT NONE
3195      REAL(wp), DIMENSION(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp), INTENT(INOUT) ::  f   !:
3196      CHARACTER(LEN=1), INTENT(IN) ::  var   !:
3197     
3198      INTEGER(iwp) ::  i     !:
3199      INTEGER(iwp) ::  j     !:
3200      INTEGER(iwp) ::  k     !:
3201      INTEGER(iwp) ::  ked   !:
3202      REAL(wp) ::  vdotnor   !:
3203     
3204
3205      IF ( var == 'w' )  THEN
3206         k    = nzt
3207         ked  = nzt - 1
3208      ELSE
3209         k    = nzt + 1 
3210         ked  = nzt
3211      ENDIF
3212
3213      DO  i = nxl, nxr
3214         DO  j = nys, nyn
3215            vdotnor = w(ked,j,i)
3216            IF ( vdotnor > 0.0_wp )  THEN   !:  Local outflow.             
3217               f(k,j,i) = f(ked,j,i)
3218            ENDIF
3219         ENDDO
3220      ENDDO
3221
3222!
3223!--   Just fill up the second ghost-node layer for w.
3224      IF ( var == 'w' )  THEN
3225         f(nzt+1,:,:) = f(nzt,:,:)
3226      ENDIF
3227
3228   END SUBROUTINE pmci_extrap_ifoutflow_t
3229
3230
3231
3232   SUBROUTINE pmci_anterp_tophat( f, fc, kce, ifl, ifu, jfl, jfu, kfl, kfu, kb, var )
3233!
3234!--   Anterpolation of internal-node values to be used as the server-domain
3235!--   values. This subroutine is based on the first-order numerical
3236!--   integration of the fine-grid values contained within the coarse-grid
3237!--   cell.
3238!     
3239!--                                        Antti Hellsten 16.9.2015.
3240!
3241      IMPLICIT NONE
3242      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN)  ::  f     !:
3243      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)   ::  fc    !:
3244      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifl   !:
3245      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifu   !:
3246      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfl   !:
3247      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfu   !:
3248      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)        ::  kb    !: may be unnecessary
3249      INTEGER(iwp), INTENT(IN)                                        ::  kce   !:
3250      INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfl   !:
3251      INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfu   !:
3252      CHARACTER(LEN=1), INTENT(IN)                                    ::  var   !:                           
3253
3254      INTEGER(iwp) ::  i         !:
3255      INTEGER(iwp) ::  icb       !:
3256      INTEGER(iwp) ::  ice       !:
3257      INTEGER(iwp) ::  ifc       !:
3258      INTEGER(iwp) ::  ijfc      !:
3259      INTEGER(iwp) ::  j         !:
3260      INTEGER(iwp) ::  jcb       !:
3261      INTEGER(iwp) ::  jce       !:
3262      INTEGER(iwp) ::  k         !:
3263      INTEGER(iwp) ::  kcb       !:
3264      INTEGER(iwp) ::  l         !:
3265      INTEGER(iwp) ::  m         !:
3266      INTEGER(iwp) ::  n         !:
3267      INTEGER(iwp) ::  nfc       !:
3268      REAL(wp)     ::  cellsum   !:
3269      REAL(wp)     ::  f1f       !:
3270      REAL(wp)     ::  fra       !:
3271 
3272
3273      icb = icl
3274      ice = icr
3275      jcb = jcs
3276      jce = jcn
3277
3278!
3279!--   Define the index bounds icb, ice, jcb and jce.
3280!--   Note that kcb is simply zero and kce enters here as a parameter and it is
3281!--   determined in init_anterp_tophat 
3282      IF ( nest_bound_l )  THEN
3283         IF ( var == 'u' )  THEN
3284            icb = icl + nhll + 1
3285         ELSE
3286            icb = icl + nhll
3287         ENDIF
3288      ENDIF
3289      IF ( nest_bound_r )  THEN
3290         ice = icr - nhlr
3291      ENDIF
3292
3293      IF ( nest_bound_s )  THEN
3294         IF ( var == 'v' )  THEN
3295            jcb = jcs + nhls + 1
3296         ELSE
3297            jcb = jcs + nhls
3298         ENDIF
3299      ENDIF
3300      IF ( nest_bound_n )  THEN
3301         jce = jcn - nhln
3302      ENDIF
3303      kcb = 0 
3304
3305!
3306!--   Note that l,m, and n are coarse-grid indices and i,j, and k are fine-grid indices.
3307      DO  l = icb, ice
3308         ifc = ifu(l) - ifl(l) + 1
3309         DO  m = jcb, jce
3310            ijfc = ifc * ( jfu(m) - jfl(m) +1 )
3311
3312!
3313!--         How to deal with the lower bound of k in case of non-flat topography?
3314            !kcb = MIN( kb(jfl(m),ifl(l)), kb(jfu(m),ifl(l)), kb(jfl(m),ifu(l)), kb(jfu(m),ifu(l)) ) ! Something wrong with this.
3315            DO  n = kcb, kce
3316               nfc =  ijfc * ( kfu(n) - kfl(n) + 1 )
3317               cellsum = 0.0
3318               DO  i = ifl(l), ifu(l)
3319                  DO  j = jfl(m), jfu(m)
3320                     DO  k = kfl(n), kfu(n)
3321                        cellsum = cellsum + f(k,j,i)
3322                     ENDDO
3323                  ENDDO
3324               ENDDO
3325
3326!
3327!--            Spatial under-relaxation.
3328               fra  = frax(l) * fray(m) * fraz(n)
3329!--            TO_DO: why not KIND=wp ?
3330               fc(n,m,l) = ( 1.0_wp - fra ) * fc(n,m,l) + fra * cellsum / REAL( nfc, KIND=KIND(cellsum) ) 
3331            ENDDO
3332         ENDDO
3333      ENDDO
3334
3335   END SUBROUTINE pmci_anterp_tophat
3336
3337#endif
3338 END SUBROUTINE pmci_client_datatrans
3339
3340
3341
3342 SUBROUTINE pmci_update_new
3343
3344#if defined( __parallel )
3345!
3346!-- Copy the interpolated/anterpolated boundary values to the _p
3347!-- arrays, too, to make sure the interpolated/anterpolated boundary
3348!-- values are carried over from one RK inner step to another.
3349!-- So far works only with the cpp-switch __nopointer.
3350!
3351!--                 Antti Hellsten 8.3.2015
3352!
3353
3354!-- Just debugging
3355    w(nzt+1,:,:) = w(nzt,:,:)
3356
3357    u_p  = u
3358    v_p  = v
3359    w_p  = w
3360    e_p  = e
3361    pt_p = pt
3362    IF ( humidity  .OR.  passive_scalar )  THEN
3363       q_p = q
3364    ENDIF
3365
3366!
3367!-- TO_DO: Find out later if nesting would work without __nopointer.
3368#endif
3369
3370 END SUBROUTINE pmci_update_new
3371
3372
3373
3374 SUBROUTINE pmci_set_array_pointer( name, client_id, nz_cl )
3375
3376    IMPLICIT NONE
3377
3378    INTEGER, INTENT(IN)                  ::  client_id   !:
3379    INTEGER, INTENT(IN)                  ::  nz_cl       !:
3380    CHARACTER(LEN=*), INTENT(IN)         ::  name        !:
3381   
3382#if defined( __parallel )
3383    REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d        !:
3384    REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d_sec    !:
3385    REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d        !:
3386    REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d_sec    !:
3387    INTEGER(iwp)                         ::  ierr        !:
3388    INTEGER(iwp)                         ::  istat       !:
3389   
3390
3391    NULLIFY( p_3d )
3392    NULLIFY( p_2d )
3393
3394!
3395!-- List of array names, which can be coupled
3396    IF ( TRIM(name) == "u" )     p_3d => u
3397    IF ( TRIM(name) == "v" )     p_3d => v
3398    IF ( TRIM(name) == "w" )     p_3d => w
3399    IF ( TRIM(name) == "e" )     p_3d => e
3400    IF ( TRIM(name) == "pt" )    p_3d => pt
3401    IF ( TRIM(name) == "q" )     p_3d => q
3402    !IF ( TRIM(name) == "z0" )    p_2d => z0
3403
3404#if defined( __nopointer )
3405    IF ( ASSOCIATED( p_3d ) )  THEN
3406       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz )
3407    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3408       CALL pmc_s_set_dataarray( client_id, p_2d )
3409    ELSE
3410!
3411!--    Give only one message for the root domain
3412       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
3413
3414          message_string = 'pointer for array "' // TRIM( name ) //            &
3415                           '" can''t be associated'
3416          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
3417       ELSE
3418!
3419!--       Avoid others to continue
3420          CALL MPI_BARRIER( comm2d, ierr )
3421       ENDIF
3422    ENDIF
3423#else
3424!-- TO_DO: Why aren't the other pointers (p_3d) not set to u_1, v_1, etc.??
3425    IF ( TRIM(name) == "u" )     p_3d_sec => u_2
3426    IF ( TRIM(name) == "v" )     p_3d_sec => v_2
3427    IF ( TRIM(name) == "w" )     p_3d_sec => w_2
3428    IF ( TRIM(name) == "e" )     p_3d_sec => e_2
3429    IF ( TRIM(name) == "pt" )    p_3d_sec => pt_2
3430    !IF ( TRIM(name) == "z0" )    p_2d_sec => z0_2
3431
3432    IF ( ASSOCIATED( p_3d ) )  THEN
3433       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz, &
3434                                 array_2 = p_3d_sec )
3435    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3436       CALL pmc_s_set_dataarray( client_id, p_2d, array_2 = p_2d_sec )
3437    ELSE
3438!
3439!--    Give only one message for the root domain
3440       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
3441
3442          message_string = 'pointer for array "' // TRIM( name ) //            &
3443                           '" can''t be associated'
3444          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
3445       ELSE
3446!
3447!--       Avoid others to continue
3448          CALL MPI_BARRIER( comm2d, ierr )
3449       ENDIF
3450
3451   ENDIF
3452#endif
3453
3454#endif
3455 END SUBROUTINE pmci_set_array_pointer
3456
3457
3458
3459 SUBROUTINE pmci_create_client_arrays( name, is, ie, js, je, nzc  )
3460
3461    IMPLICIT NONE
3462
3463    INTEGER(iwp), INTENT(IN)              ::  ie      !:
3464    INTEGER(iwp), INTENT(IN)              ::  is      !:
3465    INTEGER(iwp), INTENT(IN)              ::  je      !:
3466    INTEGER(iwp), INTENT(IN)              ::  js      !:
3467    INTEGER(iwp), INTENT(IN)              ::  nzc     !:  Note that nzc is cg%nz
3468    CHARACTER(LEN=*), INTENT(IN)          ::  name    !:
3469     
3470#if defined( __parallel )
3471    REAL(wp), POINTER,DIMENSION(:,:)      ::  p_2d    !:
3472    REAL(wp), POINTER,DIMENSION(:,:,:)    ::  p_3d    !:
3473    INTEGER(iwp)                          ::  ierr    !:
3474    INTEGER(iwp)                          ::  istat   !:
3475
3476
3477    NULLIFY( p_3d )
3478    NULLIFY( p_2d )
3479
3480!
3481!-- List of array names, which can be coupled.
3482!-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1.
3483    IF ( TRIM( name ) == "u" )  THEN
3484       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
3485       p_3d => uc
3486    ELSEIF ( TRIM( name ) == "v" )  THEN
3487       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1, js:je, is:ie) )
3488       p_3d => vc
3489    ELSEIF ( TRIM( name ) == "w" )  THEN
3490       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1, js:je, is:ie) )
3491       p_3d => wc
3492    ELSEIF ( TRIM( name ) == "e" )  THEN
3493       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1, js:je, is:ie) )
3494       p_3d => ec
3495    ELSEIF ( TRIM( name ) == "pt")  THEN
3496       IF ( .NOT. ALLOCATED( ptc ) ) ALLOCATE( ptc(0:nzc+1, js:je, is:ie) )
3497       p_3d => ptc
3498    ELSEIF ( TRIM( name ) == "q")  THEN
3499       IF ( .NOT. ALLOCATED( qc ) ) ALLOCATE( qc(0:nzc+1, js:je, is:ie) )
3500       p_3d => qc
3501    !ELSEIF (trim(name) == "z0") then
3502       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
3503       !p_2d => z0c
3504    ENDIF
3505
3506    IF ( ASSOCIATED( p_3d ) )  THEN
3507       CALL pmc_c_set_dataarray( p_3d )
3508    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3509       CALL pmc_c_set_dataarray( p_2d )
3510    ELSE
3511!
3512!--    Give only one message for the first client domain
3513       IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
3514
3515          message_string = 'pointer for array "' // TRIM( name ) //            &
3516                           '" can''t be associated'
3517          CALL message( 'pmci_create_client_arrays', 'PA0170', 3, 2, 0, 6, 0 )
3518       ELSE
3519!
3520!--       Avoid others to continue
3521          CALL MPI_BARRIER( comm2d, ierr )
3522       ENDIF
3523    ENDIF
3524
3525#endif
3526 END SUBROUTINE pmci_create_client_arrays
3527
3528
3529
3530 SUBROUTINE pmci_server_initialize
3531
3532#if defined( __parallel )
3533    IMPLICIT NONE
3534
3535    INTEGER(iwp)   ::  client_id   !:
3536    INTEGER(iwp)   ::  m           !:
3537    REAL(wp)       ::  waittime    !:
3538
3539
3540    DO  m = 1, SIZE( pmc_server_for_client ) - 1
3541       client_id = pmc_server_for_client(m)
3542       CALL pmc_s_fillbuffer( client_id, waittime=waittime )
3543    ENDDO
3544
3545#endif
3546 END SUBROUTINE pmci_server_initialize
3547
3548
3549
3550 SUBROUTINE pmci_client_initialize
3551
3552#if defined( __parallel )
3553    IMPLICIT NONE
3554
3555    INTEGER(iwp)   ::  i          !:
3556    INTEGER(iwp)   ::  icl        !:
3557    INTEGER(iwp)   ::  icr        !:
3558    INTEGER(iwp)   ::  j          !:
3559    INTEGER(iwp)   ::  jcn        !:
3560    INTEGER(iwp)   ::  jcs        !:
3561    REAL(wp)       ::  waittime   !:
3562
3563
3564    IF ( cpl_id > 1 )  THEN   !  Root id is never a client     
3565
3566!
3567!--    Client domain boundaries in the server indice space.
3568       icl = coarse_bound(1)
3569       icr = coarse_bound(2)
3570       jcs = coarse_bound(3)
3571       jcn = coarse_bound(4)
3572
3573!
3574!--    Get data from server     
3575       CALL pmc_c_getbuffer( waittime = waittime )
3576
3577!
3578!--    The interpolation.
3579       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
3580       CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
3581       CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
3582       CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
3583       CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
3584       IF ( humidity  .OR.  passive_scalar )  THEN
3585          CALL pmci_interp_tril_all ( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
3586       ENDIF
3587
3588       IF ( topography /= 'flat' )  THEN
3589
3590!
3591!--       Inside buildings set velocities and TKE back to zero.
3592!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
3593!--       maybe revise later.
3594          DO   i = nxlg, nxrg
3595             DO   j = nysg, nyng
3596                u(nzb:nzb_u_inner(j,i),j,i)     = 0.0_wp
3597                v(nzb:nzb_v_inner(j,i),j,i)     = 0.0_wp
3598                w(nzb:nzb_w_inner(j,i),j,i)     = 0.0_wp
3599                e(nzb:nzb_w_inner(j,i),j,i)     = 0.0_wp             
3600                tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
3601                tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
3602                tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
3603                te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
3604             ENDDO
3605          ENDDO
3606       ENDIF
3607    ENDIF
3608   
3609
3610 CONTAINS
3611
3612
3613   SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, var )
3614
3615!
3616!-- Interpolation of the internal values for the client-domain initialization.
3617!-- This subroutine is based on trilinear interpolation.
3618!-- Coding based on interp_tril_lr/sn/t
3619!
3620!--                                     Antti Hellsten 20.10.2015.
3621      IMPLICIT NONE
3622      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
3623      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)        ::  fc    !:
3624      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r1x   !:
3625      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r2x   !:
3626      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r1y   !:
3627      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r2y   !:
3628      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r1z   !:
3629      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r2z   !:
3630      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                    ::  ic    !:
3631      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                    ::  jc    !:
3632      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                    ::  kc    !:
3633      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb    !:
3634      CHARACTER(LEN=1), INTENT(IN) :: var  !:
3635     
3636      INTEGER(iwp) ::  i      !:
3637      INTEGER(iwp) ::  ib     !:
3638      INTEGER(iwp) ::  ie     !:
3639      INTEGER(iwp) ::  j      !:
3640      INTEGER(iwp) ::  jb     !:
3641      INTEGER(iwp) ::  je     !:
3642      INTEGER(iwp) ::  k      !:
3643      INTEGER(iwp) ::  k1     !:
3644      INTEGER(iwp) ::  kbc    !:
3645      INTEGER(iwp) ::  l      !:
3646      INTEGER(iwp) ::  m      !:
3647      INTEGER(iwp) ::  n      !:
3648      REAL(wp) ::  fk         !:
3649      REAL(wp) ::  fkj        !:
3650      REAL(wp) ::  fkjp       !:
3651      REAL(wp) ::  fkp        !:
3652      REAL(wp) ::  fkpj       !:
3653      REAL(wp) ::  fkpjp      !:
3654      REAL(wp) ::  logratio   !:
3655      REAL(wp) ::  logzuc1    !:
3656      REAL(wp) ::  zuc1       !:
3657     
3658     
3659      ib = nxl
3660      ie = nxr
3661      jb = nys   
3662      je = nyn
3663      IF ( nest_bound_l )  THEN
3664         ib = nxl - 1 
3665         IF ( var == 'u' )  THEN   !  For u, nxl is a ghost node, but not for the other variables.
3666            ib = nxl
3667         ENDIF
3668      ENDIF
3669      IF ( nest_bound_s )  THEN
3670         jb = nys - 1 
3671         IF ( var == 'v' )  THEN   !  For v, nys is a ghost node, but not for the other variables.
3672            jb = nys
3673         ENDIF
3674      ENDIF
3675      IF ( nest_bound_r )  THEN
3676         ie = nxr + 1
3677      ENDIF
3678      IF ( nest_bound_n )  THEN
3679         je = nyn + 1
3680      ENDIF
3681
3682!
3683!--   Trilinear interpolation.
3684      DO  i = ib, ie
3685         DO  j = jb, je
3686            DO  k = kb(j,i), nzt + 1 
3687               l = ic(i)
3688               m = jc(j)
3689               n = kc(k)
3690               fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
3691               fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
3692               fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
3693               fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
3694               fk       = r1y(j) * fkj  + r2y(j) * fkjp
3695               fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
3696               f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
3697            ENDDO
3698         ENDDO
3699      ENDDO
3700
3701!
3702!--   Correct the interpolated values of u and v in near-wall nodes, i.e. in
3703!--   the nodes below the coarse-grid nodes with k=1. The corrction is only made
3704!--   over horizontal wall surfaces in this phase. For the nest boundary conditions,
3705!--   a corresponding corrections is made for all vertical walls, too.
3706      IF ( var == 'u' .OR. var == 'v' )  THEN
3707         DO  i = ib, nxr
3708            DO  j = jb, nyn             
3709               kbc = 1
3710               DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )   !  kbc is the first coarse-grid point above the surface.
3711                  kbc = kbc + 1
3712               ENDDO
3713               zuc1    = cg%zu(kbc)
3714               k1      = kb(j,i) + 1 
3715               DO  WHILE ( zu(k1) < zuc1 )
3716                  k1 = k1 + 1
3717               ENDDO
3718               logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
3719               
3720               k = kb(j,i) + 1
3721               DO  WHILE ( zu(k) < zuc1 ) 
3722                  logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) / logzuc1
3723                  f(k,j,i) = logratio * f(k1,j,i)
3724                  k  = k + 1
3725               ENDDO
3726               f(kb(j,i),j,i) = 0.0_wp
3727            ENDDO
3728         ENDDO
3729      ELSEIF ( var == 'w' )  THEN
3730         DO  i = ib, nxr
3731            DO  j = jb, nyn
3732               f(kb(j,i),j,i) = 0.0_wp
3733            ENDDO
3734         ENDDO
3735      ENDIF
3736     
3737   END SUBROUTINE pmci_interp_tril_all
3738
3739#endif
3740 END SUBROUTINE pmci_client_initialize
3741
3742
3743
3744 SUBROUTINE pmci_ensure_nest_mass_conservation
3745
3746#if defined( __parallel )
3747!
3748!-- Adjust the volume-flow rate through the top boundary
3749!-- so that the net volume flow through all boundaries
3750!-- of the current nest domain becomes zero.
3751    IMPLICIT NONE
3752
3753    INTEGER(iwp)  ::  i                          !:
3754    INTEGER(iwp)  ::  ierr                       !:
3755    INTEGER(iwp)  ::  j                          !:
3756    INTEGER(iwp)  ::  k                          !:
3757    REAL(wp) ::  dxdy                            !:
3758    REAL(wp) ::  innor                           !:
3759    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !:
3760    REAL(wp) ::  w_lt                            !:
3761
3762!
3763!-- Sum up the volume flow through the left/right boundaries.
3764    volume_flow(1)   = 0.0_wp
3765    volume_flow_l(1) = 0.0_wp
3766
3767    IF ( nest_bound_l )  THEN
3768       i = 0 
3769       innor = dy
3770       DO   j = nys, nyn
3771          DO   k = nzb_u_inner(j,i) + 1, nzt
3772             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)           
3773          ENDDO
3774       ENDDO
3775    ENDIF
3776
3777    IF ( nest_bound_r )  THEN     
3778       i = nx + 1
3779       innor = -dy
3780       DO   j = nys, nyn
3781          DO   k = nzb_u_inner(j,i) + 1, nzt
3782             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)           
3783          ENDDO
3784       ENDDO
3785    ENDIF
3786     
3787#if defined( __parallel )   
3788    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3789    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
3790                        MPI_SUM, comm2d, ierr )
3791#else
3792    volume_flow(1) = volume_flow_l(1)
3793#endif
3794
3795!
3796!-- Sum up the volume flow through the south/north boundaries.
3797    volume_flow(2)   = 0.0_wp
3798    volume_flow_l(2) = 0.0_wp
3799
3800    IF ( nest_bound_s )  THEN
3801       j = 0
3802       innor = dx
3803       DO   i = nxl, nxr
3804          DO   k = nzb_v_inner(j,i) + 1, nzt
3805             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
3806          ENDDO
3807       ENDDO
3808    ENDIF
3809
3810    IF ( nest_bound_n )  THEN
3811       j = ny + 1
3812       innor = -dx
3813       DO   i = nxl, nxr
3814          DO   k = nzb_v_inner(j,i)+1, nzt
3815             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
3816          ENDDO
3817       ENDDO
3818    ENDIF
3819
3820#if defined( __parallel )   
3821    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3822    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
3823                        MPI_SUM, comm2d, ierr )   
3824#else
3825    volume_flow(2) = volume_flow_l(2)
3826#endif
3827     
3828!
3829!-- Sum up the volume flow through the top boundary.
3830    volume_flow(3)   = 0.0_wp
3831    volume_flow_l(3) = 0.0_wp
3832    dxdy = dx * dy
3833    k = nzt
3834    DO   i = nxl, nxr
3835       DO   j = nys, nyn
3836          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
3837       ENDDO
3838    ENDDO
3839       
3840#if defined( __parallel )   
3841    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3842    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL, &
3843                        MPI_SUM, comm2d, ierr )   
3844#else
3845    volume_flow(3) = volume_flow_l(3)
3846#endif
3847
3848!
3849!-- Correct the top-boundary value of w.
3850    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
3851    DO   i = nxl, nxr
3852       DO   j = nys, nyn
3853          DO  k = nzt, nzt + 1
3854             w(k,j,i) = w(k,j,i) + w_lt
3855          ENDDO
3856       ENDDO
3857    ENDDO
3858
3859#endif
3860 END SUBROUTINE pmci_ensure_nest_mass_conservation
3861
3862END MODULE pmc_interface
Note: See TracBrowser for help on using the repository browser.