source: palm/trunk/SOURCE/pmc_interface_mod.f90 @ 1882

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

Precomputation of ijfc for anterpolation added in pmc_interface_mod.f90

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