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

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

last commit documented

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