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

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

last commit documented

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