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

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

output of nesting informations of all domains; filling up redundant ghost points; renaming of variables, etc.; formatting cleanup

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