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

Last change on this file since 4313 was 4273, checked in by monakurppa, 20 months ago

Add logical switched nesting_chem and nesting_offline_chem

  • Property svn:keywords set to Id
File size: 221.3 KB
RevLine 
[2696]1!> @file pmc_interface_mod.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1762]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1762]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!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1762]19!
20! Current revisions:
21! ------------------
[2175]22!
[3946]23!
[2175]24! Former revisions:
25! -----------------
26! $Id: pmc_interface_mod.f90 4273 2019-10-24 13:40:54Z suehring $
[4273]27! Add a logical switch nesting_chem and rename nest_salsa to nesting_salsa
28!
29! 4260 2019-10-09 14:04:03Z hellstea
[4260]30! Rest of the possibly round-off-error sensitive grid-line matching tests
31! changed to round-off-error tolerant forms throughout the module.
32!
33! 4249 2019-10-01 12:27:47Z hellstea
[4249]34! Several grid-line matching tests changed to a round-off-error tolerant form
35! in pmci_setup_parent, pmci_define_index_mapping and pmci_check_grid_matching.
36!
37! 4182 2019-08-22 15:20:23Z scharf
[4182]38! Corrected "Former revisions" section
39!
40! 4168 2019-08-16 13:50:17Z suehring
[4168]41! Replace function get_topography_top_index by topo_top_ind
42!
43! 4029 2019-06-14 14:04:35Z raasch
[4029]44! nest_chemistry switch removed
45!
46! 4026 2019-06-12 16:50:15Z suehring
[4026]47! Masked topography at boundary grid points in mass conservation, in order to
48! avoid that mean velocities within topography are imposed
49!
50! 4011 2019-05-31 14:34:03Z hellstea
[4010]51! Mass (volume) flux correction included to ensure global mass conservation for child domains.
52!
53! 3987 2019-05-22 09:52:13Z kanani
[3987]54! Introduce alternative switch for debug output during timestepping
55!
56! 3984 2019-05-16 15:17:03Z hellstea
[3984]57! Commenting improved, pmci_map_fine_to_coarse_grid renamed as pmci_map_child_grid_to_parent_grid,
58! set_child_edge_coords renamed as pmci_set_child_edge_coords, some variables renamed, etc.
59!
60! 3979 2019-05-15 13:54:29Z hellstea
[3979]61! Bugfix in pmc_interp_1sto_sn. This bug had effect only in case of 1-d domain
62! decomposition with npex = 1.
63!
64! 3976 2019-05-15 11:02:34Z hellstea
[3976]65! Child initialization also for the redundant ghost points behind the nested
66! boundaries added (2nd and 3rd ghost-point layers and corners).
67!
68! 3948 2019-05-03 14:49:57Z hellstea
[3948]69! Some variables renamed, a little cleaning up and some commenting improvements
70!
71! 3947 2019-05-03 07:56:44Z hellstea
[3947]72! The checks included in 3946 are extended for the z-direction and moved into its
73! own subroutine called from pmci_define_index_mapping.
74!
75! 3946 2019-05-02 14:18:59Z hellstea
[3946]76! Check added for child domains too small in terms of number of parent-grid cells so
77! that anterpolation is not possible. Checks added for too wide anterpolation buffer
78! for the same reason. Some minor code reformatting done.
[3945]79!
[3946]80! 3945 2019-05-02 11:29:27Z raasch
81!
[3945]82! 3932 2019-04-24 17:31:34Z suehring
83! Add missing if statements for call of pmc_set_dataarray_name for TKE and
[3943]84! dissipation.
[3945]85!
86! 3888 2019-04-12 09:18:10Z hellstea
[3888]87! Variables renamed, commenting improved etc.
88!
89! 3885 2019-04-11 11:29:34Z kanani
[3885]90! Changes related to global restructuring of location messages and introduction
91! of additional debug messages
92!
93! 3883 2019-04-10 12:51:50Z hellstea
[3883]94! Checks and error messages improved and extended. All the child index bounds in the
95! parent-grid index space are made module variables. Function get_number_of_childs
96! renamed get_number_of_children. A number of variables renamed
97! and qite a lot of other code reshaping made all around the module.
98!
99! 3876 2019-04-08 18:41:49Z knoop
[3864]100! Implemented nesting for salsa variables.
101!
102! 3833 2019-03-28 15:04:04Z forkel
[3833]103! replaced USE chem_modules by USE chem_gasphase_mod
104!
105! 3822 2019-03-27 13:10:23Z hellstea
[3822]106! Temporary increase of the vertical dimension of the parent-grid arrays and
107! workarrc_t is cancelled as unnecessary.
108!
109! 3819 2019-03-27 11:01:36Z hellstea
[3819]110! Adjustable anterpolation buffer introduced on all nest boundaries, it is controlled
111! by the new nesting_parameters parameter anterpolation_buffer_width.
112!
113! 3804 2019-03-19 13:46:20Z hellstea
[3804]114! Anterpolation domain is lowered from kct-1 to kct-3 to avoid exessive       
115! kinetic energy from building up in CBL flows.
116!
117! 3803 2019-03-19 13:44:40Z hellstea
[3798]118! A bug fixed in lateral boundary interpolations. Dimension of val changed from 
119! 5 to 3 in pmci_setup_parent and pmci_setup_child.
120!
121! 3794 2019-03-15 09:36:33Z raasch
[3794]122! two remaining unused variables removed
123!
124! 3792 2019-03-14 16:50:07Z hellstea
[3792]125! Interpolations improved. Large number of obsolete subroutines removed.
126! All unused variables removed.
127!
128! 3741 2019-02-13 16:24:49Z hellstea
[3741]129! Interpolations and child initialization adjusted to handle set ups with child
130! pe-subdomain dimension not integer divisible by the grid-spacing ratio in the
131! respective direction. Set ups with pe-subdomain dimension smaller than the
132! grid-spacing ratio in the respective direction are now forbidden.
133!
134! 3708 2019-01-30 12:58:13Z hellstea
[3708]135! Checks for parent / child grid line matching introduced.
136! Interpolation of nest-boundary-tangential velocity components revised.
137!
138! 3697 2019-01-24 17:16:13Z hellstea
[3697]139! Bugfix: upper k-bound in the child initialization interpolation
140! pmci_interp_1sto_all corrected.
141! Copying of the nest boundary values into the redundant 2nd and 3rd ghost-node
142! layers is added to the pmci_interp_1sto_*-routines.
143!
144! 3681 2019-01-18 15:06:05Z hellstea
[3681]145! Linear interpolations are replaced by first order interpolations. The linear
146! interpolation routines are still included but not called. In the child
147! inititialization the interpolation is also changed to 1st order and the linear
148! interpolation is not kept.
149! Subroutine pmci_map_fine_to_coarse_grid is rewritten.
150! Several changes in pmci_init_anterp_tophat.
151! Child's parent-grid arrays (uc, vc,...) are made non-overlapping on the PE-
152! subdomain boundaries in order to allow grid-spacing ratios higher than nbgp.
153! Subroutine pmci_init_tkefactor is removed as unnecessary.
154!
155! 3655 2019-01-07 16:51:22Z knoop
[3646]156! Remove unused variable simulated_time
157!
[4182]158! 1762 2016-02-25 12:31:13Z hellstea
159! Initial revision by A. Hellsten
[2938]160!
[1762]161! Description:
162! ------------
163! Domain nesting interface routines. The low-level inter-domain communication   
164! is conducted by the PMC-library routines.
[2174]165!
166! @todo Remove array_3d variables from USE statements thate not used in the
167!       routine
168! @todo Data transfer of qc and nc is prepared but not activated
[2801]169!------------------------------------------------------------------------------!
[2696]170 MODULE pmc_interface
[1762]171
[2801]172    USE ISO_C_BINDING
173
174
[3681]175    USE arrays_3d,                                                             &
[2938]176        ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2,  &
[3182]177               pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                   &
[2938]178               u, u_p, u_2, v, v_p, v_2, w, w_p, w_2, zu, zw
[1762]179
[2801]180    USE control_parameters,                                                    &
[3182]181        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
[3274]182               bc_dirichlet_s, child_domain,                                   &
[2899]183               constant_diffusion, constant_flux_layer,                        &
[3885]184               coupling_char,                                                  &
[3987]185               debug_output_timestep,                                          &
[3885]186               dt_3d, dz, humidity, message_string,                            &
[3182]187               neutral, passive_scalar, rans_mode, rans_tke_e,                 &
[3864]188               roughness_length, salsa, topography, volume_flow
[1762]189
[3876]190    USE chem_gasphase_mod,                                                     &
[2773]191        ONLY:  nspec
192
[3876]193    USE chem_modules,                                                          &
[4273]194        ONLY:  chem_species, nesting_chem
[3876]195
[2773]196    USE chemistry_model_mod,                                                   &
[4029]197        ONLY:  spec_conc_2
[2773]198
199    USE cpulog,                                                                &
[1762]200        ONLY:  cpu_log, log_point_s
201
[2773]202    USE grid_variables,                                                        &
[1764]203        ONLY:  dx, dy
[1762]204
[2773]205    USE indices,                                                               &
206        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
[4168]207               nysv, nz, nzb, nzt, topo_top_ind, wall_flags_0
[1764]208
[3274]209    USE bulk_cloud_model_mod,                                                  &
210        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
211
[2801]212    USE particle_attributes,                                                   &
213        ONLY:  particle_advection
214
[1764]215    USE kinds
216
217#if defined( __parallel )
[2841]218#if !defined( __mpifh )
[1764]219    USE MPI
220#endif
221
[2773]222    USE pegrid,                                                                &
223        ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
[3681]224               numprocs, pdims, pleft, pnorth, pright, psouth, status
[1764]225
[2773]226    USE pmc_child,                                                             &
227        ONLY:  pmc_childinit, pmc_c_clear_next_array_list,                     &
228               pmc_c_getnextarray, pmc_c_get_2d_index_list, pmc_c_getbuffer,   &
229               pmc_c_putbuffer, pmc_c_setind_and_allocmem,                     &
[1764]230               pmc_c_set_dataarray, pmc_set_dataarray_name
231
[2773]232    USE pmc_general,                                                           &
[1900]233        ONLY:  da_namelen
[1764]234
[2773]235    USE pmc_handle_communicator,                                               &
236        ONLY:  pmc_get_model_info, pmc_init_model, pmc_is_rootmodel,           &
[2801]237               pmc_no_namelist_found, pmc_parent_for_child, m_couplers
[1764]238
[2773]239    USE pmc_mpi_wrapper,                                                       &
240        ONLY:  pmc_bcast, pmc_recv_from_child, pmc_recv_from_parent,           &
[1933]241               pmc_send_to_child, pmc_send_to_parent
[1764]242
[2773]243    USE pmc_parent,                                                            &
244        ONLY:  pmc_parentinit, pmc_s_clear_next_array_list, pmc_s_fillbuffer,  &
245               pmc_s_getdata_from_buffer, pmc_s_getnextarray,                  &
246               pmc_s_setind_and_allocmem, pmc_s_set_active_data_array,         &
[1779]247               pmc_s_set_dataarray, pmc_s_set_2d_index_list
[1764]248
249#endif
[3864]250   
251    USE salsa_mod,                                                             &
252        ONLY:  aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol,  &
[4273]253               ncomponents_mass, nconc_2, nesting_salsa, ngases_salsa,         &
254               salsa_gas, salsa_gases_from_chem
[1764]255
[2773]256    USE surface_mod,                                                           &
[4168]257        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
[2232]258
[1762]259    IMPLICIT NONE
260
[2841]261#if defined( __parallel )
262#if defined( __mpifh )
263    INCLUDE "mpif.h"
264#endif
265#endif
266
[1791]267    PRIVATE
[1762]268!
269!-- Constants
[3984]270    INTEGER(iwp), PARAMETER ::  child_to_parent = 2   !< Parameter for pmci_parent_datatrans indicating the direction of transfer
271    INTEGER(iwp), PARAMETER ::  parent_to_child = 1   !< Parameter for pmci_parent_datatrans indicating the direction of transfer
[3883]272    INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries
273    INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary
[4249]274
275    REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp                 !< Relative tolerence for grid-line matching tests and comparisons
[1762]276!
277!-- Coupler setup
[3984]278    INTEGER(iwp), SAVE      ::  comm_world_nesting    !< Global nesting communicator
279    INTEGER(iwp), SAVE      ::  cpl_id  = 1           !<
[2801]280    INTEGER(iwp), SAVE      ::  cpl_npe_total         !<
281    INTEGER(iwp), SAVE      ::  cpl_parent_id         !<
[3883]282   
283    CHARACTER(LEN=32), SAVE ::  cpl_name              !<
284
[1762]285!
[3484]286!-- Control parameters
[3883]287    INTEGER(iwp),     SAVE ::  anterpolation_buffer_width = 2       !< Boundary buffer width for anterpolation
[3792]288    CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !< steering parameter for data-transfer mode
289    CHARACTER(LEN=8), SAVE ::  nesting_mode = 'two-way'             !< steering parameter for 1- or 2-way nesting
[3819]290   
[2801]291    LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
[2938]292    LOGICAL       ::  rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode)
[1762]293!
294!-- Geometry
[3984]295    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_x            !< Array for the absolute x-coordinates
296    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_y            !< Array for the absolute y-coordinates
297    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_x !< x-coordinate of the lower left corner of the domain
298    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_y !< y-coordinate of the lower left corner of the domain
[1762]299!
[3484]300!-- Children's parent-grid arrays
[3948]301    INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC    ::  parent_bound        !< subdomain index bounds for children's parent-grid arrays
[1762]302
[3792]303    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< Parent-grid array on child domain - dissipation rate
304    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec    !< Parent-grid array on child domain - SGS TKE
305    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc   !< Parent-grid array on child domain - potential temperature
306    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uc    !< Parent-grid array on child domain - velocity component u
307    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vc    !< Parent-grid array on child domain - velocity component v
308    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc    !< Parent-grid array on child domain - velocity component w
309    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_c   !< Parent-grid array on child domain -
310    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc   !< Parent-grid array on child domain -
311    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qrc   !< Parent-grid array on child domain -
312    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nrc   !< Parent-grid array on child domain -
313    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc   !< Parent-grid array on child domain -
314    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc    !< Parent-grid array on child domain -
[2801]315    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  nr_partc    !<
316    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
[2773]317
[3883]318    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c  !< Parent-grid array on child domain - chemical species
[3864]319
[3883]320    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_mass_c    !< Aerosol mass
321    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_number_c  !< Aerosol number
322    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gas_c       !< SALSA gases
[1762]323!
[3792]324!-- Grid-spacing ratios.
[3883]325    INTEGER(iwp), SAVE ::  igsr    !< Integer grid-spacing ratio in i-direction
326    INTEGER(iwp), SAVE ::  jgsr    !< Integer grid-spacing ratio in j-direction
327    INTEGER(iwp), SAVE ::  kgsr    !< Integer grid-spacing ratio in k-direction
[3792]328!
[3819]329!-- Global parent-grid index bounds
[3883]330    INTEGER(iwp), SAVE ::  iplg    !< Leftmost parent-grid array ip index of the whole child domain
331    INTEGER(iwp), SAVE ::  iprg    !< Rightmost parent-grid array ip index of the whole child domain
332    INTEGER(iwp), SAVE ::  jpsg    !< Southmost parent-grid array jp index of the whole child domain
333    INTEGER(iwp), SAVE ::  jpng    !< Northmost parent-grid array jp index of the whole child domain
[3819]334!
[3883]335!-- Local parent-grid index bounds. Different sets of index bounds are needed for parent-grid arrays (uc, etc),
336!-- for index mapping arrays (iflu, etc) and for work arrays (workarr_lr, etc). This is because these arrays
337!-- have different dimensions depending on the location of the subdomain relative to boundaries and corners.
338    INTEGER(iwp), SAVE ::  ipl     !< Left index limit for children's parent-grid arrays
339    INTEGER(iwp), SAVE ::  ipla    !< Left index limit for allocation of index-mapping and other auxiliary arrays
340    INTEGER(iwp), SAVE ::  iplw    !< Left index limit for children's parent-grid work arrays
341    INTEGER(iwp), SAVE ::  ipr     !< Right index limit for children's parent-grid arrays
342    INTEGER(iwp), SAVE ::  ipra    !< Right index limit for allocation of index-mapping and other auxiliary arrays
343    INTEGER(iwp), SAVE ::  iprw    !< Right index limit for children's parent-grid work arrays
344    INTEGER(iwp), SAVE ::  jpn     !< North index limit for children's parent-grid arrays
345    INTEGER(iwp), SAVE ::  jpna    !< North index limit for allocation of index-mapping and other auxiliary arrays
346    INTEGER(iwp), SAVE ::  jpnw    !< North index limit for children's parent-grid work arrays
347    INTEGER(iwp), SAVE ::  jps     !< South index limit for children's parent-grid arrays
348    INTEGER(iwp), SAVE ::  jpsa    !< South index limit for allocation of index-mapping and other auxiliary arrays
349    INTEGER(iwp), SAVE ::  jpsw    !< South index limit for children's parent-grid work arrays
[3819]350!
[3792]351!-- Highest prognostic parent-grid k-indices.
[3681]352    INTEGER(iwp), SAVE ::  kcto     !< Upper bound for k in anterpolation of variables other than w.
353    INTEGER(iwp), SAVE ::  kctw     !< Upper bound for k in anterpolation of w.
[1762]354!
[1933]355!-- Child-array indices to be precomputed and stored for anterpolation.
[2997]356    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !< child index indicating left bound of parent grid box on u-grid
357    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !< child index indicating right bound of parent grid box on u-grid
358    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !< child index indicating left bound of parent grid box on scalar-grid
359    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !< child index indicating right bound of parent grid box on scalar-grid
360    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !< child index indicating south bound of parent grid box on v-grid
361    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !< child index indicating north bound of parent grid box on v-grid
362    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !< child index indicating south bound of parent grid box on scalar-grid
363    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !< child index indicating north bound of parent grid box on scalar-grid
364    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !< child index indicating lower bound of parent grid box on w-grid
365    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !< child index indicating upper bound of parent grid box on w-grid
366    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !< child index indicating lower bound of parent grid box on scalar-grid
367    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !< child index indicating upper bound of parent grid box on scalar-grid
[1882]368!
[3883]369!-- Number of child-grid nodes within anterpolation cells to be precomputed for anterpolation.
370    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid points contributing to a parent grid
371                                                                   !< node in anterpolation, u-grid
372    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid points contributing to a parent grid
373                                                                   !< node in anterpolation, v-grid
374    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid points contributing to a parent grid
375                                                                   !< node in anterpolation, w-grid
376    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid points contributing to a parent grid
377                                                                   !< node in anterpolation, scalar-grid
[3792]378!   
379!-- Work arrays for interpolation and user-defined type definitions for horizontal work-array exchange   
[3883]380    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_lr
381    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_sn
382    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_t
383
384    INTEGER(iwp) :: workarr_lr_exchange_type
385    INTEGER(iwp) :: workarr_sn_exchange_type
386    INTEGER(iwp) :: workarr_t_exchange_type_x
387    INTEGER(iwp) :: workarr_t_exchange_type_y
[3792]388 
[4010]389    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !< Array for communicating the parent-grid dimensions
390                                                                    !< to its children.
[3883]391
[4010]392    REAL(wp), DIMENSION(6)              ::  face_area               !< Surface area of each boundary face
393    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !< Array for communicating the real-type parent-grid
394                                                                    !< parameters to its children.
[1762]395
[3819]396    TYPE parentgrid_def
[2801]397       INTEGER(iwp)                        ::  nx                 !<
398       INTEGER(iwp)                        ::  ny                 !<
399       INTEGER(iwp)                        ::  nz                 !<
400       REAL(wp)                            ::  dx                 !<
401       REAL(wp)                            ::  dy                 !<
402       REAL(wp)                            ::  dz                 !<
403       REAL(wp)                            ::  lower_left_coord_x !<
404       REAL(wp)                            ::  lower_left_coord_y !<
405       REAL(wp)                            ::  xend               !<
406       REAL(wp)                            ::  yend               !<
407       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_x            !<
408       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_y            !<
409       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu                !<
410       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw                !<
411       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu                 !<
412       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw                 !<
[3819]413    END TYPE parentgrid_def
[1762]414
[3883]415    TYPE(parentgrid_def), SAVE, PUBLIC     ::  pg                 !< Parent-grid information package of type parentgrid_def
[3484]416!
417!-- Variables for particle coupling
[2801]418    TYPE, PUBLIC :: childgrid_def
419       INTEGER(iwp)                        ::  nx                   !<
420       INTEGER(iwp)                        ::  ny                   !<
421       INTEGER(iwp)                        ::  nz                   !<
422       REAL(wp)                            ::  dx                   !<
423       REAL(wp)                            ::  dy                   !<
424       REAL(wp)                            ::  dz                   !<
[3883]425       REAL(wp)                            ::  lx_coord, lx_coord_b !<   ! split onto separate lines
[2801]426       REAL(wp)                            ::  rx_coord, rx_coord_b !<
427       REAL(wp)                            ::  sy_coord, sy_coord_b !<
428       REAL(wp)                            ::  ny_coord, ny_coord_b !<
429       REAL(wp)                            ::  uz_coord, uz_coord_b !<
430    END TYPE childgrid_def
431
[3883]432    TYPE(childgrid_def), SAVE, ALLOCATABLE, DIMENSION(:), PUBLIC ::  childgrid  !<
[2801]433
[3883]434    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  nr_part  !<
435    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  part_adr !<
[3681]436
[3792]437   
[2311]438    INTERFACE pmci_boundary_conds
439       MODULE PROCEDURE pmci_boundary_conds
440    END INTERFACE pmci_boundary_conds
[2801]441   
[1797]442    INTERFACE pmci_check_setting_mismatches
443       MODULE PROCEDURE pmci_check_setting_mismatches
[1764]444    END INTERFACE
445
[1933]446    INTERFACE pmci_child_initialize
447       MODULE PROCEDURE pmci_child_initialize
[1764]448    END INTERFACE
449
[1878]450    INTERFACE pmci_synchronize
451       MODULE PROCEDURE pmci_synchronize
[1764]452    END INTERFACE
453
[1797]454    INTERFACE pmci_datatrans
455       MODULE PROCEDURE pmci_datatrans
456    END INTERFACE pmci_datatrans
457
[4010]458    INTERFACE pmci_ensure_nest_mass_conservation
459       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
460    END INTERFACE pmci_ensure_nest_mass_conservation
461
462    INTERFACE pmci_ensure_nest_mass_conservation_vertical
463       MODULE PROCEDURE pmci_ensure_nest_mass_conservation_vertical
464    END INTERFACE pmci_ensure_nest_mass_conservation_vertical
465
[1762]466    INTERFACE pmci_init
467       MODULE PROCEDURE pmci_init
468    END INTERFACE
[1764]469
[1762]470    INTERFACE pmci_modelconfiguration
471       MODULE PROCEDURE pmci_modelconfiguration
472    END INTERFACE
[1764]473
[1933]474    INTERFACE pmci_parent_initialize
475       MODULE PROCEDURE pmci_parent_initialize
[1764]476    END INTERFACE
477
[3883]478    INTERFACE get_number_of_children
479       MODULE PROCEDURE get_number_of_children
480    END  INTERFACE get_number_of_children
[2801]481
482    INTERFACE get_childid
483       MODULE PROCEDURE get_childid
484    END  INTERFACE get_childid
485
486    INTERFACE get_child_edges
487       MODULE PROCEDURE get_child_edges
488    END  INTERFACE get_child_edges
489
490    INTERFACE get_child_gridspacing
491       MODULE PROCEDURE get_child_gridspacing
492    END  INTERFACE get_child_gridspacing
493
[1766]494    INTERFACE pmci_set_swaplevel
495       MODULE PROCEDURE pmci_set_swaplevel
496    END INTERFACE pmci_set_swaplevel
497
[3792]498    PUBLIC child_to_parent, comm_world_nesting, cpl_id, nested_run,                                 &
499           nesting_datatransfer_mode, nesting_mode, parent_to_child, rans_mode_parent
[2311]500
501    PUBLIC pmci_boundary_conds
[1933]502    PUBLIC pmci_child_initialize
[1797]503    PUBLIC pmci_datatrans
[1762]504    PUBLIC pmci_init
505    PUBLIC pmci_modelconfiguration
[1933]506    PUBLIC pmci_parent_initialize
[1878]507    PUBLIC pmci_synchronize
[1766]508    PUBLIC pmci_set_swaplevel
[3883]509    PUBLIC get_number_of_children, get_childid, get_child_edges, get_child_gridspacing
[4010]510    PUBLIC pmci_ensure_nest_mass_conservation
511    PUBLIC pmci_ensure_nest_mass_conservation_vertical
512   
[1762]513 CONTAINS
514
515
516 SUBROUTINE pmci_init( world_comm )
[1764]517
[4260]518    USE control_parameters,                                                                         &
[1797]519        ONLY:  message_string
520
[1762]521    IMPLICIT NONE
522
[2801]523    INTEGER(iwp), INTENT(OUT) ::  world_comm   !<
[1762]524
[1764]525#if defined( __parallel )
[1762]526
[3241]527    INTEGER(iwp) ::  pmc_status   !<
[1762]528
[1764]529
[4260]530    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,                       &
[3819]531                         anterpolation_buffer_width, pmc_status )
[1764]532
533    IF ( pmc_status == pmc_no_namelist_found )  THEN
534!
535!--    This is not a nested run
536       world_comm = MPI_COMM_WORLD
537       cpl_id     = 1
538       cpl_name   = ""
[1779]539
[1764]540       RETURN
[1779]541
[1762]542    ENDIF
[1779]543!
[1797]544!-- Check steering parameter values
[4260]545    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                                                   &
546         TRIM( nesting_mode ) /= 'two-way'  .AND.                                                   &
547         TRIM( nesting_mode ) /= 'vertical' )                                                       &
[1797]548    THEN
549       message_string = 'illegal nesting mode: ' // TRIM( nesting_mode )
550       CALL message( 'pmci_init', 'PA0417', 3, 2, 0, 6, 0 )
551    ENDIF
552
[4260]553    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                                      &
554         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                                      &
555         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                                           &
[1797]556    THEN
[4260]557       message_string = 'illegal nesting datatransfer mode: ' // TRIM( nesting_datatransfer_mode )
[1797]558       CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
559    ENDIF
560!
[1779]561!-- Set the general steering switch which tells PALM that its a nested run
562    nested_run = .TRUE.
563!
564!-- Get some variables required by the pmc-interface (and in some cases in the
565!-- PALM code out of the pmci) out of the pmc-core
[4260]566    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,                               &
567                             cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,                        &
568                             cpl_name = cpl_name, npe_total = cpl_npe_total,                        &
569                             lower_left_x = lower_left_coord_x,                                     &
[1791]570                             lower_left_y = lower_left_coord_y )
[1764]571!
[1779]572!-- Set the steering switch which tells the models that they are nested (of
[3948]573!-- course the root domain is not nested)
574    IF ( .NOT.  pmc_is_rootmodel() )  THEN
[3182]575       child_domain = .TRUE.
[2669]576       WRITE( coupling_char, '(A2,I2.2)') '_N', cpl_id
[1779]577    ENDIF
578
579!
[1764]580!-- Message that communicators for nesting are initialized.
581!-- Attention: myid has been set at the end of pmc_init_model in order to
582!-- guarantee that only PE0 of the root domain does the output.
[3885]583    CALL location_message( 'initialize model nesting', 'finished' )
[1764]584!
585!-- Reset myid to its default value
586    myid = 0
[1762]587#else
[1764]588!
589!-- Nesting cannot be used in serial mode. cpl_id is set to root domain (1)
590!-- because no location messages would be generated otherwise.
591!-- world_comm is given a dummy value to avoid compiler warnings (INTENT(OUT)
[3948]592!-- must get an explicit value).
593!-- Note that this branch is only to avoid compiler warnings. The actual
594!-- execution never reaches here because the call of this subroutine is
595!-- already enclosed by  #if defined( __parallel ).
[1762]596    cpl_id     = 1
[1764]597    nested_run = .FALSE.
598    world_comm = 1
[1762]599#endif
600
601 END SUBROUTINE pmci_init
602
603
604
605 SUBROUTINE pmci_modelconfiguration
[1764]606
[1762]607    IMPLICIT NONE
608
[2669]609    INTEGER(iwp) ::  ncpl   !<  number of nest domains
610
[3883]611   
[2967]612#if defined( __parallel )
[3885]613    CALL location_message( 'setup the nested model configuration', 'start' )
[2951]614    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'start' )
[1791]615!
616!-- Compute absolute coordinates for all models
[3883]617    CALL pmci_setup_coordinates         ! CONTAIN THIS
[1791]618!
[3592]619!-- Determine the number of coupled arrays
[3883]620    CALL pmci_num_arrays                ! CONTAIN THIS
[3592]621!
[1933]622!-- Initialize the child (must be called before pmc_setup_parent)
[3984]623!-- Klaus, extend this comment to explain why it must be called before   
[3883]624    CALL pmci_setup_child               ! CONTAIN THIS
[1791]625!
[1933]626!-- Initialize PMC parent
[3883]627    CALL pmci_setup_parent              ! CONTAIN THIS
[1797]628!
[1933]629!-- Check for mismatches between settings of master and child variables
630!-- (e.g., all children have to follow the end_time settings of the root master)
[3883]631    CALL pmci_check_setting_mismatches  ! CONTAIN THIS
[2669]632!
633!-- Set flag file for combine_plot_fields for precessing the nest output data
634    OPEN( 90, FILE='3DNESTING', FORM='FORMATTED' )
635    CALL pmc_get_model_info( ncpl = ncpl )
636    WRITE( 90, '(I2)' )  ncpl
637    CLOSE( 90 )
[1791]638
[2951]639    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'stop' )
[3885]640    CALL location_message( 'setup the nested model configuration', 'finished' )
[2967]641#endif
[1762]642
643 END SUBROUTINE pmci_modelconfiguration
644
645
646
[1933]647 SUBROUTINE pmci_setup_parent
[1764]648
649#if defined( __parallel )
[1762]650    IMPLICIT NONE
651
[3883]652    INTEGER(iwp) ::  child_id           !< Child id-number for the child m
653    INTEGER(iwp) ::  ierr               !< MPI-error code
654    INTEGER(iwp) ::  kp                 !< Parent-grid index n the z-direction
655    INTEGER(iwp) ::  lb = 1             !< Running index for aerosol size bins
656    INTEGER(iwp) ::  lc = 1             !< Running index for aerosol mass bins
657    INTEGER(iwp) ::  lg = 1             !< Running index for SALSA gases
658    INTEGER(iwp) ::  m                  !< Loop index over all children of the current parent
659    INTEGER(iwp) ::  msib               !< Loop index over all other children than m in case of siblings (parallel children)
660    INTEGER(iwp) ::  n = 1              !< Running index for chemical species
661    INTEGER(iwp) ::  nest_overlap = 0   !< Tag for parallel child-domains' overlap situation (>0 if overlap found)
662    INTEGER(iwp) ::  nomatch = 0        !< Tag for child-domain mismatch situation (>0 if mismatch found)
663    INTEGER(iwp) ::  nx_child           !< Number of child-grid points in the x-direction
664    INTEGER(iwp) ::  ny_child           !< Number of child-grid points in the y-direction
665    INTEGER(iwp) ::  nz_child           !< Number of child-grid points in the z-direction
666   
667    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for receiving the child-grid dimensions from the children
668   
669    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_left   !< Minimum x-coordinate of the child domain including the ghost
670                                                           !< point layers
671    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_right  !< Maximum x-coordinate of the child domain including the ghost
672                                                           !< point layers   
673    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_south  !< Minimum y-coordinate of the child domain including the ghost
674                                                           !< point layers
675    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_north  !< Maximum y-coordinate of the child domain including the ghost
676                                                           !< point layers
677    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_x  !< Child domain x-coordinate array
678    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_y  !< Child domain y-coordinate array
679   
680    REAL(wp), DIMENSION(5) ::  child_grid_info  !< Array for receiving the child-grid spacings etc from the children
681   
682    REAL(wp) ::  child_height         !< Height of the child domain defined on the child side as zw(nzt+1)
683    REAL(wp) ::  dx_child             !< Child-grid spacing in the x-direction
684    REAL(wp) ::  dy_child             !< Child-grid spacing in the y-direction
685    REAL(wp) ::  dz_child             !< Child-grid spacing in the z-direction
686    REAL(wp) ::  left_limit           !< Left limit for the absolute x-coordinate of the child left boundary
687    REAL(wp) ::  north_limit          !< North limit for the absolute y-coordinate of the child north boundary
688    REAL(wp) ::  right_limit          !< Right limit for the absolute x-coordinate of the child right boundary
689    REAL(wp) ::  south_limit          !< South limit for the absolute y-coordinate of the child south boundary
690    REAL(wp) ::  upper_right_coord_x  !< Absolute x-coordinate of the upper right corner of the child domain
691    REAL(wp) ::  upper_right_coord_y  !< Absolute y-coordinate of the upper right corner of the child domain 
692    REAL(wp) ::  xez                  !< Minimum separation in the x-direction required between the child and
693                                      !< parent boundaries (left or right)
694    REAL(wp) ::  yez                  !< Minimum separation in the y-direction required between the child and
695                                      !< parent boundaries (south or north)
[4260]696    REAL(wp)     ::  tolex            !< Tolerance for grid-line matching in x-direction
697    REAL(wp)     ::  toley            !< Tolerance for grid-line matching in y-direction
[4249]698    REAL(wp)     ::  tolez            !< Tolerance for grid-line matching in z-direction   
699
[3883]700    CHARACTER(LEN=32) ::  myname      !< String for variable name such as 'u'
[3864]701
[3883]702    LOGICAL :: m_left_in_msib         !< Logical auxiliary parameter for the overlap test: true if the left border
703                                      !< of the child m is within the x-range of the child msib
704    LOGICAL :: m_right_in_msib        !< Logical auxiliary parameter for the overlap test: true if the right border
705                                      !< of the child m is within the x-range of the child msib
706    LOGICAL :: msib_left_in_m         !< Logical auxiliary parameter for the overlap test: true if the left border
707                                      !< of the child msib is within the x-range of the child m
708    LOGICAL :: msib_right_in_m        !< Logical auxiliary parameter for the overlap test: true if the right border
709                                      !< of the child msib is within the x-range of the child m
710    LOGICAL :: m_south_in_msib        !< Logical auxiliary parameter for the overlap test: true if the south border
711                                      !< of the child m is within the y-range of the child msib
712    LOGICAL :: m_north_in_msib        !< Logical auxiliary parameter for the overlap test: true if the north border
713                                      !< of the child m is within the y-range of the child msib
714    LOGICAL :: msib_south_in_m        !< Logical auxiliary parameter for the overlap test: true if the south border
715                                      !< of the child msib is within the y-range of the child m
716    LOGICAL :: msib_north_in_m        !< Logical auxiliary parameter for the overlap test: true if the north border
717                                      !< of the child msib is within the y-range of the child m
[4249]718
719
[1764]720!
[4249]721!-- Grid-line tolerances.
722    tolex = tolefac * dx
723    toley = tolefac * dy
724    tolez = tolefac * dz(1)   
725!
[3948]726!-- Initialize the current pmc parent.
[1933]727    CALL pmc_parentinit
[1762]728!
[3883]729!-- Corners of all children of the present parent. Note that
730!-- SIZE( pmc_parent_for_child ) = 1 if we have no children.
731    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 )  .AND.  myid == 0 )  THEN
732       ALLOCATE( child_x_left(1:SIZE( pmc_parent_for_child ) - 1) )
733       ALLOCATE( child_x_right(1:SIZE( pmc_parent_for_child ) - 1) )
734       ALLOCATE( child_y_south(1:SIZE( pmc_parent_for_child ) - 1) )
735       ALLOCATE( child_y_north(1:SIZE( pmc_parent_for_child ) - 1) )
[1925]736    ENDIF
[2801]737    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) )  THEN
738       ALLOCATE( childgrid(1:SIZE( pmc_parent_for_child ) - 1) )
739    ENDIF
[1925]740!
[3883]741!-- Get coordinates from all children and check that the children match the parent
742!-- domain and each others. Note that SIZE( pmc_parent_for_child ) = 1
743!-- if we have no children, thence the loop is not executed at all.
[1933]744    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
[1791]745
[1933]746       child_id = pmc_parent_for_child(m)
[1762]747
[2938]748       IF ( myid == 0 )  THEN
749
[3883]750          CALL pmc_recv_from_child( child_id, child_grid_dim,  SIZE(child_grid_dim), 0, 123, ierr )
751          CALL pmc_recv_from_child( child_id, child_grid_info, SIZE(child_grid_info), 0, 124, ierr )
[1762]752         
[3883]753          nx_child     = child_grid_dim(1)
754          ny_child     = child_grid_dim(2)
755          dx_child     = child_grid_info(3)
756          dy_child     = child_grid_info(4)
757          dz_child     = child_grid_info(5)
758          child_height = child_grid_info(1)
[1762]759!
[4260]760!--       Find the highest child-domain level in the parent grid for the reduced z transfer
[3948]761          DO  kp = 1, nzt                 
[4249]762             IF ( zw(kp) - child_height > tolez )  THEN                   
[3883]763                nz_child = kp
[1762]764                EXIT
765             ENDIF
766          ENDDO
767!   
[1925]768!--       Get absolute coordinates from the child
[3883]769          ALLOCATE( child_coord_x(-nbgp:nx_child+nbgp) )
770          ALLOCATE( child_coord_y(-nbgp:ny_child+nbgp) )
[1762]771         
[3883]772          CALL pmc_recv_from_child( child_id, child_coord_x, SIZE( child_coord_x ), 0, 11, ierr )
773          CALL pmc_recv_from_child( child_id, child_coord_y, SIZE( child_coord_y ), 0, 12, ierr )
[1762]774         
[2602]775          parent_grid_info_real(1) = lower_left_coord_x
776          parent_grid_info_real(2) = lower_left_coord_y
777          parent_grid_info_real(3) = dx
778          parent_grid_info_real(4) = dy
[3883]779
780          upper_right_coord_x      = lower_left_coord_x + ( nx + 1 ) * dx
781          upper_right_coord_y      = lower_left_coord_y + ( ny + 1 ) * dy
782          parent_grid_info_real(5) = upper_right_coord_x
783          parent_grid_info_real(6) = upper_right_coord_y
[3065]784          parent_grid_info_real(7) = dz(1)
[1762]785
[2602]786          parent_grid_info_int(1)  = nx
787          parent_grid_info_int(2)  = ny
[3883]788          parent_grid_info_int(3)  = nz_child
[1762]789!
[3883]790!--       Check that the child domain matches its parent domain.
[1933]791          IF ( nesting_mode == 'vertical' )  THEN
[3883]792!
793!--          In case of vertical nesting, the lateral boundaries must match exactly.
794             right_limit = upper_right_coord_x
795             north_limit = upper_right_coord_y
[4260]796             IF ( ( ABS( child_coord_x(nx_child+1) - right_limit ) > tolex )  .OR.                  &
797                  ( ABS( child_coord_y(ny_child+1) - north_limit ) > toley ) )  THEN               
[1933]798                nomatch = 1
799             ENDIF
[2602]800          ELSE       
[1933]801!
[3883]802!--          In case of 3-D nesting, check that the child domain is completely
803!--          inside its parent domain.
[1933]804             xez = ( nbgp + 1 ) * dx
805             yez = ( nbgp + 1 ) * dy
806             left_limit  = lower_left_coord_x + xez
[3883]807             right_limit = upper_right_coord_x - xez
[1933]808             south_limit = lower_left_coord_y + yez
[3883]809             north_limit = upper_right_coord_y - yez
[4249]810             IF ( ( left_limit - child_coord_x(0) > tolex )  .OR.                                   &
811                  ( child_coord_x(nx_child+1) - right_limit > tolex )  .OR.                         &
812                  ( south_limit - child_coord_y(0) > toley )  .OR.                                  &
813                  ( child_coord_y(ny_child+1) - north_limit > toley ) )  THEN
[1933]814                nomatch = 1
815             ENDIF
[1791]816          ENDIF
[2602]817!             
[4260]818!--       Child domain must be lower than the parent domain such that the top ghost
819!--       layer of the child grid does not exceed the parent domain top boundary.
820          IF ( child_height - zw(nzt) > tolez ) THEN
[2602]821             nomatch = 1
822          ENDIF
[1925]823!
[3883]824!--       If parallel child domains (siblings) do exist ( m > 1 ),
825!--       check that they do not overlap.
826          child_x_left(m)  = child_coord_x(-nbgp)
827          child_x_right(m) = child_coord_x(nx_child+nbgp)
828          child_y_south(m) = child_coord_y(-nbgp)
829          child_y_north(m) = child_coord_y(ny_child+nbgp)
[1925]830
[3883]831          IF ( nesting_mode /= 'vertical' )  THEN
[2801]832!
[3883]833!--          Note that the msib-loop is executed only if ( m > 1 ). 
834!--          Also note that the tests have to be made both ways (m vs msib and msib vs m)
835!--          in order to detect all the possible overlap situations.
836             DO  msib = 1, m - 1
837!
[4260]838!--             Set some logical auxiliary parameters to simplify the IF-condition.                 
839                m_left_in_msib  = ( child_x_left(m)  >= child_x_left(msib)  - tolex )  .AND.        &
840                                  ( child_x_left(m)  <= child_x_right(msib) + tolex )
841                m_right_in_msib = ( child_x_right(m) >= child_x_left(msib)  - tolex )  .AND.        &
842                                  ( child_x_right(m) <= child_x_right(msib) + tolex )
843                msib_left_in_m  = ( child_x_left(msib)  >= child_x_left(m)  - tolex )  .AND.        &
844                                  ( child_x_left(msib)  <= child_x_right(m) + tolex )
845                msib_right_in_m = ( child_x_right(msib) >= child_x_left(m)  - tolex )  .AND.        &
846                                  ( child_x_right(msib) <= child_x_right(m) + tolex )
847                m_south_in_msib = ( child_y_south(m) >= child_y_south(msib) - toley )  .AND.        &
848                                  ( child_y_south(m) <= child_y_north(msib) + toley )
849                m_north_in_msib = ( child_y_north(m) >= child_y_south(msib) - toley )  .AND.        &
850                                  ( child_y_north(m) <= child_y_north(msib) + toley )
851                msib_south_in_m = ( child_y_south(msib) >= child_y_south(m) - toley )  .AND.        &
852                                  ( child_y_south(msib) <= child_y_north(m) + toley )
853                msib_north_in_m = ( child_y_north(msib) >= child_y_south(m) - toley )  .AND.        &
854                                  ( child_y_north(msib) <= child_y_north(m) + toley )
[3883]855               
856                IF ( ( m_left_in_msib  .OR.  m_right_in_msib  .OR.                                  &
857                       msib_left_in_m  .OR.  msib_right_in_m )                                      &
858                     .AND.                                                                          &
859                     ( m_south_in_msib  .OR.  m_north_in_msib  .OR.                                 &
860                       msib_south_in_m  .OR.  msib_north_in_m ) )  THEN
861                     nest_overlap = 1
862                ENDIF
[1925]863
[3883]864             ENDDO
865          ENDIF         
866
[3984]867          CALL pmci_set_child_edge_coords
[2801]868
[3883]869          DEALLOCATE( child_coord_x )
870          DEALLOCATE( child_coord_y )
[1762]871!
[2938]872!--       Send information about operating mode (LES or RANS) to child. This will be
873!--       used to control TKE nesting and setting boundary conditions properly.
874          CALL pmc_send_to_child( child_id, rans_mode, 1, 0, 19, ierr ) 
875!
[3984]876!--       Send parent grid information to child
[3883]877          CALL pmc_send_to_child( child_id, parent_grid_info_real,                                  &
878                                  SIZE( parent_grid_info_real ), 0, 21,                             &
879                                  ierr )
880          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,                            &
881                                  22, ierr )
[1762]882!
[1925]883!--       Send local grid to child
[3883]884          CALL pmc_send_to_child( child_id, coord_x, nx+1+2*nbgp, 0, 24,                            &
885                                  ierr )
886          CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,                            &
887                                  ierr )
[1762]888!
[1764]889!--       Also send the dzu-, dzw-, zu- and zw-arrays here
[3883]890          CALL pmc_send_to_child( child_id, dzu, nz_child + 1, 0, 26, ierr )
891          CALL pmc_send_to_child( child_id, dzw, nz_child + 1, 0, 27, ierr )
892          CALL pmc_send_to_child( child_id, zu,  nz_child + 2, 0, 28, ierr )
893          CALL pmc_send_to_child( child_id, zw,  nz_child + 2, 0, 29, ierr )
[1762]894
[3883]895          IF ( nomatch /= 0 )  THEN
896             WRITE ( message_string, * ) 'nested child domain does not fit into its parent domain'
897             CALL message( 'pmci_setup_parent', 'PA0425', 3, 2, 0, 6, 0 )
898          ENDIF
[1925]899 
[3883]900          IF ( nest_overlap /= 0  .AND.  nesting_mode /= 'vertical' )  THEN
901             WRITE ( message_string, * ) 'nested parallel child domains overlap'
902             CALL message( 'pmci_setup_parent', 'PA0426', 3, 2, 0, 6, 0 )
903          ENDIF
904         
905       ENDIF  ! ( myid == 0 )
[2801]906
[3883]907       CALL MPI_BCAST( nz_child, 1, MPI_INTEGER, 0, comm2d, ierr )
908
[2809]909       CALL MPI_BCAST( childgrid(m), STORAGE_SIZE(childgrid(1))/8, MPI_BYTE, 0, comm2d, ierr )
[1791]910!
[3888]911!--    Set up the index-list which is an integer array that maps the child index space on
912!--    the parent index- and subdomain spaces.
[1762]913       CALL pmci_create_index_list
914!
[3946]915!--    Include couple arrays into parent content.
916!--    The adresses of the PALM 2D or 3D array (here parent grid) which are candidates
[2801]917!--    for coupling are stored once into the pmc context. While data transfer, the array do not
918!--    have to be specified again
[1779]919       CALL pmc_s_clear_next_array_list
[3883]920       DO WHILE ( pmc_s_getnextarray( child_id, myname ) )
[2773]921          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
[3946]922             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = n )
[3864]923             n = n + 1 
924          ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
[3946]925             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lb )
[3883]926             lb = lb + 1 
[3864]927          ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
[3946]928             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lc )
[3883]929             lc = lc + 1 
[3946]930          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.  .NOT. salsa_gases_from_chem )  THEN
931             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child, n = lg )
[3883]932             lg = lg + 1
[2773]933          ELSE
[3883]934             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child )
[2773]935          ENDIF
[1762]936       ENDDO
[2938]937
[1933]938       CALL pmc_s_setind_and_allocmem( child_id )
[3883]939       
940    ENDDO  ! m
[1762]941
[3883]942    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN
943       DEALLOCATE( child_x_left )
944       DEALLOCATE( child_x_right )
945       DEALLOCATE( child_y_south )
946       DEALLOCATE( child_y_north )
[1925]947    ENDIF
948
[3883]949   
[1762]950 CONTAINS
951
952
[3792]953    SUBROUTINE pmci_create_index_list
[1764]954
[1762]955       IMPLICIT NONE
[1764]956
[3946]957       INTEGER(iwp) ::  ilist            !< Index-list index running over the child's parent-grid jc,ic-space
958       INTEGER(iwp) ::  index_list_size  !< Dimension 2 of the array index_list
959       INTEGER(iwp) ::  ierr             !< MPI error code
960       INTEGER(iwp) ::  ip               !< Running parent-grid index on the child domain in the x-direction
961       INTEGER(iwp) ::  jp               !< Running parent-grid index on the child domain in the y-direction
962       INTEGER(iwp) ::  n                !< Running index over child subdomains
963       INTEGER(iwp) ::  nrx              !< Parent subdomain dimension in the x-direction
964       INTEGER(iwp) ::  nry              !< Parent subdomain dimension in the y-direction
965       INTEGER(iwp) ::  pex              !< Two-dimensional subdomain (pe) index in the x-direction
966       INTEGER(iwp) ::  pey              !< Two-dimensional subdomain (pe) index in the y-direction
967       INTEGER(iwp) ::  parent_pe        !< Parent subdomain index (one-dimensional)
[1791]968
[3888]969       INTEGER(iwp), DIMENSION(2) ::  pe_indices_2d                                  !< Array for two-dimensional subdomain (pe)
970                                                                                     !< indices needed for MPI_CART_RANK
971       INTEGER(iwp), DIMENSION(2) ::  size_of_childs_parent_grid_bounds_all          !< Dimensions of childs_parent_grid_bounds_all
972       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  childs_parent_grid_bounds_all  !< Array that contains the child's
973                                                                                     !< parent-grid index bounds for all its
974                                                                                     !< subdomains (pes)
975       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list                     !< Array that maps the child index space on
976                                                                                     !< the parent index- and subdomain spaces
[3883]977       
[1762]978       IF ( myid == 0 )  THEN
[3888]979         
980          CALL pmc_recv_from_child( child_id, size_of_childs_parent_grid_bounds_all,                &
981                                    2, 0, 40, ierr )
982          ALLOCATE( childs_parent_grid_bounds_all(size_of_childs_parent_grid_bounds_all(1),         &
983                                                  size_of_childs_parent_grid_bounds_all(2)) )
984          CALL pmc_recv_from_child( child_id, childs_parent_grid_bounds_all,                        &
985                                    SIZE( childs_parent_grid_bounds_all ), 0, 41, ierr )
[1762]986!
[3888]987!--       Compute size (dimension) of the index_list.
988          index_list_size = 0         
989          DO  n = 1, size_of_childs_parent_grid_bounds_all(2)
990             index_list_size = index_list_size +                                                    &
991                  ( childs_parent_grid_bounds_all(4,n) - childs_parent_grid_bounds_all(3,n) + 1 ) * &
992                  ( childs_parent_grid_bounds_all(2,n) - childs_parent_grid_bounds_all(1,n) + 1 )
[1762]993          ENDDO
994
[3888]995          ALLOCATE( index_list(6,index_list_size) )
[1762]996
[1791]997          nrx = nxr - nxl + 1
[1762]998          nry = nyn - nys + 1
[3888]999          ilist = 0
[1791]1000!
[1933]1001!--       Loop over all children PEs
[3888]1002          DO  n = 1, size_of_childs_parent_grid_bounds_all(2)           !
[1791]1003!
[3888]1004!--          Subspace along y required by actual child PE
1005             DO  jp = childs_parent_grid_bounds_all(3,n), childs_parent_grid_bounds_all(4,n)  ! jp = jps, jpn of child PE# n
[1791]1006!
[3888]1007!--             Subspace along x required by actual child PE
1008                DO  ip = childs_parent_grid_bounds_all(1,n), childs_parent_grid_bounds_all(2,n)  ! ip = ipl, ipr of child PE# n
[1791]1009
[3888]1010                   pex = ip / nrx
1011                   pey = jp / nry
1012                   pe_indices_2d(1) = pex
1013                   pe_indices_2d(2) = pey
1014                   CALL MPI_CART_RANK( comm2d, pe_indices_2d, parent_pe, ierr )
[1762]1015                 
[3888]1016                   ilist = ilist + 1
[1791]1017!
[3946]1018!--                First index in parent array  ! TO_DO: Klaus, please explain better
[3888]1019                   index_list(1,ilist) = ip - ( pex * nrx ) + 1 + nbgp
[1791]1020!
[3946]1021!--                Second index in parent array  ! TO_DO: Klaus, please explain better
[3888]1022                   index_list(2,ilist) = jp - ( pey * nry ) + 1 + nbgp
[1791]1023!
[3946]1024!--                x index of child's parent grid
[3888]1025                   index_list(3,ilist) = ip - childs_parent_grid_bounds_all(1,n) + 1
[1791]1026!
[3888]1027!--                y index of child's parent grid
1028                   index_list(4,ilist) = jp - childs_parent_grid_bounds_all(3,n) + 1
[1791]1029!
[1933]1030!--                PE number of child
[3888]1031                   index_list(5,ilist) = n - 1
[1791]1032!
[1933]1033!--                PE number of parent
[3888]1034                   index_list(6,ilist) = parent_pe
[1791]1035
[1762]1036                ENDDO
1037             ENDDO
1038          ENDDO
[1791]1039!
1040!--       TO_DO: Klaus: comment what is done here
[3888]1041          CALL pmc_s_set_2d_index_list( child_id, index_list(:,1:ilist) )
[1791]1042
[1762]1043       ELSE
[1791]1044!
[1797]1045!--       TO_DO: Klaus: comment why this dummy allocation is required
[1791]1046          ALLOCATE( index_list(6,1) )
[1933]1047          CALL pmc_s_set_2d_index_list( child_id, index_list )
[1762]1048       ENDIF
1049
1050       DEALLOCATE(index_list)
1051
1052     END SUBROUTINE pmci_create_index_list
1053
[3792]1054
1055
[3984]1056     SUBROUTINE pmci_set_child_edge_coords
[2801]1057        IMPLICIT  NONE
1058
[3984]1059        INTEGER(iwp) ::  nbgp_lpm = 1  !< Number of ghost-point layers used for lpm (Klaus, is this correct?)
[2801]1060
[3883]1061       
1062        nbgp_lpm = MIN( nbgp_lpm, nbgp )
[2801]1063
[3883]1064        childgrid(m)%nx = nx_child
1065        childgrid(m)%ny = ny_child
1066        childgrid(m)%nz = nz_child
1067        childgrid(m)%dx = dx_child
1068        childgrid(m)%dy = dy_child
1069        childgrid(m)%dz = dz_child
[2801]1070
[3883]1071        childgrid(m)%lx_coord   = child_coord_x(0)
1072        childgrid(m)%lx_coord_b = child_coord_x(-nbgp_lpm)
1073        childgrid(m)%rx_coord   = child_coord_x(nx_child) + dx_child
1074        childgrid(m)%rx_coord_b = child_coord_x(nx_child+nbgp_lpm) + dx_child
1075        childgrid(m)%sy_coord   = child_coord_y(0)
1076        childgrid(m)%sy_coord_b = child_coord_y(-nbgp_lpm)
1077        childgrid(m)%ny_coord   = child_coord_y(ny_child) + dy_child
1078        childgrid(m)%ny_coord_b = child_coord_y(ny_child+nbgp_lpm) + dy_child
[3984]1079        childgrid(m)%uz_coord   = child_grid_info(2)
1080        childgrid(m)%uz_coord_b = child_grid_info(1)
[2801]1081
[3984]1082     END SUBROUTINE pmci_set_child_edge_coords
[2801]1083
[1764]1084#endif
[1933]1085 END SUBROUTINE pmci_setup_parent
[1762]1086
1087
1088
[1933]1089 SUBROUTINE pmci_setup_child
[1764]1090
1091#if defined( __parallel )
[1762]1092    IMPLICIT NONE
[1764]1093
[3883]1094    INTEGER(iwp) ::  ierr                          !< MPI error code
1095    INTEGER(iwp) ::  lb                            !< Running index for aerosol size bins
1096    INTEGER(iwp) ::  lc                            !< Running index for aerosol mass bins
[3984]1097    INTEGER(iwp) ::  lg                            !< Running index for SALSA gases
[3883]1098    INTEGER(iwp) ::  n                             !< Running index for number of chemical species
1099    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for sending the child-grid dimensions to parent
[3864]1100
[3883]1101    REAL(wp), DIMENSION(5) ::  child_grid_info     !< Array for sending the child-grid spacings etc to parent
1102         
[3984]1103    CHARACTER( LEN=da_namelen ) ::  myname         !< Name of the variable to be coupled
1104    CHARACTER(LEN=5) ::  salsa_char                !< Name extension for the variable name in case of SALSA variable
[3883]1105   
[1791]1106!
[2801]1107!-- Child setup
1108!-- Root model does not have a parent and is not a child, therefore no child setup on root model
[1764]1109    IF ( .NOT. pmc_is_rootmodel() )  THEN
[3883]1110!
[3984]1111!--    KLaus, add a description here what pmc_childinit does       
[1933]1112       CALL pmc_childinit
[1779]1113!
[3792]1114!--    The arrays, which actually will be exchanged between child and parent
1115!--    are defined Here AND ONLY HERE.
[1892]1116!--    If a variable is removed, it only has to be removed from here.
1117!--    Please check, if the arrays are in the list of POSSIBLE exchange arrays
[1779]1118!--    in subroutines:
[1933]1119!--    pmci_set_array_pointer (for parent arrays)
[3883]1120!--    pmci_create_childs_parent_grid_arrays (for child's parent-grid arrays)
[3984]1121       CALL pmc_set_dataarray_name( 'parent', 'u', 'child', 'u', ierr )
1122       CALL pmc_set_dataarray_name( 'parent', 'v', 'child', 'v', ierr )
1123       CALL pmc_set_dataarray_name( 'parent', 'w', 'child', 'w', ierr )
[2938]1124!
1125!--    Set data array name for TKE. Please note, nesting of TKE is actually
1126!--    only done if both parent and child are in LES or in RANS mode. Due to
1127!--    design of model coupler, however, data array names must be already
1128!--    available at this point.
[3984]1129       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                              &
1130            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                               &
[3932]1131               .NOT. constant_diffusion ) )  THEN
[3984]1132          CALL pmc_set_dataarray_name( 'parent', 'e', 'child', 'e', ierr )
[3932]1133       ENDIF
[2938]1134!
1135!--    Nesting of dissipation rate only if both parent and child are in RANS
[3946]1136!--    mode and TKE-epsilon closure is applied. Please see also comment for TKE
[2938]1137!--    above.
[3932]1138       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
[3984]1139          CALL pmc_set_dataarray_name( 'parent', 'diss', 'child', 'diss', ierr )
[3932]1140       ENDIF
[2174]1141
[1892]1142       IF ( .NOT. neutral )  THEN
[3984]1143          CALL pmc_set_dataarray_name( 'parent', 'pt' ,'child', 'pt', ierr )
[1892]1144       ENDIF
[2174]1145
[2003]1146       IF ( humidity )  THEN
[2174]1147
[3984]1148          CALL pmc_set_dataarray_name( 'parent', 'q', 'child', 'q', ierr )
[2174]1149
[3274]1150          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
[3984]1151            CALL pmc_set_dataarray_name( 'parent', 'qc', 'child', 'qc', ierr ) 
1152            CALL pmc_set_dataarray_name( 'parent', 'nc', 'child', 'nc', ierr ) 
[2292]1153          ENDIF
1154
[3274]1155          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
[3984]1156             CALL pmc_set_dataarray_name( 'parent', 'qr', 'child', 'qr', ierr )
1157             CALL pmc_set_dataarray_name( 'parent', 'nr', 'child', 'nr', ierr ) 
[2174]1158          ENDIF
1159     
[1762]1160       ENDIF
[2174]1161
[2003]1162       IF ( passive_scalar )  THEN
[3984]1163          CALL pmc_set_dataarray_name( 'parent', 's', 'child', 's', ierr )
[2003]1164       ENDIF
[1762]1165
[3883]1166       IF ( particle_advection )  THEN
[3984]1167          CALL pmc_set_dataarray_name( 'parent', 'nr_part', 'child', 'nr_part', ierr )
1168          CALL pmc_set_dataarray_name( 'parent', 'part_adr', 'child', 'part_adr', ierr )
[2801]1169       ENDIF
1170       
[4273]1171       IF ( air_chemistry  .AND.  nesting_chem )  THEN
[3883]1172          DO n = 1, nspec
[3984]1173             CALL pmc_set_dataarray_name( 'parent', 'chem_' // TRIM( chem_species(n)%name ),        &
1174                                          'child',  'chem_' // TRIM( chem_species(n)%name ), ierr )
[2773]1175          ENDDO
1176       ENDIF
1177
[4273]1178       IF ( salsa  .AND.  nesting_salsa )  THEN
[3883]1179          DO  lb = 1, nbins_aerosol
1180             WRITE(salsa_char,'(i0)') lb
[3984]1181             CALL pmc_set_dataarray_name( 'parent', 'an_' // TRIM( salsa_char ),                    &
1182                                          'child',  'an_' // TRIM( salsa_char ), ierr )
[3864]1183          ENDDO
[3883]1184          DO  lc = 1, nbins_aerosol * ncomponents_mass
1185             WRITE(salsa_char,'(i0)') lc
[3984]1186             CALL pmc_set_dataarray_name( 'parent', 'am_' // TRIM( salsa_char ),                    &
1187                                          'child',  'am_' // TRIM( salsa_char ), ierr )
[3864]1188          ENDDO
1189          IF ( .NOT. salsa_gases_from_chem )  THEN
[3883]1190             DO  lg = 1, ngases_salsa
1191                WRITE(salsa_char,'(i0)') lg
[3984]1192                CALL pmc_set_dataarray_name( 'parent', 'sg_' // TRIM( salsa_char ),                 &
1193                                             'child',  'sg_' // TRIM( salsa_char ), ierr )
[3864]1194             ENDDO
1195          ENDIF
1196       ENDIF
1197
[1764]1198       CALL pmc_set_dataarray_name( lastentry = .TRUE. )
[1762]1199!
[1933]1200!--    Send grid to parent
[3883]1201       child_grid_dim(1)  = nx
1202       child_grid_dim(2)  = ny
1203       child_grid_dim(3)  = nz
1204       child_grid_info(1) = zw(nzt+1)
1205       child_grid_info(2) = zw(nzt)
1206       child_grid_info(3) = dx
1207       child_grid_info(4) = dy
1208       child_grid_info(5) = dz(1)
[1762]1209
1210       IF ( myid == 0 )  THEN
[1791]1211
[3883]1212          CALL pmc_send_to_parent( child_grid_dim, SIZE( child_grid_dim ), 0, 123, ierr )
1213          CALL pmc_send_to_parent( child_grid_info, SIZE( child_grid_info ), 0, 124, ierr )
[1933]1214          CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
1215          CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
[2938]1216
1217          CALL pmc_recv_from_parent( rans_mode_parent, 1, 0, 19, ierr )
[1762]1218!
[3948]1219!--       Receive parent-grid information.
[2801]1220          CALL pmc_recv_from_parent( parent_grid_info_real,                    &
[2602]1221                                     SIZE(parent_grid_info_real), 0, 21, ierr )
1222          CALL pmc_recv_from_parent( parent_grid_info_int,  3, 0, 22, ierr )
[3792]1223
[1762]1224       ENDIF
1225
[2801]1226       CALL MPI_BCAST( parent_grid_info_real, SIZE(parent_grid_info_real),     &
[1779]1227                       MPI_REAL, 0, comm2d, ierr )
[2602]1228       CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr )
[1762]1229
[3883]1230       pg%dx = parent_grid_info_real(3)
1231       pg%dy = parent_grid_info_real(4)
1232       pg%dz = parent_grid_info_real(7)
1233       pg%nx = parent_grid_info_int(1)
1234       pg%ny = parent_grid_info_int(2)
1235       pg%nz = parent_grid_info_int(3)
[1762]1236!
[3948]1237!--    Allocate 1-D arrays for parent-grid coordinates and grid-spacings in the z-direction
[3883]1238       ALLOCATE( pg%coord_x(-nbgp:pg%nx+nbgp) )
1239       ALLOCATE( pg%coord_y(-nbgp:pg%ny+nbgp) )
1240       ALLOCATE( pg%dzu(1:pg%nz+1) )
1241       ALLOCATE( pg%dzw(1:pg%nz+1) )
1242       ALLOCATE( pg%zu(0:pg%nz+1) )
1243       ALLOCATE( pg%zw(0:pg%nz+1) )
[1762]1244!
[3948]1245!--    Get parent-grid coordinates and grid-spacings in the z-direction from the parent
[1791]1246       IF ( myid == 0)  THEN
[3883]1247          CALL pmc_recv_from_parent( pg%coord_x, pg%nx+1+2*nbgp, 0, 24, ierr )
1248          CALL pmc_recv_from_parent( pg%coord_y, pg%ny+1+2*nbgp, 0, 25, ierr )
1249          CALL pmc_recv_from_parent( pg%dzu, pg%nz+1, 0, 26, ierr )
1250          CALL pmc_recv_from_parent( pg%dzw, pg%nz+1, 0, 27, ierr )
1251          CALL pmc_recv_from_parent( pg%zu, pg%nz+2, 0, 28, ierr )
1252          CALL pmc_recv_from_parent( pg%zw, pg%nz+2, 0, 29, ierr )
[1762]1253       ENDIF
1254!
[1791]1255!--    Broadcast this information
[3883]1256       CALL MPI_BCAST( pg%coord_x, pg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1257       CALL MPI_BCAST( pg%coord_y, pg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1258       CALL MPI_BCAST( pg%dzu, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
1259       CALL MPI_BCAST( pg%dzw, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
1260       CALL MPI_BCAST( pg%zu, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
1261       CALL MPI_BCAST( pg%zw, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
[3948]1262       CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )       
[1791]1263!
[3984]1264!--    Find the index bounds for the nest domain in the parent-grid index space
1265       CALL pmci_map_child_grid_to_parent_grid
[1797]1266!
1267!--    TO_DO: Klaus give a comment what is happening here
[1762]1268       CALL pmc_c_get_2d_index_list
1269!
[1933]1270!--    Include couple arrays into child content
1271!--    TO_DO: Klaus: better explain the above comment (what is child content?)
[1779]1272       CALL  pmc_c_clear_next_array_list
[2773]1273
[3864]1274       n  = 1
[3883]1275       lb = 1
1276       lc = 1
1277       lg = 1
[3864]1278
[1762]1279       DO  WHILE ( pmc_c_getnextarray( myname ) )
[3932]1280!
[3984]1281!--       Note that pg%nz is not the original nz of parent, but the highest
[2773]1282!--       parent-grid level needed for nesting.
[3984]1283!--       Note that in case of chemical species or SALSA variables an additional
1284!--       parameter needs to be passed. The parameter is required to set the pointer
1285!--       correctlyto the chemical-species or SALSA data structure. Hence, first check if
1286!--       the current variable is a chemical species or a SALSA variable. If so, pass
1287!--       index id of respective sub-variable (species or bin) and increment this subsequently.
[2773]1288          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
[3883]1289             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, n )
[3864]1290             n = n + 1   
1291          ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
[3883]1292             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lb )
1293             lb = lb + 1
[3864]1294          ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
[3883]1295             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lc )
1296             lc = lc + 1
1297          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.  .NOT.  salsa_gases_from_chem )  THEN
1298             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lg )
1299             lg = lg + 1
[2773]1300          ELSE
[3883]1301             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz )
[2773]1302          ENDIF
[1762]1303       ENDDO
1304       CALL pmc_c_setind_and_allocmem
1305!
[3792]1306!--    Precompute the index-mapping arrays
1307       CALL pmci_define_index_mapping
[1762]1308!
[3792]1309!--    Check that the child and parent grid lines do match
[4010]1310       CALL pmci_check_grid_matching
1311!       
1312!--    Compute surface areas of the nest-boundary faces
1313       CALL pmci_compute_face_areas
1314       
[1764]1315    ENDIF
[1762]1316
1317 CONTAINS
1318
[2602]1319
[3984]1320    SUBROUTINE pmci_map_child_grid_to_parent_grid
[1797]1321!
[3984]1322!--    Determine index bounds of interpolation/anterpolation area in the parent-grid index space
[1791]1323       IMPLICIT NONE
[1764]1324
[3948]1325       INTEGER(iwp), DIMENSION(5,numprocs) ::  parent_bound_all     !< Transfer array for parent-grid index bounds
[3883]1326
[3819]1327       INTEGER(iwp), DIMENSION(4)          ::  parent_bound_global  !< Transfer array for global parent-grid index bounds
[3984]1328       INTEGER(iwp), DIMENSION(2)          ::  size_of_array        !< For sending the dimensions of parent_bound_all to parent
[3883]1329
[3984]1330       INTEGER(iwp) ::  ip      !< Running parent-grid index in the x-direction
1331       INTEGER(iwp) ::  iauxl   !< Offset between the index bound ipl and the auxiliary index bound ipla
1332       INTEGER(iwp) ::  iauxr   !< Offset between the index bound ipr and the auxiliary index bound ipra
1333       INTEGER(iwp) ::  ijaux   !< Temporary variable for receiving the index bound from the neighbouring subdomain
1334       INTEGER(iwp) ::  jp      !< Running parent-grid index in the y-direction
1335       INTEGER(iwp) ::  jauxs   !< Offset between the index bound jps and the auxiliary index bound jpsa
1336       INTEGER(iwp) ::  jauxn   !< Offset between the index bound jpn and the auxiliary index bound jpna
[3883]1337
[4260]1338       REAL(wp) ::  tolex       !< Tolerance for grid-line matching in x-direction   
1339       REAL(wp) ::  toley       !< Tolerance for grid-line matching in y-direction   
[3681]1340       REAL(wp) ::  xexl        !< Parent-grid array exceedance behind the left edge of the child PE subdomain
1341       REAL(wp) ::  xexr        !< Parent-grid array exceedance behind the right edge of the child PE subdomain
1342       REAL(wp) ::  yexs        !< Parent-grid array exceedance behind the south edge of the child PE subdomain
1343       REAL(wp) ::  yexn        !< Parent-grid array exceedance behind the north edge of the child PE subdomain
[3984]1344       REAL(wp) ::  xpl         !< Requested left-edge x-coordinate of the parent-grid array domain (at the internal boundaries
1345                                !< the real edge may differ from this in some cases as explained in the comment block below) 
1346       REAL(wp) ::  xpr         !< Requested right-edge x-coordinate of the parent-grid array domain (at the internal boundaries
1347                                !< the real edge may differ from this in some cases as explained in the comment block below)
1348       REAL(wp) ::  yps         !< Requested south-edge y-coordinate of the parent-grid array domain (at the internal boundaries
1349                                !< the real edge may differ from this in some cases as explained in the comment block below)
1350       REAL(wp) ::  ypn         !< Requested south-edge y-coordinate of the parent-grid array domain (at the internal boundaries
1351                                !< the real edge may differ from this in some cases as explained in the comment block below)
[1762]1352
1353!
[3984]1354!--    Determine the index limits for the child's parent-grid arrays (such as uc for example).
1355!--    Note that at the outer edges of the child domain (nest boundaries) these arrays exceed
1356!--    the boundary by two parent-grid cells. At the internal boundaries, there are no
1357!--    exceedances and thus no overlaps with the neighbouring subdomain. If at least half
1358!--    of the parent-grid cell is within the current child sub-domain, then it is included
1359!--    in the current sub-domain's parent-grid array. Else the parent-grid cell is
1360!--    included in the neighbouring subdomain's parent-grid array, or not included at all if
1361!--    we are at the outer edge of the child domain. This may occur especially when a large
[4260]1362!--    grid-spacing ratio is used.
[1791]1363!
[4260]1364!--    Tolerances for grid-line matching.
1365       tolex = tolefac * dx
1366       toley = tolefac * dy
1367!
[3984]1368!--    Left boundary.
1369!--    Extension by two parent-grid cells behind the boundary, see the comment block above.   
[3883]1370       IF ( bc_dirichlet_l )  THEN
[3984]1371          xexl  = 2.0_wp * pg%dx
[3681]1372          iauxl = 0
1373       ELSE
1374          xexl  = 0.0_wp
1375          iauxl = 1
1376       ENDIF
[3984]1377       xpl     = coord_x(nxl) - xexl
1378       DO  ip = 0, pg%nx
[4260]1379          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx >= xpl - tolex )  THEN
[3984]1380             ipl = MAX( 0, ip )
[1791]1381             EXIT
1382          ENDIF
1383       ENDDO
1384!
[3984]1385!--    Right boundary.
1386!--    Extension by two parent-grid cells behind the boundary, see the comment block above.       
[3883]1387       IF ( bc_dirichlet_r )  THEN
[3984]1388          xexr  = 2.0_wp * pg%dx
[3681]1389          iauxr = 0 
1390       ELSE
1391          xexr  = 0.0_wp
1392          iauxr = 1 
1393       ENDIF
[3984]1394       xpr  = coord_x(nxr+1) + xexr
1395       DO  ip = pg%nx, 0 , -1
[4260]1396          IF ( pg%coord_x(ip) + 0.5_wp * pg%dx <= xpr + tolex )  THEN
[3984]1397             ipr = MIN( pg%nx, MAX( ipl, ip ) )
[1791]1398             EXIT
1399          ENDIF
1400       ENDDO
1401!
[3984]1402!--    South boundary.
1403!--    Extension by two parent-grid cells behind the boundary, see the comment block above.   
[3883]1404       IF ( bc_dirichlet_s )  THEN
[3984]1405          yexs  = 2.0_wp * pg%dy
[3681]1406          jauxs = 0 
1407       ELSE
1408          yexs  = 0.0_wp
1409          jauxs = 1 
1410       ENDIF
[3984]1411       yps  = coord_y(nys) - yexs
1412       DO  jp = 0, pg%ny
[4260]1413          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy >= yps - toley )  THEN             
[3984]1414             jps = MAX( 0, jp )
[1791]1415             EXIT
1416          ENDIF
1417       ENDDO
1418!
[3984]1419!--    North boundary.
1420!--    Extension by two parent-grid cells behind the boundary, see the comment block above. 
[3681]1421       IF  ( bc_dirichlet_n )  THEN
[3984]1422          yexn  = 2.0_wp * pg%dy
[3681]1423          jauxn = 0
1424       ELSE
1425          yexn  = 0.0_wp
1426          jauxn = 1
1427       ENDIF
[3984]1428       ypn  = coord_y(nyn+1) + yexn
1429       DO  jp = pg%ny, 0 , -1
[4260]1430          IF ( pg%coord_y(jp) + 0.5_wp * pg%dy <= ypn + toley )  THEN
[3984]1431             jpn = MIN( pg%ny, MAX( jps, jp ) )
[3484]1432             EXIT
1433          ENDIF
1434       ENDDO
1435!
[3883]1436!--    Make sure that the indexing is contiguous (no gaps, no overlaps).
1437!--    This is a safety measure mainly for cases with high grid-spacing
1438!--    ratio and narrow child subdomains.
[3984]1439       IF ( pdims(1) > 1 )  THEN
1440          IF ( nxl == 0 )  THEN
1441             CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
1442          ELSE IF ( nxr == nx )  THEN
1443             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr )
1444             ipl = ijaux + 1
1445          ELSE
1446             CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
1447             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 
1448             ipl = ijaux + 1
1449          ENDIF
[3484]1450       ENDIF
[3984]1451
1452       IF ( pdims(2) > 1 )  THEN
1453          IF ( nys == 0 )  THEN
1454             CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
1455          ELSE IF ( nyn == ny )  THEN
1456             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr )
1457             jps = ijaux + 1
1458          ELSE
1459             CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
1460             CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 
1461             jps = ijaux + 1
1462          ENDIF
[3484]1463       ENDIF
[3984]1464         
1465       WRITE(9,"('pmci_map_child_grid_to_parent_grid. Parent-grid array bounds: ',4(i4,2x))")             &
[3883]1466            ipl, ipr, jps, jpn
[3681]1467       FLUSH(9)
1468
[3948]1469       parent_bound(1) = ipl
1470       parent_bound(2) = ipr
1471       parent_bound(3) = jps
1472       parent_bound(4) = jpn
1473       parent_bound(5) = myid
[3484]1474!
[3984]1475!--    The following auxiliary index bounds are used for allocating index mapping and
1476!--    some other auxiliary arrays.
[3883]1477       ipla = ipl - iauxl
1478       ipra = ipr + iauxr
1479       jpsa = jps - jauxs
1480       jpna = jpn + jauxn
[3681]1481!
[3984]1482!--    The index-bounds parent_bound of all subdomains of the current child domain
1483!--    must be sent to the parent in order for the parent to create the index list.
1484!--    For this reason, the parent_bound arrays are packed together in single
1485!--    array parent_bound_all using MPI_GATHER.       
[1791]1486!--    Note that MPI_Gather receives data from all processes in the rank order
[3984]1487!--    This fact is exploited in creating the index list in pmci_create_index_list.
[3948]1488       CALL MPI_GATHER( parent_bound, 5, MPI_INTEGER, parent_bound_all, 5,                          &
[1791]1489                        MPI_INTEGER, 0, comm2d, ierr )
[1762]1490
[1791]1491       IF ( myid == 0 )  THEN
[3948]1492          size_of_array(1) = SIZE( parent_bound_all, 1 )
1493          size_of_array(2) = SIZE( parent_bound_all, 2 )
[1933]1494          CALL pmc_send_to_parent( size_of_array, 2, 0, 40, ierr )
[3948]1495          CALL pmc_send_to_parent( parent_bound_all, SIZE( parent_bound_all ), 0, 41, ierr )
[3819]1496!
1497!--       Determine the global parent-grid index bounds       
[3948]1498          parent_bound_global(1) = MINVAL( parent_bound_all(1,:) )
1499          parent_bound_global(2) = MAXVAL( parent_bound_all(2,:) )
1500          parent_bound_global(3) = MINVAL( parent_bound_all(3,:) )
1501          parent_bound_global(4) = MAXVAL( parent_bound_all(4,:) )
[1791]1502       ENDIF
[3819]1503!
[3984]1504!--    Broadcast the global parent-grid index bounds to all current child processes
[3819]1505       CALL MPI_BCAST( parent_bound_global, 4, MPI_INTEGER, 0, comm2d, ierr )
1506       iplg = parent_bound_global(1)
1507       iprg = parent_bound_global(2)
1508       jpsg = parent_bound_global(3)
1509       jpng = parent_bound_global(4)
[3984]1510       WRITE( 9, "('pmci_map_child_grid_to_parent_grid. Global parent-grid index bounds iplg, iprg, jpsg, jpng: ',4(i4,2x))" ) &
[3883]1511            iplg, iprg, jpsg, jpng
1512       FLUSH( 9 )
[3819]1513       
[3984]1514    END SUBROUTINE pmci_map_child_grid_to_parent_grid
[1762]1515
[3792]1516     
1517     
1518    SUBROUTINE pmci_define_index_mapping
[1762]1519!
[3792]1520!--    Precomputation of the mapping of the child- and parent-grid indices.
[1762]1521
[1791]1522       IMPLICIT NONE
[1762]1523
[3883]1524       INTEGER(iwp) ::  i         !< Child-grid index in the x-direction
1525       INTEGER(iwp) ::  ii        !< Parent-grid index in the x-direction
[3792]1526       INTEGER(iwp) ::  istart    !<
1527       INTEGER(iwp) ::  ir        !<
[3883]1528       INTEGER(iwp) ::  iw        !< Child-grid index limited to -1 <= iw <= nx+1 for wall_flags_0
1529       INTEGER(iwp) ::  j         !< Child-grid index in the y-direction
1530       INTEGER(iwp) ::  jj        !< Parent-grid index in the y-direction
[3792]1531       INTEGER(iwp) ::  jstart    !<
1532       INTEGER(iwp) ::  jr        !<
[3883]1533       INTEGER(iwp) ::  jw        !< Child-grid index limited to -1 <= jw <= ny+1 for wall_flags_0
1534       INTEGER(iwp) ::  k         !< Child-grid index in the z-direction
1535       INTEGER(iwp) ::  kk        !< Parent-grid index in the z-direction
[3792]1536       INTEGER(iwp) ::  kstart    !<
[3883]1537       INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1 for wall_flags_0
[4249]1538
1539       REAL(wp)     ::  tolex     !< Tolerance for grid-line matching in x-direction   
1540       REAL(wp)     ::  toley     !< Tolerance for grid-line matching in y-direction   
1541       REAL(wp)     ::  tolez     !< Tolerance for grid-line matching in z-direction   
1542
[1762]1543!
[4249]1544!--    Grid-line tolerances.
1545       tolex = tolefac * dx
1546       toley = tolefac * dy
1547       tolez = tolefac * dz(1)
1548!
[3681]1549!--    Allocate child-grid work arrays for interpolation.
[3883]1550       igsr = NINT( pg%dx / dx, iwp )
1551       jgsr = NINT( pg%dy / dy, iwp )
1552       kgsr = NINT( pg%dzw(1) / dzw(1), iwp )
[3792]1553       WRITE(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
1554       FLUSH(9)
[3681]1555!       
1556!--    Determine index bounds for the parent-grid work arrays for
1557!--    interpolation and allocate them.
[3792]1558       CALL pmci_allocate_workarrays
[3681]1559!       
1560!--    Define the MPI-datatypes for parent-grid work array
1561!--    exchange between the PE-subdomains.
[3792]1562       CALL pmci_create_workarray_exchange_datatypes
[3681]1563!
1564!--    First determine kcto and kctw which refer to the uppermost
[3984]1565!--    parent-grid levels below the child top-boundary level.
[4260]1566!--    Note that these comparison tests are not round-off-error
1567!--    sensitive and therefore tolerance buffering is not needed here.
[1797]1568       kk = 0
[3883]1569       DO WHILE ( pg%zu(kk) <= zu(nzt) )
[1797]1570          kk = kk + 1
1571       ENDDO
[3484]1572       kcto = kk - 1
[1762]1573
[1797]1574       kk = 0
[3883]1575       DO WHILE ( pg%zw(kk) <= zw(nzt-1) )
[1797]1576          kk = kk + 1
1577       ENDDO
1578       kctw = kk - 1
[1762]1579
[3883]1580       WRITE( 9, "('kcto, kctw = ', 2(i3,2x))" ) kcto, kctw
1581       FLUSH( 9 )
[3947]1582!       
1583!--    In case of two-way coupling, check that the child domain is sufficiently
1584!--    large in terms of the number of parent-grid cells covered. Otherwise
1585!--    anterpolation is not possible.
1586       IF ( nesting_mode == 'two-way')  THEN
1587          CALL pmci_check_child_domain_size
1588       ENDIF
[3792]1589       
[3883]1590       ALLOCATE( iflu(ipla:ipra) )
1591       ALLOCATE( iflo(ipla:ipra) )
1592       ALLOCATE( ifuu(ipla:ipra) )
1593       ALLOCATE( ifuo(ipla:ipra) )
1594       ALLOCATE( jflv(jpsa:jpna) )
1595       ALLOCATE( jflo(jpsa:jpna) )
1596       ALLOCATE( jfuv(jpsa:jpna) )
1597       ALLOCATE( jfuo(jpsa:jpna) )       
1598       ALLOCATE( kflw(0:pg%nz+1) )
1599       ALLOCATE( kflo(0:pg%nz+1) )
1600       ALLOCATE( kfuw(0:pg%nz+1) )
1601       ALLOCATE( kfuo(0:pg%nz+1) )
1602       ALLOCATE( ijkfc_u(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1603       ALLOCATE( ijkfc_v(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1604       ALLOCATE( ijkfc_w(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
1605       ALLOCATE( ijkfc_s(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
[1762]1606
[2997]1607       ijkfc_u = 0
1608       ijkfc_v = 0
1609       ijkfc_w = 0
1610       ijkfc_s = 0
[1762]1611!
[1797]1612!--    i-indices of u for each ii-index value
1613       istart = nxlg
[3883]1614       DO  ii = ipla, ipra
[3484]1615!
[3741]1616!--       The parent and child grid lines do always match in x, hence we
1617!--       use only the local k,j-child-grid plane for the anterpolation.
[3883]1618!--       However, icru still has to be stored separately as these index bounds
1619!--       are passed as arguments to the interpolation and anterpolation
1620!--       subroutines.
[4260]1621!--       Note that this comparison test is round-off-error sensitive
1622!--       and therefore tolerance buffering is needed here.
[1797]1623          i = istart
[4249]1624          DO WHILE ( pg%coord_x(ii) - coord_x(i) > tolex  .AND. i < nxrg )
[3484]1625             i = i + 1
[1797]1626          ENDDO
[3741]1627          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
1628          ifuu(ii) = iflu(ii)
1629          istart   = iflu(ii)
[3681]1630!
1631!--       Print out the index bounds for checking and debugging purposes
[3883]1632          WRITE( 9, "('pmci_define_index_mapping, ii, iflu, ifuu: ', 3(i4,2x))" )                   &
[3484]1633               ii, iflu(ii), ifuu(ii)
[3883]1634          FLUSH( 9 )
[1797]1635       ENDDO
[3883]1636       WRITE( 9, * )
[1762]1637!
[4260]1638!--    i-indices of others for each ii-index value.
1639!--    Note that these comparison tests are not round-off-error
1640!--    sensitive and therefore tolerance buffering is not needed here.
[1797]1641       istart = nxlg
[3883]1642       DO  ii = ipla, ipra
[1797]1643          i = istart
[4249]1644          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.  ( i < nxrg ) )
[2229]1645             i  = i + 1
[1797]1646          ENDDO
1647          iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
[2229]1648          ir = i
[4249]1649          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx < pg%coord_x(ii) + pg%dx )  .AND.  ( i < nxrg+1 ) )
[2229]1650             i  = i + 1
1651             ir = MIN( i, nxrg )
[1797]1652          ENDDO
[2019]1653          ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg )
[1797]1654          istart = iflo(ii)
[3681]1655!
1656!--       Print out the index bounds for checking and debugging purposes
[3883]1657          WRITE( 9, "('pmci_define_index_mapping, ii, iflo, ifuo: ', 3(i4,2x))" )                   &
[3681]1658               ii, iflo(ii), ifuo(ii)
[3883]1659          FLUSH( 9 )
[1797]1660       ENDDO
[3883]1661       WRITE( 9, * )
[1762]1662!
[1797]1663!--    j-indices of v for each jj-index value
1664       jstart = nysg
[3883]1665       DO  jj = jpsa, jpna
[3484]1666!
[3741]1667!--       The parent and child grid lines do always match in y, hence we
1668!--       use only the local k,i-child-grid plane for the anterpolation.
[3883]1669!--       However, jcnv still has to be stored separately as these index bounds
1670!--       are passed as arguments to the interpolation and anterpolation
1671!--       subroutines.
[4260]1672!--       Note that this comparison test is round-off-error sensitive
1673!--       and therefore tolerance buffering is needed here.
[1797]1674          j = jstart
[4249]1675          DO WHILE ( pg%coord_y(jj) - coord_y(j) > toley  .AND. j < nyng )
[3484]1676             j = j + 1
[1797]1677          ENDDO
[3741]1678          jflv(jj) = MIN( MAX( j, nysg ), nyng )
1679          jfuv(jj) = jflv(jj)
1680          jstart   = jflv(jj)
[3681]1681!
1682!--       Print out the index bounds for checking and debugging purposes
[3883]1683          WRITE( 9, "('pmci_define_index_mapping, jj, jflv, jfuv: ', 3(i4,2x))" )                   &
[3484]1684               jj, jflv(jj), jfuv(jj)
[3792]1685          FLUSH(9)
[1797]1686       ENDDO
[3883]1687       WRITE( 9, * )
[1762]1688!
[1797]1689!--    j-indices of others for each jj-index value
[4260]1690!--    Note that these comparison tests are not round-off-error
1691!--    sensitive and therefore tolerance buffering is not needed here.
[1797]1692       jstart = nysg
[3883]1693       DO  jj = jpsa, jpna
[1797]1694          j = jstart
[3883]1695          DO WHILE ( ( coord_y(j) + 0.5_wp * dy < pg%coord_y(jj) ) .AND. ( j < nyng ) )
[2229]1696             j  = j + 1
[1797]1697          ENDDO
1698          jflo(jj) = MIN( MAX( j, nysg ), nyng )
[2229]1699          jr = j
[4249]1700          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy < pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
[2229]1701             j  = j + 1
1702             jr = MIN( j, nyng )
[1797]1703          ENDDO
[2019]1704          jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng )
1705          jstart = jflo(jj)
[3681]1706!
1707!--       Print out the index bounds for checking and debugging purposes
[3883]1708          WRITE( 9, "('pmci_define_index_mapping, jj, jflo, jfuo: ', 3(i4,2x))" )                   &
[3484]1709               jj, jflo(jj), jfuo(jj)
[3883]1710          FLUSH( 9 )
[1797]1711       ENDDO
[3883]1712       WRITE( 9, * )
[1762]1713!
[1797]1714!--    k-indices of w for each kk-index value
[3484]1715!--    Note that anterpolation index limits are needed also for the top boundary
[3741]1716!--    ghost cell level because they are used also in the interpolation.
[1797]1717       kstart  = 0
1718       kflw(0) = 0
1719       kfuw(0) = 0
[3883]1720       DO  kk = 1, pg%nz+1
[3484]1721!
[3741]1722!--       The parent and child grid lines do always match in z, hence we
1723!--       use only the local j,i-child-grid plane for the anterpolation.
[3883]1724!--       However, kctw still has to be stored separately as these index bounds
1725!--       are passed as arguments to the interpolation and anterpolation
1726!--       subroutines.
[4260]1727!--       Note that this comparison test is round-off-error sensitive
1728!--       and therefore tolerance buffering is needed here.
[1797]1729          k = kstart
[4249]1730          DO WHILE ( ( pg%zw(kk) - zw(k) > tolez )  .AND.  ( k < nzt+1 ) )
[1797]1731             k = k + 1
1732          ENDDO
[3741]1733          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
1734          kfuw(kk) = kflw(kk)
1735          kstart   = kflw(kk)
1736!
1737!--       Print out the index bounds for checking and debugging purposes
[3883]1738          WRITE( 9, "('pmci_define_index_mapping, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))" )      &
1739               kk, kflw(kk), kfuw(kk), nzt,  pg%zu(kk), pg%zw(kk)
1740          FLUSH( 9 )
[1797]1741       ENDDO
[3883]1742       WRITE( 9, * )
[1762]1743!
[1797]1744!--    k-indices of others for each kk-index value
1745       kstart  = 0
1746       kflo(0) = 0
1747       kfuo(0) = 0
[3484]1748!
1749!--    Note that anterpolation index limits are needed also for the top boundary
[3741]1750!--    ghost cell level because they are used also in the interpolation.
[4260]1751!--    Note that these comparison tests are not round-off-error
1752!--    sensitive and therefore tolerance buffering is not needed here.
[3883]1753       DO  kk = 1, pg%nz+1
[1797]1754          k = kstart
[3883]1755          DO WHILE ( ( zu(k) < pg%zw(kk-1) )  .AND.  ( k <= nzt ) )
[1797]1756             k = k + 1
1757          ENDDO
1758          kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
[4249]1759          DO WHILE ( ( zu(k) < pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
[1797]1760             k = k + 1
[3484]1761             IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
[1797]1762          ENDDO
[2019]1763          kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 )
[1797]1764          kstart = kflo(kk)
[3681]1765       ENDDO
1766!
[3792]1767!--    Print out the index bounds for checking and debugging purposes
[3883]1768       DO  kk = 1, pg%nz+1
1769          WRITE( 9, "('pmci_define_index_mapping, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))" )      &
1770               kk, kflo(kk), kfuo(kk), nzt,  pg%zu(kk), pg%zw(kk)
1771          FLUSH( 9 )
[1797]1772       ENDDO
[3883]1773       WRITE( 9, * )
[1762]1774!
[3984]1775!--    Precomputation of number of child-grid nodes inside parent-grid cells.
[3484]1776!--    Note that ii, jj, and kk are parent-grid indices.
[3883]1777!--    This information is needed in the anterpolation.
1778!--    The indices for wall_flags_0 (kw,jw,iw) must be limited to the range
1779!--    [-1,...,nx/ny/nzt+1] in order to avoid zero values on the outer ghost nodes.
1780       DO  ii = ipla, ipra
1781          DO  jj = jpsa, jpna
1782             DO  kk = 0, pg%nz+1
[2997]1783!
1784!--             u-component
1785                DO  i = iflu(ii), ifuu(ii)
[3681]1786                   iw = MAX( MIN( i, nx+1 ), -1 )
[2997]1787                   DO  j = jflo(jj), jfuo(jj)
[3681]1788                      jw = MAX( MIN( j, ny+1 ), -1 )
[2997]1789                      DO  k = kflo(kk), kfuo(kk)
[3681]1790                         kw = MIN( k, nzt+1 )               
[3883]1791                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii)                                      &
[3681]1792                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 1 ) )
[2997]1793                      ENDDO
1794                   ENDDO
1795                ENDDO
1796!
1797!--             v-component
1798                DO  i = iflo(ii), ifuo(ii)
[3681]1799                   iw = MAX( MIN( i, nx+1 ), -1 )
[2997]1800                   DO  j = jflv(jj), jfuv(jj)
[3681]1801                      jw = MAX( MIN( j, ny+1 ), -1 )
[2997]1802                      DO  k = kflo(kk), kfuo(kk)
[3681]1803                         kw = MIN( k, nzt+1 )                                       
[3883]1804                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii)                                      &
[3681]1805                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 2 ) )
[2997]1806                      ENDDO
1807                   ENDDO
1808                ENDDO
1809!
1810!--             scalars
1811                DO  i = iflo(ii), ifuo(ii)
[3681]1812                   iw = MAX( MIN( i, nx+1 ), -1 )
[2997]1813                   DO  j = jflo(jj), jfuo(jj)
[3681]1814                      jw = MAX( MIN( j, ny+1 ), -1 )
[2997]1815                      DO  k = kflo(kk), kfuo(kk)
[3681]1816                         kw = MIN( k, nzt+1 )
[3883]1817                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii)                                      &
[3681]1818                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 0 ) )
[2997]1819                      ENDDO
1820                   ENDDO
1821                ENDDO
1822!
1823!--             w-component
1824                DO  i = iflo(ii), ifuo(ii)
[3681]1825                   iw = MAX( MIN( i, nx+1 ), -1 )
[2997]1826                   DO  j = jflo(jj), jfuo(jj)
[3681]1827                      jw = MAX( MIN( j, ny+1 ), -1 )
[2997]1828                      DO  k = kflw(kk), kfuw(kk)
[3681]1829                         kw = MIN( k, nzt+1 )
[3883]1830                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii)                                      &
1831                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 3 ) )
[2997]1832                      ENDDO
1833                   ENDDO
1834                ENDDO
[3681]1835
1836             ENDDO  ! kk       
1837          ENDDO  ! jj
1838       ENDDO  ! ii
[3792]1839
1840    END SUBROUTINE pmci_define_index_mapping
1841
1842
[3947]1843
1844    SUBROUTINE pmci_check_child_domain_size
1845!       
1846!--    Check if the child domain is too small in terms of number of parent-grid cells
1847!--    covered so that anterpolation buffers fill the whole domain so that anterpolation
1848!--    not possible. Also, check that anterpolation_buffer_width is not too large to 
1849!--    prevent anterpolation.
1850       IMPLICIT NONE
1851     
1852!
1853!--    First x-direction
1854       IF ( iplg + 3 + anterpolation_buffer_width > iprg - 3 - anterpolation_buffer_width )  THEN
1855          IF ( iprg - iplg + 1 < 7 )  THEN
1856!
1857!--          Error
1858             WRITE( message_string, * ) 'child domain too narrow for anterpolation in x-direction'
[3984]1859             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
[3947]1860          ELSE IF ( iprg - iplg + 1 < 11 )  THEN
1861!               
1862!--          Warning
1863             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
[3984]1864             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
[3947]1865             anterpolation_buffer_width = 0
1866          ELSE
1867!               
1868!--          Informative message
1869             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
[3984]1870             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
[3947]1871             anterpolation_buffer_width = 2
1872          ENDIF
1873       ENDIF
1874!
1875!--    Then y-direction         
1876       IF ( jpsg + 3 + anterpolation_buffer_width > jpng - 3 - anterpolation_buffer_width )  THEN
1877          IF ( jpng - jpsg + 1 < 7 )  THEN
1878!
1879!--          Error
1880             WRITE( message_string, * ) 'child domain too narrow for anterpolation in y-direction'
[3984]1881             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
[3947]1882          ELSE IF ( jpng - jpsg + 1 < 11 )  THEN
1883!               
1884!--          Warning
1885             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
[3984]1886             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
[3947]1887             anterpolation_buffer_width = 0
1888          ELSE
1889!               
1890!--          Informative message
1891             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
[3984]1892             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
[3947]1893             anterpolation_buffer_width = 2
1894          ENDIF
1895       ENDIF
1896!
1897!--    Finally z-direction               
1898       IF ( kctw - 1 - anterpolation_buffer_width < 1 )  THEN
1899          IF ( kctw - 1 < 1 )  THEN
1900!
1901!--          Error
1902             WRITE( message_string, * ) 'child domain too shallow for anterpolation in z-direction'
[3984]1903             CALL message( 'pmci_check_child_domain_size', 'PA0652', 3, 2, 0, 6, 0 )
[3947]1904          ELSE IF ( kctw - 3 < 1 )  THEN
1905!               
1906!--          Warning
1907             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to 0'
[3984]1908             CALL message( 'pmci_check_child_domain_size', 'PA0653', 0, 1, 0, 6, 0 )
[3947]1909             anterpolation_buffer_width = 0
1910          ELSE
1911!               
1912!--          Informative message
1913             WRITE( message_string, * ) 'anterpolation_buffer_width value too high, reset to default value 2'
[3984]1914             CALL message( 'pmci_check_child_domain_size', 'PA0654', 0, 0, 0, 6, 0 )
[3947]1915             anterpolation_buffer_width = 2 
1916          ENDIF
1917       ENDIF
1918
1919    END SUBROUTINE pmci_check_child_domain_size
1920
1921
1922   
[3792]1923    SUBROUTINE pmci_allocate_workarrays
[1882]1924!
[3792]1925!--    Allocate parent-grid work-arrays for interpolation
1926       IMPLICIT NONE
[1791]1927
[3792]1928!
1929!--    Determine and store the PE-subdomain dependent index bounds
[3883]1930       IF ( bc_dirichlet_l )  THEN
1931          iplw = ipl + 1
[3792]1932       ELSE
[3883]1933          iplw = ipl - 1
[3792]1934       ENDIF
[1762]1935
[3883]1936       IF ( bc_dirichlet_r )  THEN
1937          iprw = ipr - 1
[3792]1938       ELSE
[3883]1939          iprw = ipr + 1
[3792]1940       ENDIF
[1762]1941
[3883]1942       IF ( bc_dirichlet_s )  THEN
1943          jpsw = jps + 1
[3792]1944       ELSE
[3883]1945          jpsw = jps - 1
[3792]1946       ENDIF
[3708]1947
[3883]1948       IF ( bc_dirichlet_n )  THEN
1949          jpnw = jpn - 1
[3792]1950       ELSE
[3883]1951          jpnw = jpn + 1
[3792]1952       ENDIF
1953!
1954!--    Left and right boundaries.
[3883]1955       ALLOCATE( workarr_lr(0:pg%nz+1,jpsw:jpnw,0:2) )
[3792]1956!
1957!--    South and north boundaries.
[3883]1958       ALLOCATE( workarr_sn(0:pg%nz+1,0:2,iplw:iprw) )
[3792]1959!
1960!--    Top boundary.
[3883]1961       ALLOCATE( workarr_t(0:2,jpsw:jpnw,iplw:iprw) )
[3708]1962
[3792]1963    END SUBROUTINE pmci_allocate_workarrays
1964
1965
1966
1967    SUBROUTINE pmci_create_workarray_exchange_datatypes
1968!
[3883]1969!--    Define specific MPI types for workarr-exchange.
[3792]1970       IMPLICIT NONE
1971
1972!
1973!--    For the left and right boundaries
[3883]1974       CALL MPI_TYPE_VECTOR( 3, pg%nz+2, (jpnw-jpsw+1)*(pg%nz+2), MPI_REAL,                         &
1975            workarr_lr_exchange_type, ierr )
1976       CALL MPI_TYPE_COMMIT( workarr_lr_exchange_type, ierr )
[3792]1977!
1978!--    For the south and north boundaries
[3883]1979       CALL MPI_TYPE_VECTOR( 1, 3*(pg%nz+2), 3*(pg%nz+2), MPI_REAL,                                 &
1980            workarr_sn_exchange_type, ierr )
1981       CALL MPI_TYPE_COMMIT( workarr_sn_exchange_type, ierr )
[3792]1982!
1983!--    For the top-boundary x-slices
[3883]1984       CALL MPI_TYPE_VECTOR( iprw-iplw+1, 3, 3*(jpnw-jpsw+1), MPI_REAL,                             &
1985            workarr_t_exchange_type_x, ierr )
1986       CALL MPI_TYPE_COMMIT( workarr_t_exchange_type_x, ierr )
[3792]1987!
1988!--    For the top-boundary y-slices
[3883]1989       CALL MPI_TYPE_VECTOR( 1, 3*(jpnw-jpsw+1), 3*(jpnw-jpsw+1), MPI_REAL,                         &
1990            workarr_t_exchange_type_y, ierr )
1991       CALL MPI_TYPE_COMMIT( workarr_t_exchange_type_y, ierr )
[3792]1992       
1993    END SUBROUTINE pmci_create_workarray_exchange_datatypes
1994
1995
1996
[3708]1997    SUBROUTINE pmci_check_grid_matching
1998!
[3741]1999!--    Check that the grid lines of child and parent do match.
2000!--    Also check that the child subdomain width is not smaller than
2001!--    the parent grid spacing in the respective direction.
[3708]2002       IMPLICIT NONE
2003
[3883]2004       INTEGER(iwp) ::  non_matching_height = 0              !< Flag for non-matching child-domain height
2005       INTEGER(iwp) ::  non_matching_lower_left_corner = 0   !< Flag for non-matching lower left corner
2006       INTEGER(iwp) ::  non_matching_upper_right_corner = 0  !< Flag for non-matching upper right corner
2007       INTEGER(iwp) ::  non_int_gsr_x = 0                    !< Flag for non-integer grid-spacing ration in x-direction
2008       INTEGER(iwp) ::  non_int_gsr_y = 0                    !< Flag for non-integer grid-spacing ration in y-direction
2009       INTEGER(iwp) ::  non_int_gsr_z = 0                    !< Flag for non-integer grid-spacing ration in z-direction
2010       INTEGER(iwp) ::  too_narrow_pesd_x = 0                !< Flag for too narrow pe-subdomain in x-direction
2011       INTEGER(iwp) ::  too_narrow_pesd_y = 0                !< Flag for too narrow pe-subdomain in y-direction
[4249]2012                                                                                                                 
[3883]2013       REAL(wp) ::  child_ngp_x_l                            !< Number of gridpoints in child subdomain in x-direction
2014                                                             !< converted to REAL(wp)
2015       REAL(wp) ::  child_ngp_y_l                            !< Number of gridpoints in child subdomain in y-direction
2016                                                             !< converted to REAL(wp)
[4249]2017       REAL(wp) ::  gridline_mismatch_x                      !< Mismatch between the parent and child gridlines in the x-direction
2018       REAL(wp) ::  gridline_mismatch_y                      !< Mismatch between the parent and child gridlines in the y-direction
2019       REAL(wp) ::  gsr_mismatch_x                           !< Deviation of the grid-spacing ratio from the nearest integer value, the x-direction
2020       REAL(wp) ::  gsr_mismatch_y                           !< Deviation of the grid-spacing ratio from the nearest integer value, the y-direction
2021       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
2022       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
[3883]2023       REAL(wp) ::  tolex                                    !< Tolerance for grid-line matching in x-direction
2024       REAL(wp) ::  toley                                    !< Tolerance for grid-line matching in y-direction
2025       REAL(wp) ::  tolez                                    !< Tolerance for grid-line matching in z-direction
[3708]2026
[3883]2027       
2028       IF ( myid == 0 )  THEN
[3708]2029
2030          tolex = tolefac * dx
2031          toley = tolefac * dy
[3883]2032          tolez = tolefac * dz(1)
[3708]2033!
[3883]2034!--       First check that the child domain lower left corner matches the parent grid lines.
[4249]2035          gridline_mismatch_x = ABS( NINT( lower_left_coord_x / pg%dx ) * pg%dx - lower_left_coord_x )
2036          gridline_mismatch_y = ABS( NINT( lower_left_coord_y / pg%dy ) * pg%dy - lower_left_coord_y )
2037          IF ( gridline_mismatch_x > tolex ) non_matching_lower_left_corner = 1
2038          IF ( gridline_mismatch_y > toley ) non_matching_lower_left_corner = 1
[3708]2039!
[3883]2040!--       Then check that the child doman upper right corner matches the parent grid lines.
2041          upper_right_coord_x = lower_left_coord_x + ( nx + 1 ) * dx
2042          upper_right_coord_y = lower_left_coord_y + ( ny + 1 ) * dy
[4249]2043          gridline_mismatch_x = ABS( NINT( upper_right_coord_x / pg%dx ) * pg%dx - upper_right_coord_x )
2044          gridline_mismatch_y = ABS( NINT( upper_right_coord_y / pg%dy ) * pg%dy - upper_right_coord_y )
2045          IF ( gridline_mismatch_x > tolex ) non_matching_upper_right_corner = 1
2046          IF ( gridline_mismatch_y > toley ) non_matching_upper_right_corner = 1
[3883]2047!
2048!--       Also check that the cild domain height matches the parent grid lines.
2049          IF ( MOD( zw(nzt), pg%dz ) > tolez ) non_matching_height = 1
2050!
2051!--       Check that the grid-spacing ratios in each direction are integer valued.   
[4249]2052          gsr_mismatch_x = ABS( NINT( pg%dx / dx ) * dx - pg%dx )
2053          gsr_mismatch_y = ABS( NINT( pg%dy / dy ) * dy - pg%dy )
2054          IF ( gsr_mismatch_x > tolex )  non_int_gsr_x = 1
2055          IF ( gsr_mismatch_y > toley )  non_int_gsr_y = 1
[3883]2056!
2057!--       In the z-direction, all levels need to be checked separately against grid stretching 
2058!--       which is not allowed.
[3708]2059          DO  n = 0, kctw+1
[3883]2060             IF ( ABS( pg%zw(n) - zw(kflw(n)) ) > tolez )  non_int_gsr_z = 1
[3708]2061          ENDDO
2062
[3883]2063          child_ngp_x_l = REAL( nxr - nxl + 1, KIND=wp )
2064          IF ( child_ngp_x_l / REAL( igsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_x = 1
2065          child_ngp_y_l = REAL( nyn - nys + 1, KIND=wp )
2066          IF ( child_ngp_y_l / REAL( jgsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_y = 1
2067         
2068          IF ( non_matching_height > 0 )  THEN
2069             WRITE( message_string, * ) 'nested child domain height must match ',                   &
2070                                        'its parent grid lines'
2071             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2072          ENDIF
[3708]2073
[3883]2074          IF ( non_matching_lower_left_corner > 0 )  THEN
2075             WRITE( message_string, * ) 'nested child domain lower left ',                          &
2076                                        'corner must match its parent grid lines'
2077             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2078          ENDIF
[3708]2079
[3883]2080          IF ( non_matching_upper_right_corner > 0 )  THEN
2081             WRITE( message_string, * ) 'nested child domain upper right ',                         &
2082                                        'corner must match its parent grid lines'
2083             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
2084          ENDIF
[3708]2085
[3883]2086          IF ( non_int_gsr_x > 0 )  THEN
2087             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dx / child dx ) ',     &
2088                                        'must have an integer value'
2089             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2090          ENDIF
[3708]2091
[3883]2092          IF ( non_int_gsr_y > 0 )  THEN
2093             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dy / child dy ) ',     &
2094                                        'must have an integer value'
2095             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2096          ENDIF
[3741]2097
[3883]2098          IF ( non_int_gsr_z > 0 )  THEN
2099             WRITE( message_string, * ) 'nesting grid-spacing ratio ( parent dz / child dz ) ',     &
2100                                        'must have an integer value for each z-level'
2101             CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
2102          ENDIF
2103
2104          IF ( too_narrow_pesd_x > 0 )  THEN
[4249]2105            WRITE( message_string, * ) 'child subdomain width in x-direction must not be ',        &
[3883]2106                                        'smaller than its parent grid dx. Change the PE-grid ',     &
2107                                        'setting (npex, npey) to satisfy this requirement.' 
2108             CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
2109          ENDIF
[3741]2110 
[3883]2111          IF ( too_narrow_pesd_y > 0 )  THEN
2112             WRITE( message_string, * ) 'child subdomain width in y-direction must not be ',        &
2113                                        'smaller than its parent grid dy. Change the PE-grid ',     &
2114                                        'setting (npex, npey) to satisfy this requirement.' 
2115             CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
2116          ENDIF
2117                 
2118       ENDIF  !  ( myid == 0 )
2119       
2120    END SUBROUTINE pmci_check_grid_matching
[4010]2121
2122
2123
2124    SUBROUTINE pmci_compute_face_areas
2125
2126       IMPLICIT NONE
2127       REAL(wp)  :: face_area_local   !< Local (for the current pe) integral face area of the left boundary
2128       REAL(wp)  :: sub_sum           !< Intermediate sum in order to improve the accuracy of the summation
2129
2130       INTEGER(iwp) :: i              !< Running index in the x-direction
2131       INTEGER(iwp) :: j              !< Running index in the y-direction
2132       INTEGER(iwp) :: k              !< Running index in the z-direction
2133       INTEGER(iwp) :: k_wall         !< Local topography top k-index
2134       INTEGER(iwp) :: n              !< Running index over boundary faces
2135
2136       
2137!
2138!--    Sum up the volume flow through the left boundary
2139       face_area(1) = 0.0_wp
2140       face_area_local = 0.0_wp
2141       IF ( bc_dirichlet_l )  THEN
2142          i = 0
2143          DO  j = nys, nyn
2144             sub_sum = 0.0_wp
[4168]2145             k_wall = topo_top_ind(j,i,1)
[4010]2146             DO   k = k_wall + 1, nzt
2147                sub_sum = sub_sum + dzw(k)
2148             ENDDO
2149             face_area_local =  face_area_local + dy * sub_sum
2150          ENDDO
2151       ENDIF
2152       
2153#if defined( __parallel )
2154       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2155       CALL MPI_ALLREDUCE( face_area_local, face_area(1), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2156#else
2157       face_area(1) = face_area_local
2158#endif
2159!
2160!--    Sum up the volume flow through the right boundary
2161       face_area(2) = 0.0_wp
2162       face_area_local = 0.0_wp
2163       IF ( bc_dirichlet_r )  THEN
2164          i = nx
2165          DO  j = nys, nyn
2166             sub_sum = 0.0_wp
[4168]2167             k_wall = topo_top_ind(j,i,1)
[4010]2168             DO   k = k_wall + 1, nzt
2169                sub_sum = sub_sum + dzw(k)
2170             ENDDO
2171             face_area_local =  face_area_local + dy * sub_sum
2172          ENDDO
2173       ENDIF
2174       
2175#if defined( __parallel )
2176       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2177       CALL MPI_ALLREDUCE( face_area_local, face_area(2), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2178#else
2179       face_area(2) = face_area_local
2180#endif
2181!
2182!--    Sum up the volume flow through the south boundary
2183       face_area(3) = 0.0_wp
2184       face_area_local = 0.0_wp
2185       IF ( bc_dirichlet_s )  THEN
2186          j = 0
2187          DO  i = nxl, nxr
2188             sub_sum = 0.0_wp
[4168]2189             k_wall = topo_top_ind(j,i,2)
[4010]2190             DO  k = k_wall + 1, nzt
2191                sub_sum = sub_sum + dzw(k)
2192             ENDDO
2193             face_area_local = face_area_local + dx * sub_sum
2194          ENDDO
2195       ENDIF
2196       
2197#if defined( __parallel )
2198       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2199       CALL MPI_ALLREDUCE( face_area_local, face_area(3), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2200#else
2201       face_area(3) = face_area_local
2202#endif
2203!
2204!--    Sum up the volume flow through the north boundary
2205       face_area(4) = 0.0_wp
2206       face_area_local = 0.0_wp
2207       IF ( bc_dirichlet_n )  THEN
2208          j = ny
2209          DO  i = nxl, nxr
2210             sub_sum = 0.0_wp
[4168]2211             k_wall = topo_top_ind(j,i,2)
[4010]2212             DO  k = k_wall + 1, nzt
2213                sub_sum = sub_sum + dzw(k)
2214             ENDDO
2215             face_area_local = face_area_local + dx * sub_sum
2216          ENDDO
2217       ENDIF
2218       
2219#if defined( __parallel )
2220       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2221       CALL MPI_ALLREDUCE( face_area_local, face_area(4), 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2222#else
2223       face_area(4) = face_area_local
2224#endif
2225!
2226!--    The top face area does not depend on the topography at all.       
2227       face_area(5) = ( nx + 1 ) * ( ny + 1 ) * dx * dy
2228!
2229!--    The 6th element is used for the total area
2230       face_area(6) = 0.0_wp
2231       DO  n = 1, 5
2232          face_area(6) = face_area(6) + face_area(n)
2233       ENDDO
2234
[4249]2235!       write( 9, "(6(e12.5,2x))") ( face_area(n), n = 1, 6 )
2236!       flush( 9 )
[4010]2237       
2238    END SUBROUTINE pmci_compute_face_areas
2239#endif
2240   
[1933]2241 END SUBROUTINE pmci_setup_child
[1762]2242
2243
2244
2245 SUBROUTINE pmci_setup_coordinates
[1764]2246
2247#if defined( __parallel )
[1762]2248    IMPLICIT NONE
[1764]2249
[2801]2250    INTEGER(iwp) ::  i   !<
2251    INTEGER(iwp) ::  j   !<
[1762]2252
2253!
2254!-- Create coordinate arrays.
2255    ALLOCATE( coord_x(-nbgp:nx+nbgp) )
2256    ALLOCATE( coord_y(-nbgp:ny+nbgp) )
2257     
2258    DO  i = -nbgp, nx + nbgp
2259       coord_x(i) = lower_left_coord_x + i * dx
2260    ENDDO
[3484]2261
[1762]2262    DO  j = -nbgp, ny + nbgp
2263       coord_y(j) = lower_left_coord_y + j * dy
2264    ENDDO
[1764]2265
2266#endif
[1762]2267 END SUBROUTINE pmci_setup_coordinates
[3681]2268
[3592]2269!------------------------------------------------------------------------------!
2270! Description:
2271! ------------
2272!> In this subroutine the number of coupled arrays is determined.
2273!------------------------------------------------------------------------------!
2274  SUBROUTINE pmci_num_arrays 
2275               
2276#if defined( __parallel )
2277    USE pmc_general,                                                           &
2278        ONLY:  pmc_max_array
[1762]2279
[3592]2280    IMPLICIT NONE
2281!
2282!-- The number of coupled arrays depends on the model settings. At least
2283!-- 5 arrays need to be coupled (u, v, w, e, diss).  Please note, actually
2284!-- e and diss (TKE and dissipation rate) are only required if RANS-RANS
2285!-- nesting is applied, but memory is allocated nevertheless. This is because
2286!-- the information whether they are needed or not is retrieved at a later
2287!-- point in time. In case e and diss are not needed, they are also not
2288!-- exchanged between parent and child.
2289    pmc_max_array = 5
2290!
2291!-- pt
2292    IF ( .NOT. neutral )  pmc_max_array = pmc_max_array + 1
2293   
2294    IF ( humidity )  THEN
2295!
2296!--    q
2297       pmc_max_array = pmc_max_array + 1
2298!
2299!--    qc, nc
2300       IF ( bulk_cloud_model  .AND.  microphysics_morrison )                   &
2301          pmc_max_array = pmc_max_array + 2
2302!
2303!--    qr, nr
2304       IF ( bulk_cloud_model  .AND.  microphysics_seifert )                    &
2305          pmc_max_array = pmc_max_array + 2
2306    ENDIF
2307!
2308!-- s
2309    IF ( passive_scalar )  pmc_max_array = pmc_max_array + 1
2310!
2311!-- nr_part, part_adr
2312    IF ( particle_advection )  pmc_max_array = pmc_max_array + 2
2313!
2314!-- Chemistry, depends on number of species
[4273]2315    IF ( air_chemistry  .AND.  nesting_chem )  pmc_max_array = pmc_max_array + nspec
[3864]2316!
[3883]2317!-- SALSA, depens on the number aerosol size bins and chemical components +
[3864]2318!-- the number of default gases
[4273]2319    IF ( salsa  .AND.  nesting_salsa )  pmc_max_array = pmc_max_array + nbins_aerosol +            &
2320                                                        nbins_aerosol * ncomponents_mass
[3883]2321    IF ( .NOT. salsa_gases_from_chem )  pmc_max_array = pmc_max_array + ngases_salsa
[3864]2322
[3592]2323#endif
[3792]2324   
[3592]2325 END SUBROUTINE pmci_num_arrays
2326
2327
[3883]2328 SUBROUTINE pmci_set_array_pointer( name, child_id, nz_child, n )
[3948]2329   
[1791]2330    IMPLICIT NONE
[3883]2331   
2332    INTEGER(iwp), INTENT(IN) ::  child_id  !<
2333    INTEGER(iwp), INTENT(IN) ::  nz_child  !<
2334   
2335    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n          !< index of chemical species
2336   
[3792]2337    CHARACTER(LEN=*), INTENT(IN) ::  name             !<
[3948]2338
2339#if defined( __parallel )     
[3883]2340!
2341!-- Local variables:       
2342    INTEGER(iwp) ::  ierr                             !< MPI error code
[1791]2343
[3883]2344    INTEGER(idp), POINTER, DIMENSION(:,:) ::  i_2d    !<
2345       
2346    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d      !<
2347    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d      !<
2348    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d_sec  !<
2349   
[1791]2350
2351    NULLIFY( p_3d )
2352    NULLIFY( p_2d )
[2801]2353    NULLIFY( i_2d )
[1791]2354!
2355!-- List of array names, which can be coupled.
2356!-- In case of 3D please change also the second array for the pointer version
[2938]2357    IF ( TRIM(name) == "u"          )  p_3d => u
2358    IF ( TRIM(name) == "v"          )  p_3d => v
2359    IF ( TRIM(name) == "w"          )  p_3d => w
2360    IF ( TRIM(name) == "e"          )  p_3d => e
2361    IF ( TRIM(name) == "pt"         )  p_3d => pt
2362    IF ( TRIM(name) == "q"          )  p_3d => q
2363    IF ( TRIM(name) == "qc"         )  p_3d => qc
2364    IF ( TRIM(name) == "qr"         )  p_3d => qr
2365    IF ( TRIM(name) == "nr"         )  p_3d => nr
2366    IF ( TRIM(name) == "nc"         )  p_3d => nc
2367    IF ( TRIM(name) == "s"          )  p_3d => s
2368    IF ( TRIM(name) == "diss"       )  p_3d => diss   
[3209]2369    IF ( TRIM(name) == "nr_part"    )  i_2d => nr_part
[2938]2370    IF ( TRIM(name) == "part_adr"   )  i_2d => part_adr
[3864]2371    IF ( INDEX( TRIM(name), "chem_" ) /= 0      )  p_3d => chem_species(n)%conc
2372    IF ( INDEX( TRIM(name), "an_" ) /= 0  )  p_3d => aerosol_number(n)%conc
2373    IF ( INDEX( TRIM(name), "am_" ) /= 0 )  p_3d => aerosol_mass(n)%conc
2374    IF ( INDEX( TRIM(name), "sg_" ) /= 0  .AND.  .NOT. salsa_gases_from_chem ) &
2375       p_3d => salsa_gas(n)%conc
[1791]2376!
[3883]2377!-- Next line is just an example for a 2D array (not active for coupling!)
2378!-- Please note, that z0 has to be declared as TARGET array in modules.f90.
[1791]2379!    IF ( TRIM(name) == "z0" )    p_2d => z0
[2938]2380    IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
2381    IF ( TRIM(name) == "v"    )  p_3d_sec => v_2
2382    IF ( TRIM(name) == "w"    )  p_3d_sec => w_2
2383    IF ( TRIM(name) == "e"    )  p_3d_sec => e_2
2384    IF ( TRIM(name) == "pt"   )  p_3d_sec => pt_2
2385    IF ( TRIM(name) == "q"    )  p_3d_sec => q_2
2386    IF ( TRIM(name) == "qc"   )  p_3d_sec => qc_2
2387    IF ( TRIM(name) == "qr"   )  p_3d_sec => qr_2
2388    IF ( TRIM(name) == "nr"   )  p_3d_sec => nr_2
2389    IF ( TRIM(name) == "nc"   )  p_3d_sec => nc_2
2390    IF ( TRIM(name) == "s"    )  p_3d_sec => s_2
2391    IF ( TRIM(name) == "diss" )  p_3d_sec => diss_2
[4029]2392    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d_sec => spec_conc_2(:,:,:,n)
2393    IF ( INDEX( TRIM(name), "an_" )   /= 0 )  p_3d_sec => nconc_2(:,:,:,n)
2394    IF ( INDEX( TRIM(name), "am_" )   /= 0 )  p_3d_sec => mconc_2(:,:,:,n)
2395    IF ( INDEX( TRIM(name), "sg_" )   /= 0  .AND.  .NOT. salsa_gases_from_chem ) &
2396                                 p_3d_sec => gconc_2(:,:,:,n)
[1791]2397
2398    IF ( ASSOCIATED( p_3d ) )  THEN
[3883]2399       CALL pmc_s_set_dataarray( child_id, p_3d, nz_child, nz, array_2 = p_3d_sec )
2400    ELSEIF  ( ASSOCIATED( p_2d ) )  THEN
[1933]2401       CALL pmc_s_set_dataarray( child_id, p_2d )
[3883]2402    ELSEIF  ( ASSOCIATED( i_2d ) )  THEN
[2801]2403       CALL pmc_s_set_dataarray( child_id, i_2d )
[1791]2404    ELSE
2405!
2406!--    Give only one message for the root domain
[3948]2407       IF ( pmc_is_rootmodel()  .AND.  myid == 0 )  THEN
[3883]2408          message_string = 'pointer for array "' // TRIM( name ) // '" can''t be associated'
[1791]2409          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
2410       ELSE
2411!
2412!--       Avoid others to continue
2413          CALL MPI_BARRIER( comm2d, ierr )
2414       ENDIF
[3883]2415       
2416    ENDIF
2417   
[3948]2418#endif
2419   
[3883]2420 END SUBROUTINE pmci_set_array_pointer
[1791]2421
2422
[3883]2423     
2424 INTEGER FUNCTION get_number_of_children()
[1791]2425
[3883]2426    IMPLICIT NONE
[2967]2427
[3883]2428   
[2967]2429#if defined( __parallel )
[3883]2430    get_number_of_children = SIZE( pmc_parent_for_child ) - 1
[2967]2431#else
[3883]2432    get_number_of_children = 0
[2967]2433#endif
[1791]2434
[3883]2435    RETURN
[2967]2436
[3883]2437 END FUNCTION get_number_of_children
[1791]2438
[2967]2439
[3883]2440 
2441 INTEGER FUNCTION get_childid( id_index )
[2967]2442
[3883]2443    IMPLICIT NONE
[2801]2444
[3883]2445    INTEGER, INTENT(IN) ::  id_index   !<
[2801]2446
[3883]2447   
[2967]2448#if defined( __parallel )
[3883]2449    get_childid = pmc_parent_for_child(id_index)
[2967]2450#else
[3883]2451    get_childid = 0
[2967]2452#endif
[2801]2453
[3883]2454    RETURN
[2967]2455
[3883]2456 END FUNCTION get_childid
[2801]2457
[2967]2458
[3484]2459
[3883]2460 SUBROUTINE get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b, sy_coord, sy_coord_b,   &
2461      ny_coord, ny_coord_b, uz_coord, uz_coord_b )
2462   
2463    IMPLICIT NONE
[2801]2464
[3883]2465    INTEGER,INTENT(IN)   ::  m                     !<
[2801]2466
[3883]2467    REAL(wp),INTENT(OUT) ::  lx_coord, lx_coord_b  !<
2468    REAL(wp),INTENT(OUT) ::  rx_coord, rx_coord_b  !<
[3948]2469    REAL(wp),INTENT(OUT) ::  ny_coord, ny_coord_b  !<
[3883]2470    REAL(wp),INTENT(OUT) ::  sy_coord, sy_coord_b  !<
2471    REAL(wp),INTENT(OUT) ::  uz_coord, uz_coord_b  !<
[2801]2472
[3883]2473   
[3948]2474#if defined( __parallel )
2475   
[3883]2476    lx_coord = childgrid(m)%lx_coord
2477    rx_coord = childgrid(m)%rx_coord
2478    sy_coord = childgrid(m)%sy_coord
2479    ny_coord = childgrid(m)%ny_coord
2480    uz_coord = childgrid(m)%uz_coord
2481   
2482    lx_coord_b = childgrid(m)%lx_coord_b
2483    rx_coord_b = childgrid(m)%rx_coord_b
2484    sy_coord_b = childgrid(m)%sy_coord_b
2485    ny_coord_b = childgrid(m)%ny_coord_b
2486    uz_coord_b = childgrid(m)%uz_coord_b
2487   
[3948]2488#endif
2489   
[3883]2490 END SUBROUTINE get_child_edges
[2801]2491
[3484]2492
2493
[3883]2494 SUBROUTINE  get_child_gridspacing( m, dx, dy,dz )
[2801]2495
[3883]2496    IMPLICIT NONE
2497   
2498    INTEGER, INTENT(IN)             ::  m      !<
[2801]2499
[3883]2500    REAL(wp), INTENT(OUT)           ::  dx,dy  !<
[2801]2501
[3883]2502    REAL(wp), INTENT(OUT), OPTIONAL ::  dz     !<
[2801]2503
[3948]2504
2505#if defined( __parallel )
[3883]2506   
2507    dx = childgrid(m)%dx
2508    dy = childgrid(m)%dy
2509    IF ( PRESENT( dz ) )  THEN
2510       dz = childgrid(m)%dz
2511    ENDIF
2512   
[3948]2513#endif
2514   
[3883]2515 END SUBROUTINE get_child_gridspacing
[2801]2516
[3484]2517
2518
[3883]2519 SUBROUTINE pmci_create_childs_parent_grid_arrays( name, is, ie, js, je, nzc, n  )
2520
[1791]2521    IMPLICIT NONE
2522
[3883]2523    INTEGER(iwp), INTENT(IN) ::  ie      !<  RENAME ie, is, je, js?
[2801]2524    INTEGER(iwp), INTENT(IN) ::  is      !<
2525    INTEGER(iwp), INTENT(IN) ::  je      !<
2526    INTEGER(iwp), INTENT(IN) ::  js      !<
[3883]2527    INTEGER(iwp), INTENT(IN) ::  nzc     !<  nzc is pg%nz, but note that pg%nz is not the original nz of parent,
[3948]2528                                         !<  but the highest parent-grid level needed for nesting.
[3883]2529    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species / salsa variables
2530   
2531    CHARACTER(LEN=*), INTENT(IN) ::  name    !<
[3948]2532
2533#if defined( __parallel )
[3883]2534!       
2535!-- Local variables:
[2801]2536    INTEGER(iwp) ::  ierr    !<
[3883]2537   
2538    INTEGER(idp), POINTER,DIMENSION(:,:)   ::  i_2d    !<
2539   
[2801]2540    REAL(wp), POINTER,DIMENSION(:,:)       ::  p_2d    !<
2541    REAL(wp), POINTER,DIMENSION(:,:,:)     ::  p_3d    !<
[3883]2542   
[1791]2543    NULLIFY( p_3d )
2544    NULLIFY( p_2d )
[2801]2545    NULLIFY( i_2d )
[1791]2546!
[1797]2547!-- List of array names, which can be coupled
[1791]2548    IF ( TRIM( name ) == "u" )  THEN
[2773]2549       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1,js:je,is:ie) )
[1791]2550       p_3d => uc
2551    ELSEIF ( TRIM( name ) == "v" )  THEN
[2773]2552       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1,js:je,is:ie) )
[1791]2553       p_3d => vc
2554    ELSEIF ( TRIM( name ) == "w" )  THEN
[2773]2555       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1,js:je,is:ie) )
[1791]2556       p_3d => wc
2557    ELSEIF ( TRIM( name ) == "e" )  THEN
[2773]2558       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1,js:je,is:ie) )
[1791]2559       p_3d => ec
[2938]2560    ELSEIF ( TRIM( name ) == "diss" )  THEN
2561       IF ( .NOT. ALLOCATED( dissc ) )  ALLOCATE( dissc(0:nzc+1,js:je,is:ie) )
2562       p_3d => dissc
[1791]2563    ELSEIF ( TRIM( name ) == "pt")  THEN
[2773]2564       IF ( .NOT. ALLOCATED( ptc ) )  ALLOCATE( ptc(0:nzc+1,js:je,is:ie) )
[1791]2565       p_3d => ptc
2566    ELSEIF ( TRIM( name ) == "q")  THEN
[2773]2567       IF ( .NOT. ALLOCATED( q_c ) )  ALLOCATE( q_c(0:nzc+1,js:je,is:ie) )
[2174]2568       p_3d => q_c
[2292]2569    ELSEIF ( TRIM( name ) == "qc")  THEN
[2773]2570       IF ( .NOT. ALLOCATED( qcc ) )  ALLOCATE( qcc(0:nzc+1,js:je,is:ie) )
[2292]2571       p_3d => qcc
[2174]2572    ELSEIF ( TRIM( name ) == "qr")  THEN
[2773]2573       IF ( .NOT. ALLOCATED( qrc ) )  ALLOCATE( qrc(0:nzc+1,js:je,is:ie) )
[2174]2574       p_3d => qrc
2575    ELSEIF ( TRIM( name ) == "nr")  THEN
[2773]2576       IF ( .NOT. ALLOCATED( nrc ) )  ALLOCATE( nrc(0:nzc+1,js:je,is:ie) )
[2174]2577       p_3d => nrc
[2292]2578    ELSEIF ( TRIM( name ) == "nc")  THEN
[2773]2579       IF ( .NOT. ALLOCATED( ncc ) )  ALLOCATE( ncc(0:nzc+1,js:je,is:ie) )
[2292]2580       p_3d => ncc
[2003]2581    ELSEIF ( TRIM( name ) == "s")  THEN
[2773]2582       IF ( .NOT. ALLOCATED( sc ) )  ALLOCATE( sc(0:nzc+1,js:je,is:ie) )
[2003]2583       p_3d => sc
[2938]2584    ELSEIF ( TRIM( name ) == "nr_part") THEN
2585       IF ( .NOT. ALLOCATED( nr_partc ) )  ALLOCATE( nr_partc(js:je,is:ie) )
[2801]2586       i_2d => nr_partc
[2938]2587    ELSEIF ( TRIM( name ) == "part_adr") THEN
2588       IF ( .NOT. ALLOCATED( part_adrc ) )  ALLOCATE( part_adrc(js:je,is:ie) )
[2801]2589       i_2d => part_adrc
[2773]2590    ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
[3883]2591       IF ( .NOT. ALLOCATED( chem_spec_c ) ) ALLOCATE( chem_spec_c(0:nzc+1,js:je,is:ie,1:nspec) )
[2773]2592       p_3d => chem_spec_c(:,:,:,n)
[3864]2593    ELSEIF ( TRIM( name(1:3) ) == "an_" )  THEN
2594       IF ( .NOT. ALLOCATED( aerosol_number_c ) )                              &
2595          ALLOCATE( aerosol_number_c(0:nzc+1,js:je,is:ie,1:nbins_aerosol) )
2596       p_3d => aerosol_number_c(:,:,:,n)
2597    ELSEIF ( TRIM( name(1:3) ) == "am_" )  THEN
2598       IF ( .NOT. ALLOCATED( aerosol_mass_c ) )                                &
2599          ALLOCATE( aerosol_mass_c(0:nzc+1,js:je,is:ie,1:(nbins_aerosol*ncomponents_mass) ) )
2600       p_3d => aerosol_mass_c(:,:,:,n)
2601    ELSEIF ( TRIM( name(1:3) ) == "sg_"  .AND.  .NOT. salsa_gases_from_chem )  &
2602    THEN
2603       IF ( .NOT. ALLOCATED( salsa_gas_c ) )                                   &
2604          ALLOCATE( salsa_gas_c(0:nzc+1,js:je,is:ie,1:ngases_salsa) )
2605       p_3d => salsa_gas_c(:,:,:,n)
[1791]2606    !ELSEIF (trim(name) == "z0") then
2607       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
2608       !p_2d => z0c
2609    ENDIF
2610
2611    IF ( ASSOCIATED( p_3d ) )  THEN
2612       CALL pmc_c_set_dataarray( p_3d )
2613    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
2614       CALL pmc_c_set_dataarray( p_2d )
[2801]2615    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
2616       CALL pmc_c_set_dataarray( i_2d )
[1791]2617    ELSE
2618!
[3948]2619!--    Give only one message for the first child domain.
2620       IF ( cpl_id == 2  .AND.  myid == 0 )  THEN
[2801]2621          message_string = 'pointer for array "' // TRIM( name ) //            &
[3883]2622               '" can''t be associated'
2623          CALL message( 'pmci_create_childs_parent_grid_arrays', 'PA0170', 3, 2, 0, 6, 0 )
[1791]2624       ELSE
2625!
[3883]2626!--          Prevent others from continuing in case the abort is to come.
[1791]2627          CALL MPI_BARRIER( comm2d, ierr )
2628       ENDIF
[3883]2629
[1791]2630    ENDIF
2631
[3948]2632#endif
[3883]2633 END SUBROUTINE pmci_create_childs_parent_grid_arrays
[1791]2634
2635
[1933]2636 SUBROUTINE pmci_parent_initialize
2637
2638!
2639!-- Send data for the children in order to let them create initial
2640!-- conditions by interpolating the parent-domain fields.
[1791]2641#if defined( __parallel )
2642    IMPLICIT NONE
2643
[2801]2644    INTEGER(iwp) ::  child_id    !<
2645    INTEGER(iwp) ::  m           !<
2646    REAL(wp) ::  waittime        !<
[1791]2647
2648
[1933]2649    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
2650       child_id = pmc_parent_for_child(m)
2651       CALL pmc_s_fillbuffer( child_id, waittime=waittime )
[1791]2652    ENDDO
2653
2654#endif
[1933]2655 END SUBROUTINE pmci_parent_initialize
[1791]2656
2657
2658
[1933]2659 SUBROUTINE pmci_child_initialize
2660
2661!
2662!-- Create initial conditions for the current child domain by interpolating
2663!-- the parent-domain fields.
[1791]2664#if defined( __parallel )
2665    IMPLICIT NONE
2666
[3883]2667    INTEGER(iwp) ::  ic         !< Child-grid index in x-direction
2668    INTEGER(iwp) ::  jc         !< Child-grid index in y-direction
2669    INTEGER(iwp) ::  kc         !< Child-grid index in z-direction
2670    INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
2671    INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
2672    INTEGER(iwp) ::  lg         !< Running index for salsa gases
2673    INTEGER(iwp) ::  n          !< Running index for chemical species
2674    REAL(wp) ::  waittime       !< Waiting time
[3864]2675
[1791]2676!
[2602]2677!-- Root model is never anyone's child
[3948]2678    IF ( .NOT.  pmc_is_rootmodel() )  THEN
[1791]2679!
[1933]2680!--    Get data from the parent
[1791]2681       CALL pmc_c_getbuffer( waittime = waittime )
2682!
2683!--    The interpolation.
[3792]2684       CALL pmci_interp_1sto_all ( u, uc, kcto, iflu, ifuu, jflo, jfuo, kflo, kfuo, 'u' )
2685       CALL pmci_interp_1sto_all ( v, vc, kcto, iflo, ifuo, jflv, jfuv, kflo, kfuo, 'v' )
2686       CALL pmci_interp_1sto_all ( w, wc, kctw, iflo, ifuo, jflo, jfuo, kflw, kfuw, 'w' )
[2174]2687
[3792]2688       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                              &
2689            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                               &
[2938]2690               .NOT. constant_diffusion ) )  THEN
[3792]2691          CALL pmci_interp_1sto_all ( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 'e' )
[2899]2692       ENDIF
[2938]2693
2694       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
[3792]2695          CALL pmci_interp_1sto_all ( diss, dissc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[2938]2696       ENDIF
2697
[1894]2698       IF ( .NOT. neutral )  THEN
[3792]2699          CALL pmci_interp_1sto_all ( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[1894]2700       ENDIF
[2174]2701
[2003]2702       IF ( humidity )  THEN
[2174]2703
[3792]2704          CALL pmci_interp_1sto_all ( q, q_c, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[2174]2705
[3274]2706          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
[3792]2707             CALL pmci_interp_1sto_all ( qc, qcc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' ) 
2708             CALL pmci_interp_1sto_all ( nc, ncc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )   
[2292]2709          ENDIF
2710
[3274]2711          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
[3792]2712             CALL pmci_interp_1sto_all ( qr, qrc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2713             CALL pmci_interp_1sto_all ( nr, nrc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[2174]2714          ENDIF
2715
[1791]2716       ENDIF
[2174]2717
[2003]2718       IF ( passive_scalar )  THEN
[3792]2719          CALL pmci_interp_1sto_all ( s, sc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[2003]2720       ENDIF
[1791]2721
[4273]2722       IF ( air_chemistry  .AND.  nesting_chem )  THEN
[2773]2723          DO  n = 1, nspec
[3792]2724             CALL pmci_interp_1sto_all ( chem_species(n)%conc, chem_spec_c(:,:,:,n),                &
2725                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
[2773]2726          ENDDO
2727       ENDIF
2728
[4273]2729       IF ( salsa  .AND.  nesting_salsa )  THEN
[3883]2730          DO  lb = 1, nbins_aerosol
2731             CALL pmci_interp_1sto_all ( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),       &
[3864]2732                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2733          ENDDO
[3883]2734          DO  lc = 1, nbins_aerosol * ncomponents_mass
2735             CALL pmci_interp_1sto_all ( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),           &
[3864]2736                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2737          ENDDO
2738          IF ( .NOT. salsa_gases_from_chem )  THEN
[3883]2739             DO  lg = 1, ngases_salsa
2740                CALL pmci_interp_1sto_all ( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),              &
[3864]2741                                            kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
2742             ENDDO
2743          ENDIF
2744       ENDIF
2745
[1791]2746       IF ( topography /= 'flat' )  THEN
2747!
[3883]2748!--       Inside buildings set velocities back to zero.
2749          DO  ic = nxlg, nxrg
2750             DO  jc = nysg, nyng
2751                DO  kc = nzb, nzt
2752                   u(kc,jc,ic)   = MERGE( u(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 1 ) )
2753                   v(kc,jc,ic)   = MERGE( v(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 2 ) )
2754                   w(kc,jc,ic)   = MERGE( w(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 3 ) )
2755                   u_p(kc,jc,ic) = MERGE( u_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 1 ) )
2756                   v_p(kc,jc,ic) = MERGE( v_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 2 ) )
2757                   w_p(kc,jc,ic) = MERGE( w_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 3 ) )
[2232]2758                ENDDO
[1791]2759             ENDDO
2760          ENDDO
2761       ENDIF
2762    ENDIF
2763
2764
2765 CONTAINS
2766
2767
[3883]2768    SUBROUTINE pmci_interp_1sto_all( child_array, parent_array, kct, ifl, ifu, jfl, jfu, kfl, kfu,  &
2769         var )
[1791]2770!
[1933]2771!--    Interpolation of the internal values for the child-domain initialization
[1791]2772       IMPLICIT NONE
2773
[3883]2774       INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
2775
2776       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !<  Indicates start index of child cells belonging to certain
2777                                                               !<  parent cell - x direction
2778       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !<  Indicates end index of child cells belonging to certain
2779                                                               !<  parent cell - x direction
2780       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !<  Indicates start index of child cells belonging to certain
2781                                                               !<  parent cell - y direction
2782       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !<  Indicates end index of child cells belonging to certain
2783                                                               !<  parent cell - y direction
2784       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !<  Indicates start index of child cells belonging to certain
2785                                                               !<  parent cell - z direction
2786       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !<  Indicates end index of child cells belonging to certain
2787                                                               !<  parent cell - z direction
2788       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array  !<  Child-grid array
2789       REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN) ::  parent_array        !<  Parent-grid array
2790
2791       CHARACTER(LEN=1), INTENT(IN) ::  var  !<  Variable symbol: 'u', 'v', 'w' or 's'
[3792]2792!
2793!--    Local variables:
[3883]2794       INTEGER(iwp) ::  ic        !< Running child-grid index in the x-direction
[3976]2795       INTEGER(iwp) ::  icb       !< Index pointing to the first redundant ghost point layer behind the actual boundary
2796                                  !< ghost point layer in the x-direction
2797       INTEGER(iwp) ::  icbc      !< Index pointing to the boundary ghost point layer in the x-direction
[3883]2798       INTEGER(iwp) ::  icfirst   !< Leftmost child-grid index initialized by the main loops (usually icfirst == icl_init)
2799       INTEGER(iwp) ::  iclast    !< Rightmost child-grid index initialized by the main loops (usually iclast == icr_init)
2800       INTEGER(iwp) ::  icl_init  !< Left child-grid index bound for initialization in the x-direction
2801       INTEGER(iwp) ::  icr_init  !< Right child-grid index bound for initialization in the x-direction
2802       INTEGER(iwp) ::  jc        !< Running child-grid index in the y-direction
[3976]2803       INTEGER(iwp) ::  jcb       !< Index pointing to the first redundant ghost point layer behind the actual boundary
2804                                  !< ghost point layer in the y-direction
2805       INTEGER(iwp) ::  jcbc      !< Index pointing to the boundary ghost point layer in the y-direction
[3883]2806       INTEGER(iwp) ::  jcfirst   !< Southmost child-grid index initialized by the main loops (usually jcfirst == jcs_init)
2807       INTEGER(iwp) ::  jclast    !< Northmost child-grid index initialized by the main loops (usually jclast == jcn_init)
2808       INTEGER(iwp) ::  jcs_init  !< South child-grid index bound for initialization in the y-direction
2809       INTEGER(iwp) ::  jcn_init  !< North child-grid index bound for initialization in the y-direction
2810       INTEGER(iwp) ::  kc        !< Running child-grid index in the z-direction
2811       INTEGER(iwp) ::  ip        !< Running parent-grid index in the x-direction
2812       INTEGER(iwp) ::  ipl_init  !< Left parent-grid index bound for initialization in the x-direction
2813       INTEGER(iwp) ::  ipr_init  !< Right parent-grid index bound for initialization in the x-direction
2814       INTEGER(iwp) ::  jp        !< Running parent-grid index in the y-direction
2815       INTEGER(iwp) ::  jps_init  !< South parent-grid index bound for initialization in the y-direction
2816       INTEGER(iwp) ::  jpn_init  !< North parent-grid index bound for initialization in the y-direction
2817       INTEGER(iwp) ::  kp        !< Running parent-grid index in the z-direction
[1791]2818
2819
[3883]2820       ipl_init = ipl
2821       ipr_init = ipr
2822       jps_init = jps
2823       jpn_init = jpn
2824       icl_init = nxl
2825       icr_init = nxr
2826       jcs_init = nys
2827       jcn_init = nyn
[1791]2828
[3976]2829       icbc = -1
2830       icb  = -2
2831       jcbc = -1
2832       jcb  = -2
2833       IF ( var == 'u' )  THEN
2834          icbc =  0
2835          icb  = -1
2836       ELSE IF ( var == 'v' )  THEN
2837          jcbc =  0
2838          jcb  = -1
2839       ENDIF
2840       
[1933]2841       IF ( nesting_mode /= 'vertical' )  THEN
[3182]2842          IF ( bc_dirichlet_l )  THEN
[3883]2843             ipl_init = ipl + 1
2844             icl_init = nxl - 1
[1791]2845!
[1933]2846!--          For u, nxl is a ghost node, but not for the other variables
2847             IF ( var == 'u' )  THEN
[3883]2848                ipl_init = ipl + 2
2849                icl_init = nxl
[1933]2850             ENDIF
[1791]2851          ENDIF
[3182]2852          IF ( bc_dirichlet_s )  THEN
[3883]2853             jps_init = jps + 1
2854             jcs_init = nys - 1 
[1791]2855!
[1933]2856!--          For v, nys is a ghost node, but not for the other variables
2857             IF ( var == 'v' )  THEN
[3883]2858                jps_init = jps + 2
2859                jcs_init = nys
[1933]2860             ENDIF
[1791]2861          ENDIF
[3182]2862          IF ( bc_dirichlet_r )  THEN
[3883]2863             ipr_init = ipr -