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

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