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

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

last commit documented

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