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

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago

added _mod string to several filenames to meet the naming convection for modules

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