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

Last change on this file since 1834 was 1818, checked in by maronga, 5 years ago

last commit documented / copyright update

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