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

Last change on this file since 1771 was 1767, checked in by raasch, 8 years ago

last commit documented

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