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

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

last commit documented

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