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

Last change on this file since 1895 was 1895, checked in by raasch, 6 years ago

last commit documented

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