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

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

last commit documented

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