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

Last change on this file since 1917 was 1901, checked in by raasch, 8 years ago

last commit documented

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