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

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

update of the nested domain system + some bugfixes

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