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

Last change on this file since 1808 was 1808, checked in by raasch, 5 years ago

preprocessor directives using machine dependent flags (lc, ibm, etc.) mostly removed from the code

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