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

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

last commit documented

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