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

Last change on this file since 1763 was 1763, checked in by hellstea, 6 years ago

last commit documented

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