source: palm/trunk/SOURCE/nesting_offl_mod.f90 @ 4227

Last change on this file since 4227 was 4227, checked in by gronemeier, 5 years ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

  • Property svn:keywords set to Id
File size: 107.7 KB
RevLine 
[3347]1!> @file nesting_offl_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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.
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
[3347]18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
[4227]22! implement new palm_date_time_mod
[3347]23!
24! Former revisions:
25! -----------------
[3413]26! $Id: nesting_offl_mod.f90 4227 2019-09-10 18:04:34Z gronemeier $
[4226]27! - Data input moved into nesting_offl_mod
28! - check rephrased
29! - time variable is now relative to time_utc_init
30! - Define module specific data type for offline nesting in nesting_offl_mod
31!
32! 4182 2019-08-22 15:20:23Z scharf
[4182]33! Corrected "Former revisions" section
34!
35! 4169 2019-08-19 13:54:35Z suehring
[4169]36! Additional check added.
37!
38! 4168 2019-08-16 13:50:17Z suehring
[4168]39! Replace function get_topography_top_index by topo_top_ind
40!
41! 4125 2019-07-29 13:31:44Z suehring
[4125]42! In order to enable netcdf parallel access, allocate dummy arrays for the
43! lateral boundary data on cores that actually do not belong to these
44! boundaries.
45!
46! 4079 2019-07-09 18:04:41Z suehring
[4079]47! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
48!   though these won't enter the numerical solution. However, due to the mass
49!   conservation these values might some up to very large values which will
50!   occur in the run-control file
51! - Bugfix in offline nesting of chemical species
52! - Do not set Neumann conditions for TKE and passive scalar
53!
54! 4022 2019-06-12 11:52:39Z suehring
[4022]55! Detection of boundary-layer depth in stable boundary layer on basis of
56! boundary data improved
57! Routine for boundary-layer depth calculation renamed and made public
58!
59! 3987 2019-05-22 09:52:13Z kanani
[3987]60! Introduce alternative switch for debug output during timestepping
61!
62! 3964 2019-05-09 09:48:32Z suehring
[3964]63! Ensure that veloctiy term in calculation of bulk Richardson number does not
64! become zero
65!
66! 3937 2019-04-29 15:09:07Z suehring
[3937]67! Set boundary conditon on upper-left and upper-south grid point for the u- and
68! v-component, respectively.
69!
70! 3891 2019-04-12 17:52:01Z suehring
[3891]71! Bugfix, do not overwrite lateral and top boundary data in case of restart
72! runs.
73!
74! 3885 2019-04-11 11:29:34Z kanani
[3885]75! Changes related to global restructuring of location messages and introduction
76! of additional debug messages
77!
78!
[3858]79! Do local data exchange for chemistry variables only when boundary data is 
80! coming from dynamic file
81!
82! 3737 2019-02-12 16:57:06Z suehring
[3737]83! Introduce mesoscale nesting for chemical species
84!
85! 3705 2019-01-29 19:56:39Z suehring
[3705]86! Formatting adjustments
87!
88! 3704 2019-01-29 19:51:41Z suehring
[3579]89! Check implemented for offline nesting in child domain
90!
[4182]91! Initial Revision:
92! - separate offline nesting from large_scale_nudging_mod
93! - revise offline nesting, adjust for usage of synthetic turbulence generator
94! - adjust Rayleigh damping depending on the time-depending atmospheric
95!   conditions
[3413]96!
[4182]97!
[3347]98! Description:
99! ------------
100!> Offline nesting in larger-scale models. Boundary conditions for the simulation
101!> are read from NetCDF file and are prescribed onto the respective arrays.
102!> Further, a mass-flux correction is performed to maintain the mass balance.
103!--------------------------------------------------------------------------------!
104 MODULE nesting_offl_mod
105
106    USE arrays_3d,                                                             &
[4226]107        ONLY:  dzw,                                                            &
108               e,                                                              &
109               diss,                                                           &
110               pt,                                                             &
111               pt_init,                                                        &
112               q,                                                              &
113               q_init,                                                         &
114               rdf,                                                            &
115               rdf_sc,                                                         &
116               s,                                                              &
117               u,                                                              &
118               u_init,                                                         &
119               ug,                                                             &
120               v,                                                              &
121               v_init,                                                         &
122               vg,                                                             &
123               w,                                                              &
124               zu,                                                             &
125               zw
126
127    USE basic_constants_and_equations_mod,                                     &
128           ONLY:  g,                                                           &
129                  pi
[3737]130                 
[3876]131    USE chem_modules,                                                          &
[3737]132        ONLY:  chem_species
[3347]133
134    USE control_parameters,                                                    &
[4169]135        ONLY:  air_chemistry,                                                  & 
136               bc_dirichlet_l,                                                 & 
137               bc_dirichlet_n,                                                 &
138               bc_dirichlet_r,                                                 &
139               bc_dirichlet_s,                                                 &
[4226]140               coupling_char,                                                  &
[4169]141               dt_3d,                                                          &
142               dz,                                                             &
143               constant_diffusion,                                             &
144               child_domain,                                                   &
[3987]145               debug_output_timestep,                                          &
[4169]146               end_time,                                                       &
[3987]147               humidity,                                                       &
148               initializing_actions,                                           &
[4169]149               message_string,                                                 &
150               nesting_offline,                                                &
151               neutral,                                                        &
152               passive_scalar,                                                 &
153               rans_mode,                                                      &
154               rans_tke_e,                                                     &
155               rayleigh_damping_factor,                                        & 
156               rayleigh_damping_height,                                        &
157               spinup_time,                                                    &
158               time_since_reference_point,                                     &
159               volume_flow
[3347]160               
161    USE cpulog,                                                                &
[4226]162        ONLY:  cpu_log,                                                        &
163               log_point,                                                      &
164               log_point_s
[3347]165
166    USE grid_variables
167
168    USE indices,                                                               &
169        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
[4168]170               nysv, nysg, nyn, nyng, nzb, nz, nzt,                            &
171               topo_top_ind,                                                   &
172               wall_flags_0
[3347]173
174    USE kinds
175
176    USE netcdf_data_input_mod,                                                 &
[4226]177        ONLY:  check_existence,                                                &
178               close_input_file,                                               &
179               get_dimension_length,                                           &
180               get_variable,                                                   &
181               get_variable_pr,                                                &
182               input_pids_dynamic,                                             &
183               inquire_num_variables,                                          &
184               inquire_variable_names,                                         &
185               input_file_dynamic,                                             &
186               num_var_pids,                                                   &
187               open_read_file,                                                 &
188               pids_id
189
[4227]190    USE palm_date_time_mod,                                                    &
191        ONLY:  get_date_time
192
[3347]193    USE pegrid
194
[4226]195    IMPLICIT NONE
196
197!
198!-- Define data type for nesting in larger-scale models like COSMO.
199!-- Data type comprises u, v, w, pt, and q at lateral and top boundaries.
200    TYPE nest_offl_type
201
202       CHARACTER(LEN=16) ::  char_l = 'ls_forcing_left_'  !< leading substring for variables at left boundary
203       CHARACTER(LEN=17) ::  char_n = 'ls_forcing_north_' !< leading substring for variables at north boundary 
204       CHARACTER(LEN=17) ::  char_r = 'ls_forcing_right_' !< leading substring for variables at right boundary 
205       CHARACTER(LEN=17) ::  char_s = 'ls_forcing_south_' !< leading substring for variables at south boundary
206       CHARACTER(LEN=15) ::  char_t = 'ls_forcing_top_'   !< leading substring for variables at top boundary
207
208       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names         !< list of variable in dynamic input file
209       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_l  !< names of mesoscale nested chemistry variables at left boundary
210       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_n  !< names of mesoscale nested chemistry variables at north boundary
211       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_r  !< names of mesoscale nested chemistry variables at right boundary
212       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_s  !< names of mesoscale nested chemistry variables at south boundary
213       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_t  !< names of mesoscale nested chemistry variables at top boundary
214
215       INTEGER(iwp) ::  nt     !< number of time levels in dynamic input file
216       INTEGER(iwp) ::  nzu    !< number of vertical levels on scalar grid in dynamic input file
217       INTEGER(iwp) ::  nzw    !< number of vertical levels on w grid in dynamic input file
218       INTEGER(iwp) ::  tind   !< time index for reference time in mesoscale-offline nesting
219       INTEGER(iwp) ::  tind_p !< time index for following time in mesoscale-offline nesting
220
221       LOGICAL      ::  init         = .FALSE. !< flag indicating that offline nesting is already initialized
222
223       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_l !< flags inidicating whether left boundary data for chemistry is in dynamic input file 
224       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_n !< flags inidicating whether north boundary data for chemistry is in dynamic input file
225       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_r !< flags inidicating whether right boundary data for chemistry is in dynamic input file
226       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_s !< flags inidicating whether south boundary data for chemistry is in dynamic input file
227       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_t !< flags inidicating whether top boundary data for chemistry is in dynamic input file
228
229       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surface_pressure !< time dependent surface pressure
230       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time             !< time levels in dynamic input file
231       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos         !< vertical levels at scalar grid in dynamic input file
232       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos         !< vertical levels at w grid in dynamic input file
233
234       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug         !< domain-averaged geostrophic component
235       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg         !< domain-averaged geostrophic component
236
237       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_left   !< u-component at left boundary
238       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_left   !< v-component at left boundary
239       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_left   !< w-component at left boundary
240       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_left   !< mixing ratio at left boundary
241       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_left  !< potentital temperautre at left boundary
242
243       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_north  !< u-component at north boundary
244       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_north  !< v-component at north boundary
245       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_north  !< w-component at north boundary
246       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_north  !< mixing ratio at north boundary
247       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_north !< potentital temperautre at north boundary
248
249       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_right  !< u-component at right boundary
250       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_right  !< v-component at right boundary
251       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_right  !< w-component at right boundary
252       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_right  !< mixing ratio at right boundary
253       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_right !< potentital temperautre at right boundary
254
255       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_south  !< u-component at south boundary
256       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_south  !< v-component at south boundary
257       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_south  !< w-component at south boundary
258       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_south  !< mixing ratio at south boundary
259       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_south !< potentital temperautre at south boundary
260
261       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_top    !< u-component at top boundary
262       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_top    !< v-component at top boundary
263       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_top    !< w-component at top boundary
264       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_top    !< mixing ratio at top boundary
265       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_top   !< potentital temperautre at top boundary
266       
267       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_left   !< chemical species at left boundary
268       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_north  !< chemical species at left boundary
269       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_right  !< chemical species at left boundary
270       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_south  !< chemical species at left boundary
271       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_top    !< chemical species at left boundary
272
273    END TYPE nest_offl_type
274
275    REAL(wp) ::  fac_dt              !< interpolation factor
[4227]276    REAL(wp) ::  time_utc_init       !< time in seconds-of-day of origin_date_time
[3347]277    REAL(wp) ::  zi_ribulk = 0.0_wp  !< boundary-layer depth according to bulk Richardson criterion, i.e. the height where Ri_bulk exceeds the critical
[4226]278                                     !< bulk Richardson number of 0.2
279
280    TYPE(nest_offl_type) ::  nest_offl  !< data structure for data input at lateral and top boundaries (provided by Inifor)
[3347]281   
282    SAVE
283    PRIVATE
284!
285!-- Public subroutines
[4022]286    PUBLIC nesting_offl_bc,                                                    &
287           nesting_offl_calc_zi,                                               &
288           nesting_offl_check_parameters,                                      &
[4226]289           nesting_offl_geostrophic_wind,                                      &
[4022]290           nesting_offl_header,                                                &
291           nesting_offl_init,                                                  &
[4226]292           nesting_offl_input,                                                 &
293           nesting_offl_interpolation_factor,                                  &
[4022]294           nesting_offl_mass_conservation,                                     &
295           nesting_offl_parin 
[3347]296!
297!-- Public variables
298    PUBLIC zi_ribulk   
299
300    INTERFACE nesting_offl_bc
301       MODULE PROCEDURE nesting_offl_bc
302    END INTERFACE nesting_offl_bc
303   
[4022]304    INTERFACE nesting_offl_calc_zi
305       MODULE PROCEDURE nesting_offl_calc_zi
306    END INTERFACE nesting_offl_calc_zi
307   
[3579]308    INTERFACE nesting_offl_check_parameters
309       MODULE PROCEDURE nesting_offl_check_parameters
310    END INTERFACE nesting_offl_check_parameters
[4226]311
312    INTERFACE nesting_offl_geostrophic_wind
313       MODULE PROCEDURE nesting_offl_geostrophic_wind
314    END INTERFACE nesting_offl_geostrophic_wind
[3579]315   
[3347]316    INTERFACE nesting_offl_header
317       MODULE PROCEDURE nesting_offl_header
318    END INTERFACE nesting_offl_header
319   
320    INTERFACE nesting_offl_init
321       MODULE PROCEDURE nesting_offl_init
322    END INTERFACE nesting_offl_init
[4226]323
324    INTERFACE nesting_offl_input
325       MODULE PROCEDURE nesting_offl_input
326    END INTERFACE nesting_offl_input
327
328    INTERFACE nesting_offl_interpolation_factor
329       MODULE PROCEDURE nesting_offl_interpolation_factor
330    END INTERFACE nesting_offl_interpolation_factor
[3347]331           
332    INTERFACE nesting_offl_mass_conservation
333       MODULE PROCEDURE nesting_offl_mass_conservation
334    END INTERFACE nesting_offl_mass_conservation
335   
336    INTERFACE nesting_offl_parin
337       MODULE PROCEDURE nesting_offl_parin
338    END INTERFACE nesting_offl_parin
339
340 CONTAINS
341
[4226]342!------------------------------------------------------------------------------!
343! Description:
344! ------------
345!> Reads data at lateral and top boundaries derived from larger-scale model.
346!------------------------------------------------------------------------------!
347    SUBROUTINE nesting_offl_input
[3347]348
[4226]349       INTEGER(iwp) ::  n   !< running index for chemistry variables
350       INTEGER(iwp) ::  t   !< running index time dimension
351
352!
353!--    Initialize INIFOR forcing in first call.
354       IF ( .NOT. nest_offl%init )  THEN
355#if defined ( __netcdf )
356!
357!--       Open file in read-only mode
358          CALL open_read_file( TRIM( input_file_dynamic ) //                   &
359                               TRIM( coupling_char ), pids_id )
360!
361!--       At first, inquire all variable names.
362          CALL inquire_num_variables( pids_id, num_var_pids )
363!
364!--       Allocate memory to store variable names.
365          ALLOCATE( nest_offl%var_names(1:num_var_pids) )
366          CALL inquire_variable_names( pids_id, nest_offl%var_names )
367!
368!--       Read time dimension, allocate memory and finally read time array
369          CALL get_dimension_length( pids_id, nest_offl%nt, 'time' )
370
371          IF ( check_existence( nest_offl%var_names, 'time' ) )  THEN
372             ALLOCATE( nest_offl%time(0:nest_offl%nt-1) )
373             CALL get_variable( pids_id, 'time', nest_offl%time )
374          ENDIF
375!
376!--       Read vertical dimension of scalar und w grid
377          CALL get_dimension_length( pids_id, nest_offl%nzu, 'z' )
378          CALL get_dimension_length( pids_id, nest_offl%nzw, 'zw' )
379
380          IF ( check_existence( nest_offl%var_names, 'z' ) )  THEN
381             ALLOCATE( nest_offl%zu_atmos(1:nest_offl%nzu) )
382             CALL get_variable( pids_id, 'z', nest_offl%zu_atmos )
383          ENDIF
384          IF ( check_existence( nest_offl%var_names, 'zw' ) )  THEN
385             ALLOCATE( nest_offl%zw_atmos(1:nest_offl%nzw) )
386             CALL get_variable( pids_id, 'zw', nest_offl%zw_atmos )
387          ENDIF
388!
389!--       Read surface pressure
390          IF ( check_existence( nest_offl%var_names,                           &
391                                'surface_forcing_surface_pressure' ) )  THEN
392             ALLOCATE( nest_offl%surface_pressure(0:nest_offl%nt-1) )
393             CALL get_variable( pids_id,                                       &
394                                'surface_forcing_surface_pressure',            &
395                                nest_offl%surface_pressure )
396          ENDIF
397!
398!--       Close input file
399          CALL close_input_file( pids_id )
400#endif
401       ENDIF
402!
403!--    Check if dynamic driver data input is required.
404       IF ( nest_offl%time(nest_offl%tind_p) <=                                &
405            MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.                   &
406            .NOT.  nest_offl%init )  THEN
407          CONTINUE
408!
409!--    Return otherwise
410       ELSE
411          RETURN
412       ENDIF
413!
414!--    CPU measurement
415       CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'start' )
416
417!
418!--    Obtain time index for current point in time. Note, the time coordinate
419!--    in the input file is relative to time_utc_init. Since time_since_...
420!--    is negativ when spinup is used, use MAX function to obtain correct
421!--    time at the beginning.
422       nest_offl%tind = MINLOC( ABS( nest_offl%time - (                        &
423                                     time_utc_init +                           &
424                                     MAX( time_since_reference_point, 0.0_wp) )&
425                                   ), DIM = 1 ) - 1
426       nest_offl%tind_p = nest_offl%tind + 1
427!
428!--    Open file in read-only mode
429#if defined ( __netcdf )
430       CALL open_read_file( TRIM( input_file_dynamic ) //                      &
431                            TRIM( coupling_char ), pids_id )
432!
433!--    Read geostrophic wind components
434       DO  t = nest_offl%tind, nest_offl%tind_p
435          CALL get_variable_pr( pids_id, 'ls_forcing_ug', t+1,                 &
436                                nest_offl%ug(t-nest_offl%tind,nzb+1:nzt) )
437          CALL get_variable_pr( pids_id, 'ls_forcing_vg', t+1,                 &
438                                nest_offl%vg(t-nest_offl%tind,nzb+1:nzt) )
439       ENDDO
440!
441!--    Read data at lateral and top boundaries. Please note, at left and
442!--    right domain boundary, yz-layers are read for u, v, w, pt and q.
443!--    For the v-component, the data starts at nysv, while for the other
444!--    quantities the data starts at nys. This is equivalent at the north
445!--    and south domain boundary for the u-component.
446!--    Note, lateral data is also accessed by parallel IO, which is the reason
447!--    why different arguments are passed depending on the boundary control
448!--    flags. Cores that do not belong to the respective boundary just make
449!--    a dummy read with count = 0, just in order to participate the collective
450!--    operation.
451!--    Read data for western boundary   
452       CALL get_variable( pids_id, 'ls_forcing_left_u',                        &
453                          nest_offl%u_left,                                    & ! array to be read
454                          MERGE( nys+1, 1, bc_dirichlet_l),                    & ! start index y direction
455                          MERGE( nzb+1, 1, bc_dirichlet_l),                    & ! start index z direction
456                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         & ! start index time dimension
457                          MERGE( nyn-nys+1, 0, bc_dirichlet_l),                & ! number of elements along y
458                          MERGE( nest_offl%nzu, 0, bc_dirichlet_l),            & ! number of elements alogn z
459                          MERGE( 2, 0, bc_dirichlet_l),                        & ! number of time steps (2 or 0)
460                          .TRUE. )                                               ! parallel IO when compiled accordingly
461     
462       CALL get_variable( pids_id, 'ls_forcing_left_v',                        &
463                          nest_offl%v_left,                                    &
464                          MERGE( nysv, 1, bc_dirichlet_l),                     &
465                          MERGE( nzb+1, 1, bc_dirichlet_l),                    &
466                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         &
467                          MERGE( nyn-nysv+1, 0, bc_dirichlet_l),               &
468                          MERGE( nest_offl%nzu, 0, bc_dirichlet_l),            &
469                          MERGE( 2, 0, bc_dirichlet_l),                        &
470                          .TRUE. )                                       
471
472       CALL get_variable( pids_id, 'ls_forcing_left_w',                        &
473                          nest_offl%w_left,                                    &
474                          MERGE( nys+1, 1, bc_dirichlet_l),                    &
475                          MERGE( nzb+1, 1, bc_dirichlet_l),                    &
476                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         &
477                          MERGE( nyn-nys+1, 0, bc_dirichlet_l),                &
478                          MERGE( nest_offl%nzw, 0, bc_dirichlet_l),            &
479                          MERGE( 2, 0, bc_dirichlet_l),                        &
480                          .TRUE. )   
481
482       IF ( .NOT. neutral )  THEN
483          CALL get_variable( pids_id, 'ls_forcing_left_pt',                    &
484                             nest_offl%pt_left,                                &
485                             MERGE( nys+1, 1, bc_dirichlet_l),                 &
486                             MERGE( nzb+1, 1, bc_dirichlet_l),                 &
487                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),      &
488                             MERGE( nyn-nys+1, 0, bc_dirichlet_l),             &
489                             MERGE( nest_offl%nzu, 0, bc_dirichlet_l),         &
490                             MERGE( 2, 0, bc_dirichlet_l),                     &
491                             .TRUE. )
492       ENDIF
493
494       IF ( humidity )  THEN
495          CALL get_variable( pids_id, 'ls_forcing_left_qv',                    &
496                             nest_offl%q_left,                                 &
497                             MERGE( nys+1, 1, bc_dirichlet_l),                 &
498                             MERGE( nzb+1, 1, bc_dirichlet_l),                 &
499                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),      &
500                             MERGE( nyn-nys+1, 0, bc_dirichlet_l),             &
501                             MERGE( nest_offl%nzu, 0, bc_dirichlet_l),         &
502                             MERGE( 2, 0, bc_dirichlet_l),                     &
503                             .TRUE. )
504       ENDIF
505       
506       IF ( air_chemistry )  THEN
507          DO  n = 1, UBOUND(nest_offl%var_names_chem_l, 1)
508             IF ( check_existence( nest_offl%var_names,                        &
509                                   nest_offl%var_names_chem_l(n) ) )  THEN 
510                CALL get_variable( pids_id,                                    &
511                           TRIM( nest_offl%var_names_chem_l(n) ),              &
512                           nest_offl%chem_left(:,:,:,n),                       &
513                           MERGE( nys+1, 1, bc_dirichlet_l),                   &
514                           MERGE( nzb+1, 1, bc_dirichlet_l),                   &
515                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),        &
516                           MERGE( nyn-nys+1, 0, bc_dirichlet_l),               &
517                           MERGE( nest_offl%nzu, 0, bc_dirichlet_l),           &
518                           MERGE( 2, 0, bc_dirichlet_l),                       &
519                           .TRUE. )
520                nest_offl%chem_from_file_l(n) = .TRUE.
521             ENDIF
522          ENDDO
523       ENDIF
524!
525!--    Read data for eastern boundary   
526       CALL get_variable( pids_id, 'ls_forcing_right_u',                       &
527                          nest_offl%u_right,                                   &
528                          MERGE( nys+1, 1, bc_dirichlet_r),                    &
529                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
530                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
531                          MERGE( nyn-nys+1, 0, bc_dirichlet_r),                &
532                          MERGE( nest_offl%nzu, 0, bc_dirichlet_r),            &
533                          MERGE( 2, 0, bc_dirichlet_r),                        &
534                          .TRUE. )                                             
535     
536       CALL get_variable( pids_id, 'ls_forcing_right_v',                       &
537                          nest_offl%v_right,                                   &
538                          MERGE( nysv, 1, bc_dirichlet_r),                     &
539                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
540                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
541                          MERGE( nyn-nysv+1, 0, bc_dirichlet_r),               &
542                          MERGE( nest_offl%nzu, 0, bc_dirichlet_r),            &
543                          MERGE( 2, 0, bc_dirichlet_r),                        &
544                          .TRUE. )                                             
545
546       CALL get_variable( pids_id, 'ls_forcing_right_w',                       &
547                          nest_offl%w_right,                                   &
548                          MERGE( nys+1, 1, bc_dirichlet_r),                    &
549                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
550                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
551                          MERGE( nyn-nys+1, 0, bc_dirichlet_r),                &
552                          MERGE( nest_offl%nzw, 0, bc_dirichlet_r),            &
553                          MERGE( 2, 0, bc_dirichlet_r),                        &
554                          .TRUE. )   
555
556       IF ( .NOT. neutral )  THEN
557          CALL get_variable( pids_id, 'ls_forcing_right_pt',                   &
558                             nest_offl%pt_right,                               &
559                             MERGE( nys+1, 1, bc_dirichlet_r),                 &
560                             MERGE( nzb+1, 1, bc_dirichlet_r),                 &
561                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),      &
562                             MERGE( nyn-nys+1, 0, bc_dirichlet_r),             &
563                             MERGE( nest_offl%nzu, 0, bc_dirichlet_r),         &
564                             MERGE( 2, 0, bc_dirichlet_r),                     &
565                             .TRUE. )
566       ENDIF
567
568       IF ( humidity )  THEN
569          CALL get_variable( pids_id, 'ls_forcing_right_qv',                   &
570                             nest_offl%q_right,                                &
571                             MERGE( nys+1, 1, bc_dirichlet_r),                 &
572                             MERGE( nzb+1, 1, bc_dirichlet_r),                 &
573                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),      &
574                             MERGE( nyn-nys+1, 0, bc_dirichlet_r),             &
575                             MERGE( nest_offl%nzu, 0, bc_dirichlet_r),         &
576                             MERGE( 2, 0, bc_dirichlet_r),                     &
577                             .TRUE. )
578       ENDIF
579       
580       IF ( air_chemistry )  THEN
581          DO  n = 1, UBOUND(nest_offl%var_names_chem_r, 1)
582             IF ( check_existence( nest_offl%var_names,                        &
583                                   nest_offl%var_names_chem_r(n) ) )  THEN     
584                CALL get_variable( pids_id,                                    &
585                           TRIM( nest_offl%var_names_chem_r(n) ),              &
586                           nest_offl%chem_right(:,:,:,n),                      &
587                           MERGE( nys+1, 1, bc_dirichlet_r),                   &
588                           MERGE( nzb+1, 1, bc_dirichlet_r),                   &
589                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),        &
590                           MERGE( nyn-nys+1, 0, bc_dirichlet_r),               &
591                           MERGE( nest_offl%nzu, 0, bc_dirichlet_r),           &
592                           MERGE( 2, 0, bc_dirichlet_r),                       &
593                           .TRUE. )
594                nest_offl%chem_from_file_r(n) = .TRUE.
595             ENDIF
596          ENDDO
597       ENDIF
598!
599!--    Read data for northern boundary
600       CALL get_variable( pids_id, 'ls_forcing_north_u',                       & ! array to be read
601                          nest_offl%u_north,                                   & ! start index x direction
602                          MERGE( nxlu, 1, bc_dirichlet_n ),                    & ! start index z direction
603                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
604                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
605                          MERGE( nxr-nxlu+1, 0, bc_dirichlet_n ),              & ! number of elements alogn z
606                          MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
607                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
608                          .TRUE. )                                             
609                                                                               
610       CALL get_variable( pids_id, 'ls_forcing_north_v',                       & ! array to be read
611                          nest_offl%v_north,                                   & ! start index x direction
612                          MERGE( nxl+1, 1, bc_dirichlet_n ),                   & ! start index z direction
613                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
614                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
615                          MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),               & ! number of elements alogn z
616                          MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
617                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
618                          .TRUE. )                                             
619                                                                               
620       CALL get_variable( pids_id, 'ls_forcing_north_w',                       & ! array to be read
621                          nest_offl%w_north,                                   & ! start index x direction
622                          MERGE( nxl+1, 1, bc_dirichlet_n ),                   & ! start index z direction
623                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
624                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
625                          MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),               & ! number of elements alogn z
626                          MERGE( nest_offl%nzw, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
627                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
628                          .TRUE. )                                             
629                                                                               
630       IF ( .NOT. neutral )  THEN                                             
631          CALL get_variable( pids_id, 'ls_forcing_north_pt',                   & ! array to be read
632                             nest_offl%pt_north,                               & ! start index x direction
633                             MERGE( nxl+1, 1, bc_dirichlet_n ),                & ! start index z direction
634                             MERGE( nzb+1, 1, bc_dirichlet_n ),                & ! start index time dimension
635                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),     & ! number of elements along x
636                             MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),            & ! number of elements alogn z
637                             MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),        & ! number of time steps (2 or 0)
638                             MERGE( 2, 0, bc_dirichlet_n ),                    & ! parallel IO when compiled accordingly
639                             .TRUE. )                                             
640       ENDIF                                                                   
641       IF ( humidity )  THEN                                                   
642          CALL get_variable( pids_id, 'ls_forcing_north_qv',                   & ! array to be read
643                             nest_offl%q_north,                                & ! start index x direction
644                             MERGE( nxl+1, 1, bc_dirichlet_n ),                & ! start index z direction
645                             MERGE( nzb+1, 1, bc_dirichlet_n ),                & ! start index time dimension
646                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),     & ! number of elements along x
647                             MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),            & ! number of elements alogn z
648                             MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),        & ! number of time steps (2 or 0)
649                             MERGE( 2, 0, bc_dirichlet_n ),                    & ! parallel IO when compiled accordingly
650                             .TRUE. )                                             
651       ENDIF                                                                   
652                                                                               
653       IF ( air_chemistry )  THEN                                             
654          DO  n = 1, UBOUND(nest_offl%var_names_chem_n, 1)                     
655             IF ( check_existence( nest_offl%var_names,                        &
656                                   nest_offl%var_names_chem_n(n) ) )  THEN     
657                CALL get_variable( pids_id,                                    &
658                           TRIM( nest_offl%var_names_chem_n(n) ),              &
659                           nest_offl%chem_north(:,:,:,n),                      &
660                           MERGE( nxl+1, 1, bc_dirichlet_n ),                  &
661                           MERGE( nzb+1, 1, bc_dirichlet_n ),                  &
662                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),       &
663                           MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),              &
664                           MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),          &
665                           MERGE( 2, 0, bc_dirichlet_n ),                      &
666                           .TRUE. )
667                nest_offl%chem_from_file_n(n) = .TRUE.
668             ENDIF
669          ENDDO
670       ENDIF
671!
672!--    Read data for southern boundary
673       CALL get_variable( pids_id, 'ls_forcing_south_u',                       & ! array to be read
674                          nest_offl%u_south,                                   & ! start index x direction
675                          MERGE( nxlu, 1, bc_dirichlet_s ),                    & ! start index z direction
676                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
677                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
678                          MERGE( nxr-nxlu+1, 0, bc_dirichlet_s ),              & ! number of elements alogn z
679                          MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
680                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
681                          .TRUE. )                                             
682                                                                               
683       CALL get_variable( pids_id, 'ls_forcing_south_v',                       & ! array to be read
684                          nest_offl%v_south,                                   & ! start index x direction
685                          MERGE( nxl+1, 1, bc_dirichlet_s ),                   & ! start index z direction
686                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
687                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
688                          MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),               & ! number of elements alogn z
689                          MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
690                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
691                          .TRUE. )                                             
692                                                                               
693       CALL get_variable( pids_id, 'ls_forcing_south_w',                       & ! array to be read
694                          nest_offl%w_south,                                   & ! start index x direction
695                          MERGE( nxl+1, 1, bc_dirichlet_s ),                   & ! start index z direction
696                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
697                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
698                          MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),               & ! number of elements alogn z
699                          MERGE( nest_offl%nzw, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
700                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
701                          .TRUE. )                                             
702                                                                               
703       IF ( .NOT. neutral )  THEN                                             
704          CALL get_variable( pids_id, 'ls_forcing_south_pt',                   & ! array to be read
705                             nest_offl%pt_south,                               & ! start index x direction
706                             MERGE( nxl+1, 1, bc_dirichlet_s ),                & ! start index z direction
707                             MERGE( nzb+1, 1, bc_dirichlet_s ),                & ! start index time dimension
708                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),     & ! number of elements along x
709                             MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),            & ! number of elements alogn z
710                             MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),        & ! number of time steps (2 or 0)
711                             MERGE( 2, 0, bc_dirichlet_s ),                    & ! parallel IO when compiled accordingly
712                             .TRUE. )                                             
713       ENDIF                                                                   
714       IF ( humidity )  THEN                                                   
715          CALL get_variable( pids_id, 'ls_forcing_south_qv',                   & ! array to be read
716                             nest_offl%q_south,                                & ! start index x direction
717                             MERGE( nxl+1, 1, bc_dirichlet_s ),                & ! start index z direction
718                             MERGE( nzb+1, 1, bc_dirichlet_s ),                & ! start index time dimension
719                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),     & ! number of elements along x
720                             MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),            & ! number of elements alogn z
721                             MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),        & ! number of time steps (2 or 0)
722                             MERGE( 2, 0, bc_dirichlet_s ),                    & ! parallel IO when compiled accordingly
723                             .TRUE. )                                             
724       ENDIF                                                                   
725                                                                               
726       IF ( air_chemistry )  THEN                                             
727          DO  n = 1, UBOUND(nest_offl%var_names_chem_s, 1)                     
728             IF ( check_existence( nest_offl%var_names,                        &
729                                   nest_offl%var_names_chem_s(n) ) )  THEN     
730                CALL get_variable( pids_id,                                    &
731                           TRIM( nest_offl%var_names_chem_s(n) ),              &
732                           nest_offl%chem_south(:,:,:,n),                      &
733                           MERGE( nxl+1, 1, bc_dirichlet_s ),                  &
734                           MERGE( nzb+1, 1, bc_dirichlet_s ),                  &
735                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),       &
736                           MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),              &
737                           MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),          &
738                           MERGE( 2, 0, bc_dirichlet_s ),                      &
739                           .TRUE. )
740                nest_offl%chem_from_file_s(n) = .TRUE.
741             ENDIF
742          ENDDO
743       ENDIF
744!
745!--    Top boundary
746       CALL get_variable( pids_id, 'ls_forcing_top_u',                         &
747                             nest_offl%u_top(0:1,nys:nyn,nxlu:nxr),            &
748                             nxlu, nys+1, nest_offl%tind+1,                    &
749                             nxr-nxlu+1, nyn-nys+1, 2, .TRUE. )
750
751       CALL get_variable( pids_id, 'ls_forcing_top_v',                         &
752                             nest_offl%v_top(0:1,nysv:nyn,nxl:nxr),            &
753                             nxl+1, nysv, nest_offl%tind+1,                    &
754                             nxr-nxl+1, nyn-nysv+1, 2, .TRUE. )
755                             
756       CALL get_variable( pids_id, 'ls_forcing_top_w',                         &
757                             nest_offl%w_top(0:1,nys:nyn,nxl:nxr),             &
758                             nxl+1, nys+1, nest_offl%tind+1,                   &
759                             nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
760                             
761       IF ( .NOT. neutral )  THEN
762          CALL get_variable( pids_id, 'ls_forcing_top_pt',                     &
763                                nest_offl%pt_top(0:1,nys:nyn,nxl:nxr),         &
764                                nxl+1, nys+1, nest_offl%tind+1,                &
765                                nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
766       ENDIF
767       IF ( humidity )  THEN
768          CALL get_variable( pids_id, 'ls_forcing_top_qv',                     &
769                                nest_offl%q_top(0:1,nys:nyn,nxl:nxr),          &
770                                nxl+1, nys+1, nest_offl%tind+1,                &
771                                nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
772       ENDIF
773       
774       IF ( air_chemistry )  THEN
775          DO  n = 1, UBOUND(nest_offl%var_names_chem_t, 1)
776             IF ( check_existence( nest_offl%var_names,                        &
777                                   nest_offl%var_names_chem_t(n) ) )  THEN     
778                CALL get_variable( pids_id,                                    &
779                              TRIM( nest_offl%var_names_chem_t(n) ),           &
780                              nest_offl%chem_top(0:1,nys:nyn,nxl:nxr,n),       &
781                              nxl+1, nys+1, nest_offl%tind+1,                  &
782                              nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
783                nest_offl%chem_from_file_t(n) = .TRUE.
784             ENDIF
785          ENDDO
786       ENDIF
787
788!
789!--    Close input file
790       CALL close_input_file( pids_id )
791#endif
792!
793!--    Set control flag to indicate that boundary data has been initially
794!--    input.
795       nest_offl%init = .TRUE.
796!
797!--    End of CPU measurement
798       CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
799
800    END SUBROUTINE nesting_offl_input
801
802
[3347]803!------------------------------------------------------------------------------!
804! Description:
805! ------------
806!> In this subroutine a constant mass within the model domain is guaranteed.
807!> Larger-scale models may be based on a compressible equation system, which is
808!> not consistent with PALMs incompressible equation system. In order to avoid
809!> a decrease or increase of mass during the simulation, non-divergent flow
810!> through the lateral and top boundaries is compensated by the vertical wind
811!> component at the top boundary.
812!------------------------------------------------------------------------------!
813    SUBROUTINE nesting_offl_mass_conservation
814
815       INTEGER(iwp) ::  i !< grid index in x-direction
816       INTEGER(iwp) ::  j !< grid index in y-direction
817       INTEGER(iwp) ::  k !< grid index in z-direction
818
819       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
820       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
821       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
822
[3987]823
824       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
825
[3347]826       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
827       
828       volume_flow   = 0.0_wp
829       volume_flow_l = 0.0_wp
830
831       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
832
833       IF ( bc_dirichlet_l )  THEN
834          i = nxl
835          DO  j = nys, nyn
836             DO  k = nzb+1, nzt
837                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
838                                     * MERGE( 1.0_wp, 0.0_wp,                  &
839                                              BTEST( wall_flags_0(k,j,i), 1 ) )
840             ENDDO
841          ENDDO
842       ENDIF
843       IF ( bc_dirichlet_r )  THEN
844          i = nxr+1
845          DO  j = nys, nyn
846             DO  k = nzb+1, nzt
847                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
848                                     * MERGE( 1.0_wp, 0.0_wp,                  &
849                                              BTEST( wall_flags_0(k,j,i), 1 ) )
850             ENDDO
851          ENDDO
852       ENDIF
853       IF ( bc_dirichlet_s )  THEN
854          j = nys
855          DO  i = nxl, nxr
856             DO  k = nzb+1, nzt
857                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
858                                     * MERGE( 1.0_wp, 0.0_wp,                  &
859                                              BTEST( wall_flags_0(k,j,i), 2 ) )
860             ENDDO
861          ENDDO
862       ENDIF
863       IF ( bc_dirichlet_n )  THEN
864          j = nyn+1
865          DO  i = nxl, nxr
866             DO  k = nzb+1, nzt
867                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
868                                     * MERGE( 1.0_wp, 0.0_wp,                  &
869                                              BTEST( wall_flags_0(k,j,i), 2 ) )
870             ENDDO
871          ENDDO
872       ENDIF
873!
874!--    Top boundary
875       k = nzt
876       DO  i = nxl, nxr
877          DO  j = nys, nyn
878             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
879          ENDDO
880       ENDDO
881
882#if defined( __parallel )
883       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[4226]884       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,   &
[3347]885                           comm2d, ierr )
886#else
887       volume_flow = volume_flow_l
888#endif
889
890       w_correct = SUM( volume_flow ) * d_area_t
891
892       DO  i = nxl, nxr
893          DO  j = nys, nyn
894             DO  k = nzt, nzt + 1
895                w(k,j,i) = w(k,j,i) + w_correct
896             ENDDO
897          ENDDO
898       ENDDO
899       
900       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
901
[3987]902       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
903
[3347]904    END SUBROUTINE nesting_offl_mass_conservation
905
906
907!------------------------------------------------------------------------------!
908! Description:
909! ------------
910!> Set the lateral and top boundary conditions in case the PALM domain is
911!> nested offline in a mesoscale model. Further, average boundary data and
912!> determine mean profiles, further used for correct damping in the sponge
913!> layer.
914!------------------------------------------------------------------------------!
915    SUBROUTINE nesting_offl_bc                     
916
917       INTEGER(iwp) ::  i !< running index x-direction
918       INTEGER(iwp) ::  j !< running index y-direction
919       INTEGER(iwp) ::  k !< running index z-direction
[3737]920       INTEGER(iwp) ::  n !< running index for chemical species
[3347]921       
922       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
923       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
924       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
925       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
926       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
927       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
928       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
929       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
930       
[3885]931
[3987]932       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
933
[3347]934       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
935!
936!--    Set mean profiles, derived from boundary data, to zero
937       pt_ref   = 0.0_wp
938       q_ref    = 0.0_wp
939       u_ref    = 0.0_wp
940       v_ref    = 0.0_wp
941       
942       pt_ref_l = 0.0_wp
943       q_ref_l  = 0.0_wp
944       u_ref_l  = 0.0_wp
945       v_ref_l  = 0.0_wp
946!
947!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
948!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
949!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
950!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
951!--    Please note, at the left (for u) and south (for v) boundary, values
952!--    for u and v are set also at i/j=-1, since these values are used in
953!--    boundary_conditions() to restore prognostic values.
954!--    Further, sum up data to calculate mean profiles from boundary data,
955!--    used for Rayleigh damping.
956       IF ( bc_dirichlet_l )  THEN
957
958          DO  j = nys, nyn
959             DO  k = nzb+1, nzt
960                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
961                                                nest_offl%u_left(1,k,j),       &
962                                                fac_dt ) *                     &
963                             MERGE( 1.0_wp, 0.0_wp,                            &
964                                    BTEST( wall_flags_0(k,j,0), 1 ) )
965                u(k,j,-1) = u(k,j,0)
966             ENDDO
967             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
968          ENDDO
969
970          DO  j = nys, nyn
971             DO  k = nzb+1, nzt-1
972                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
973                                                 nest_offl%w_left(1,k,j),      &
974                                                 fac_dt ) *                    &
975                            MERGE( 1.0_wp, 0.0_wp,                             &
976                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
977             ENDDO
[4079]978             w(nzt,j,-1) = w(nzt-1,j,-1)
[3347]979          ENDDO
980
981          DO  j = nysv, nyn
982             DO  k = nzb+1, nzt
983                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
984                                                 nest_offl%v_left(1,k,j),      &
985                                                 fac_dt ) *                    &
986                               MERGE( 1.0_wp, 0.0_wp,                          &
987                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
988             ENDDO
989             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
990          ENDDO
991
992          IF ( .NOT. neutral )  THEN
993             DO  j = nys, nyn
994                DO  k = nzb+1, nzt
995                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
996                                                     nest_offl%pt_left(1,k,j), &
997                                                     fac_dt )
998 
999                ENDDO
1000                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
1001             ENDDO
1002          ENDIF
1003
1004          IF ( humidity )  THEN
1005             DO  j = nys, nyn
1006                DO  k = nzb+1, nzt
1007                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
1008                                                    nest_offl%q_left(1,k,j),   &
1009                                                    fac_dt )
1010 
1011                ENDDO
1012                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
1013             ENDDO
1014          ENDIF
[3737]1015         
1016          IF ( air_chemistry )  THEN
1017             DO  n = 1, UBOUND( chem_species, 1 )
1018                IF ( nest_offl%chem_from_file_l(n) )  THEN                   
1019                   DO  j = nys, nyn
1020                      DO  k = nzb+1, nzt
1021                         chem_species(n)%conc(k,j,-1) = interpolate_in_time(   &
1022                                                  nest_offl%chem_left(0,k,j,n),&
1023                                                  nest_offl%chem_left(1,k,j,n),&
1024                                                  fac_dt                   )
1025                      ENDDO
1026                   ENDDO
1027                ENDIF
1028             ENDDO
1029          ENDIF
[3347]1030
1031       ENDIF
1032
1033       IF ( bc_dirichlet_r )  THEN
1034
1035          DO  j = nys, nyn
1036             DO  k = nzb+1, nzt
1037                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
1038                                                    nest_offl%u_right(1,k,j),  &
1039                                                    fac_dt ) *                 &
1040                             MERGE( 1.0_wp, 0.0_wp,                            &
1041                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1042             ENDDO
1043             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
1044          ENDDO
1045          DO  j = nys, nyn
1046             DO  k = nzb+1, nzt-1
1047                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
1048                                                    nest_offl%w_right(1,k,j),  &
1049                                                    fac_dt ) *                 &
1050                             MERGE( 1.0_wp, 0.0_wp,                            &
1051                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1052             ENDDO
[4079]1053             w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
[3347]1054          ENDDO
1055
1056          DO  j = nysv, nyn
1057             DO  k = nzb+1, nzt
1058                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
1059                                                    nest_offl%v_right(1,k,j),  &
1060                                                    fac_dt ) *                 &
1061                             MERGE( 1.0_wp, 0.0_wp,                            &
1062                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1063             ENDDO
1064             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
1065          ENDDO
1066
1067          IF ( .NOT. neutral )  THEN
1068             DO  j = nys, nyn
1069                DO  k = nzb+1, nzt
1070                   pt(k,j,nxr+1) = interpolate_in_time(                        &
1071                                                  nest_offl%pt_right(0,k,j),   &
1072                                                  nest_offl%pt_right(1,k,j),   &
1073                                                  fac_dt )
1074                ENDDO
1075                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
1076             ENDDO
1077          ENDIF
1078
1079          IF ( humidity )  THEN
1080             DO  j = nys, nyn
1081                DO  k = nzb+1, nzt
1082                   q(k,j,nxr+1) = interpolate_in_time(                         &
1083                                                  nest_offl%q_right(0,k,j),    &
1084                                                  nest_offl%q_right(1,k,j),    &
1085                                                  fac_dt ) 
1086 
1087                ENDDO
1088                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
1089             ENDDO
1090          ENDIF
[3737]1091         
1092          IF ( air_chemistry )  THEN
1093             DO  n = 1, UBOUND( chem_species, 1 )
1094                IF ( nest_offl%chem_from_file_r(n) )  THEN 
1095                   DO  j = nys, nyn
1096                      DO  k = nzb+1, nzt
1097                         chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
1098                                                 nest_offl%chem_right(0,k,j,n),&
1099                                                 nest_offl%chem_right(1,k,j,n),&
1100                                                 fac_dt                       )
1101                      ENDDO
1102                   ENDDO
1103                ENDIF
1104             ENDDO
1105          ENDIF
[3347]1106
1107       ENDIF
1108
1109       IF ( bc_dirichlet_s )  THEN
1110
1111          DO  i = nxl, nxr
1112             DO  k = nzb+1, nzt
1113                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
1114                                                nest_offl%v_south(1,k,i),      &
1115                                                fac_dt ) *                     &
1116                           MERGE( 1.0_wp, 0.0_wp,                              &
1117                                  BTEST( wall_flags_0(k,0,i), 2 ) )
1118                v(k,-1,i) = v(k,0,i) 
1119             ENDDO
1120             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
1121          ENDDO
1122
1123          DO  i = nxl, nxr
1124             DO  k = nzb+1, nzt-1
1125                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
1126                                                 nest_offl%w_south(1,k,i),     &
1127                                                 fac_dt ) *                    &
1128                           MERGE( 1.0_wp, 0.0_wp,                              &
1129                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
1130             ENDDO
[4079]1131             w(nzt,-1,i) = w(nzt-1,-1,i)
[3347]1132          ENDDO
1133
1134          DO  i = nxlu, nxr
1135             DO  k = nzb+1, nzt
1136                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
1137                                                 nest_offl%u_south(1,k,i),     &
1138                                                 fac_dt ) *                    &
1139                           MERGE( 1.0_wp, 0.0_wp,                              &
1140                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
1141             ENDDO
1142             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
1143          ENDDO
1144
1145          IF ( .NOT. neutral )  THEN
1146             DO  i = nxl, nxr
1147                DO  k = nzb+1, nzt
1148                   pt(k,-1,i) = interpolate_in_time(                           &
1149                                                 nest_offl%pt_south(0,k,i),    &
1150                                                 nest_offl%pt_south(1,k,i),    &
1151                                                 fac_dt )
1152 
1153                ENDDO
1154                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
1155             ENDDO
1156          ENDIF
1157
1158          IF ( humidity )  THEN
1159             DO  i = nxl, nxr
1160                DO  k = nzb+1, nzt
1161                   q(k,-1,i) = interpolate_in_time(                            &
1162                                                 nest_offl%q_south(0,k,i),     &
1163                                                 nest_offl%q_south(1,k,i),     &
1164                                                 fac_dt )
1165 
1166                ENDDO
1167                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
1168             ENDDO
1169          ENDIF
[3737]1170         
1171          IF ( air_chemistry )  THEN
1172             DO  n = 1, UBOUND( chem_species, 1 )
1173                IF ( nest_offl%chem_from_file_s(n) )  THEN 
1174                   DO  i = nxl, nxr
1175                      DO  k = nzb+1, nzt
1176                         chem_species(n)%conc(k,-1,i) = interpolate_in_time(   &
1177                                                 nest_offl%chem_south(0,k,i,n),&
1178                                                 nest_offl%chem_south(1,k,i,n),&
1179                                                 fac_dt                    )
1180                      ENDDO
1181                   ENDDO
1182                ENDIF
1183             ENDDO
1184          ENDIF
[3347]1185
1186       ENDIF
1187
1188       IF ( bc_dirichlet_n )  THEN
1189
1190          DO  i = nxl, nxr
1191             DO  k = nzb+1, nzt
1192                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
1193                                                    nest_offl%v_north(1,k,i),  &
1194                                                    fac_dt ) *                 &
1195                               MERGE( 1.0_wp, 0.0_wp,                          &
1196                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1197             ENDDO
1198             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
1199          ENDDO
1200          DO  i = nxl, nxr
1201             DO  k = nzb+1, nzt-1
1202                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
1203                                                    nest_offl%w_north(1,k,i),  &
1204                                                    fac_dt ) *                 &
1205                               MERGE( 1.0_wp, 0.0_wp,                          &
1206                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1207             ENDDO
[4079]1208             w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
[3347]1209          ENDDO
1210
1211          DO  i = nxlu, nxr
1212             DO  k = nzb+1, nzt
1213                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
1214                                                    nest_offl%u_north(1,k,i),  &
1215                                                    fac_dt ) *                 &
1216                               MERGE( 1.0_wp, 0.0_wp,                          &
1217                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1218
1219             ENDDO
1220             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
1221          ENDDO
1222
1223          IF ( .NOT. neutral )  THEN
1224             DO  i = nxl, nxr
1225                DO  k = nzb+1, nzt
1226                   pt(k,nyn+1,i) = interpolate_in_time(                        &
1227                                                    nest_offl%pt_north(0,k,i), &
1228                                                    nest_offl%pt_north(1,k,i), &
1229                                                    fac_dt )
1230 
1231                ENDDO
1232                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
1233             ENDDO
1234          ENDIF
1235
1236          IF ( humidity )  THEN
1237             DO  i = nxl, nxr
1238                DO  k = nzb+1, nzt
1239                   q(k,nyn+1,i) = interpolate_in_time(                         &
1240                                                    nest_offl%q_north(0,k,i),  &
1241                                                    nest_offl%q_north(1,k,i),  &
1242                                                    fac_dt )
1243 
1244                ENDDO
1245                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
1246             ENDDO
1247          ENDIF
[3737]1248         
1249          IF ( air_chemistry )  THEN
1250             DO  n = 1, UBOUND( chem_species, 1 )
1251                IF ( nest_offl%chem_from_file_n(n) )  THEN 
1252                   DO  i = nxl, nxr
1253                      DO  k = nzb+1, nzt
1254                         chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
1255                                                 nest_offl%chem_north(0,k,i,n),&
1256                                                 nest_offl%chem_north(1,k,i,n),&
1257                                                 fac_dt                       )
1258                      ENDDO
1259                   ENDDO
1260                ENDIF
1261             ENDDO
1262          ENDIF
[3347]1263
1264       ENDIF
1265!
1266!--    Top boundary
1267       DO  i = nxlu, nxr
1268          DO  j = nys, nyn
1269             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
1270                                                 nest_offl%u_top(1,j,i),       &
1271                                                 fac_dt ) *                    &
1272                           MERGE( 1.0_wp, 0.0_wp,                              &
1273                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
1274             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
1275          ENDDO
1276       ENDDO
[3937]1277!
1278!--    For left boundary set boundary condition for u-component also at top
1279!--    grid point.
1280!--    Note, this has no effect on the numeric solution, only for data output.
1281       IF ( bc_dirichlet_l )  u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
[3347]1282
1283       DO  i = nxl, nxr
1284          DO  j = nysv, nyn
1285             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
1286                                                 nest_offl%v_top(1,j,i),       &
1287                                                 fac_dt ) *                    &
1288                           MERGE( 1.0_wp, 0.0_wp,                              &
1289                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
1290             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
1291          ENDDO
1292       ENDDO
[3937]1293!
1294!--    For south boundary set boundary condition for v-component also at top
1295!--    grid point.
1296!--    Note, this has no effect on the numeric solution, only for data output.
1297       IF ( bc_dirichlet_s )  v(nzt+1,nys,:) = v(nzt+1,nysv,:)
[3347]1298
1299       DO  i = nxl, nxr
1300          DO  j = nys, nyn
1301             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
1302                                               nest_offl%w_top(1,j,i),         &
1303                                               fac_dt ) *                      &
1304                           MERGE( 1.0_wp, 0.0_wp,                              &
1305                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
1306             w(nzt+1,j,i) = w(nzt,j,i)
1307          ENDDO
1308       ENDDO
1309
1310
1311       IF ( .NOT. neutral )  THEN
1312          DO  i = nxl, nxr
1313             DO  j = nys, nyn
1314                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
1315                                                     nest_offl%pt_top(1,j,i),  &
1316                                                     fac_dt )
1317                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
1318             ENDDO
1319          ENDDO
1320       ENDIF
1321
1322       IF ( humidity )  THEN
1323          DO  i = nxl, nxr
1324             DO  j = nys, nyn
1325                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
1326                                                    nest_offl%q_top(1,j,i),    &
1327                                                    fac_dt )
1328                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
1329             ENDDO
1330          ENDDO
1331       ENDIF
[3737]1332       
1333       IF ( air_chemistry )  THEN
1334          DO  n = 1, UBOUND( chem_species, 1 )
1335             IF ( nest_offl%chem_from_file_t(n) )  THEN 
1336                DO  i = nxl, nxr
1337                   DO  j = nys, nyn
1338                      chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
[4079]1339                                              nest_offl%chem_top(0,j,i,n),   &
1340                                              nest_offl%chem_top(1,j,i,n),   &
[3737]1341                                              fac_dt                       )
1342                   ENDDO
1343                ENDDO
1344             ENDIF
1345          ENDDO
1346       ENDIF
[3347]1347!
1348!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
1349!--    passive scalar, dissipation, and chemical species if required
1350       IF ( rans_mode  .AND.  rans_tke_e )  THEN
1351          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
1352          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
1353          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
1354          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
1355       ENDIF
[4079]1356!        IF ( .NOT. constant_diffusion )  THEN
1357!           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
1358!           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
1359!           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
1360!           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
1361!           e(nzt+1,:,:) = e(nzt,:,:)
1362!        ENDIF
1363!        IF ( passive_scalar )  THEN
1364!           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
1365!           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
1366!           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
1367!           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
1368!        ENDIF
[3347]1369
1370       CALL exchange_horiz( u, nbgp )
1371       CALL exchange_horiz( v, nbgp )
1372       CALL exchange_horiz( w, nbgp )
1373       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
1374       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
[3737]1375       IF ( air_chemistry )  THEN
1376          DO  n = 1, UBOUND( chem_species, 1 )
[3858]1377!
1378!--          Do local exchange only when necessary, i.e. when data is coming
1379!--          from dynamic file.
1380             IF ( nest_offl%chem_from_file_t(n) )                              &
1381                CALL exchange_horiz( chem_species(n)%conc, nbgp )
[3737]1382          ENDDO
1383       ENDIF
[3347]1384!
[4079]1385!--    Set top boundary condition at all horizontal grid points, also at the
1386!--    lateral boundary grid points.
1387       w(nzt+1,:,:) = w(nzt,:,:)       
1388!
[3347]1389!--    In case of Rayleigh damping, where the profiles u_init, v_init
1390!--    q_init and pt_init are still used, update these profiles from the
1391!--    averaged boundary data.
1392!--    But first, average these data.
1393#if defined( __parallel )
1394       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
1395                           comm2d, ierr )
1396       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
1397                           comm2d, ierr )
1398       IF ( humidity )  THEN
1399          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
1400                              comm2d, ierr )
1401       ENDIF
1402       IF ( .NOT. neutral )  THEN
1403          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
1404                              comm2d, ierr )
1405       ENDIF
1406#else
[3704]1407       u_ref  = u_ref_l
1408       v_ref  = v_ref_l
1409       IF ( humidity )       q_ref  = q_ref_l
1410       IF ( .NOT. neutral )  pt_ref = pt_ref_l
[3347]1411#endif
1412!
[3704]1413!--    Average data. Note, reference profiles up to nzt are derived from lateral
1414!--    boundaries, at the model top it is derived from the top boundary. Thus,
1415!--    number of input data is different from nzb:nzt compared to nzt+1.
1416!--    Derived from lateral boundaries.
1417       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
1418                                               KIND = wp ) 
1419       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
1420                                               KIND = wp )
1421       IF ( humidity )                                                         &
1422          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
1423                                                          ( ny + 1 + nx + 1 ), &
1424                                                    KIND = wp )
1425       IF ( .NOT. neutral )                                                    &
1426          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
1427                                                          ( ny + 1 + nx + 1 ), &
1428                                              KIND = wp )
[3347]1429!
[3704]1430!--    Derived from top boundary.   
1431       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
1432       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
1433       IF ( humidity )                                                         &
1434          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
1435                                                KIND = wp )
1436       IF ( .NOT. neutral )                                                    &
1437          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
1438                                                KIND = wp )
[3347]1439!
[3704]1440!--    Write onto init profiles, which are used for damping
1441       u_init = u_ref
1442       v_init = v_ref
1443       IF ( humidity      )  q_init  = q_ref
1444       IF ( .NOT. neutral )  pt_init = pt_ref
[3347]1445!
[3704]1446!--    Set bottom boundary condition
1447       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
1448       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
[3347]1449!
[3704]1450!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
1451!--    Therefore, calculate boundary-layer depth first.
[4022]1452       CALL nesting_offl_calc_zi
[3704]1453       CALL adjust_sponge_layer 
[4226]1454       
1455       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
1456
1457       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'end' )
1458
1459
1460    END SUBROUTINE nesting_offl_bc
1461
1462!------------------------------------------------------------------------------!
1463! Description:
1464!------------------------------------------------------------------------------!
1465!>  Update of the geostrophic wind components.
1466!>  @todo: update geostrophic wind also in the child domains (should be done
1467!>         in the nesting.
1468!------------------------------------------------------------------------------!
1469    SUBROUTINE nesting_offl_geostrophic_wind
1470
1471       INTEGER(iwp) ::  k
[3347]1472!
[3704]1473!--    Update geostrophic wind components from dynamic input file.
1474       DO  k = nzb+1, nzt
1475          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
1476                                       fac_dt )
1477          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
1478                                       fac_dt )
1479       ENDDO
1480       ug(nzt+1) = ug(nzt)
1481       vg(nzt+1) = vg(nzt)
[3347]1482
[4226]1483    END SUBROUTINE nesting_offl_geostrophic_wind
[3987]1484
[4226]1485!------------------------------------------------------------------------------!
1486! Description:
1487!------------------------------------------------------------------------------!
1488!>  Determine the interpolation constant for time interpolation. The
1489!>  calculation is separated from the nesting_offl_bc and
1490!>  nesting_offl_geostrophic_wind in order to be independent on the order
1491!>  of calls.
1492!------------------------------------------------------------------------------!
1493    SUBROUTINE nesting_offl_interpolation_factor
1494!
1495!--    Determine interpolation factor and limit it to 1. This is because
1496!--    t+dt can slightly exceed time(tind_p) before boundary data is updated
1497!--    again.
1498       fac_dt = ( time_utc_init + time_since_reference_point                   &
1499                - nest_offl%time(nest_offl%tind) + dt_3d ) /                   &
1500           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
[3987]1501
[4226]1502       fac_dt = MIN( 1.0_wp, fac_dt )
[3347]1503
[4226]1504    END SUBROUTINE nesting_offl_interpolation_factor
1505
[3347]1506!------------------------------------------------------------------------------!
1507! Description:
1508!------------------------------------------------------------------------------!
1509!> Calculates the boundary-layer depth from the boundary data, according to
1510!> bulk-Richardson criterion.
1511!------------------------------------------------------------------------------!
[4022]1512    SUBROUTINE nesting_offl_calc_zi
[3347]1513
[4022]1514       INTEGER(iwp) :: i                            !< loop index in x-direction
1515       INTEGER(iwp) :: j                            !< loop index in y-direction
1516       INTEGER(iwp) :: k                            !< loop index in z-direction
1517       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
1518       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
1519       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
1520       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
[3347]1521   
1522       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
1523       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
1524       REAL(wp) ::  topo_max                !< maximum topography level in model domain
1525       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
1526       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
1527       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
1528       REAL(wp) ::  zi_local                !< local boundary-layer depth 
1529   
1530       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
[4022]1531       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
[3347]1532
1533       
1534!
1535!--    Calculate mean boundary-layer height from boundary data.
1536!--    Start with the left and right boundaries.
1537       zi_l      = 0.0_wp
[4022]1538       num_boundary_gp_non_cyclic_l = 0
[3347]1539       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
1540!
[4022]1541!--       Sum-up and store number of boundary grid points used for averaging
1542!--       ABL depth
1543          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
1544                                         nxr - nxl + 1
1545!
[3347]1546!--       Determine index along x. Please note, index indicates boundary
1547!--       grid point for scalars.
1548          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
1549   
1550          DO  j = nys, nyn
1551!
1552!--          Determine topography top index at current (j,i) index
[4168]1553             k_surface = topo_top_ind(j,i,0)
[3347]1554!
1555!--          Pre-compute surface virtual temperature. Therefore, use 2nd
1556!--          prognostic level according to Heinze et al. (2017).
1557             IF ( humidity )  THEN
1558                vpt_surface = pt(k_surface+2,j,i) *                            &
1559                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1560                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1561             ELSE
1562                vpt_surface = pt(k_surface+2,j,i)
1563                vpt_col     = pt(:,j,i)
1564             ENDIF
1565!
1566!--          Calculate local boundary layer height from bulk Richardson number,
1567!--          i.e. the height where the bulk Richardson number exceeds its
1568!--          critical value of 0.25 (according to Heinze et al., 2017).
1569!--          Note, no interpolation of u- and v-component is made, as both
1570!--          are mainly mean inflow profiles with very small spatial variation.
[3964]1571!--          Add a safety factor in case the velocity term becomes zero. This
1572!--          may happen if overhanging 3D structures are directly located at
1573!--          the boundary, where velocity inside the building is zero
1574!--          (k_surface is the index of the lowest upward-facing surface).
[4022]1575             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
1576                                      bc_dirichlet_l )**2 +                   &
1577                               v(:,j,i)**2 )
1578!
1579!--          Determine index of the maximum wind speed                               
1580             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1581
[3347]1582             zi_local = 0.0_wp
1583             DO  k = k_surface+1, nzt
1584                ri_bulk = zu(k) * g / vpt_surface *                            &
1585                          ( vpt_col(k) - vpt_surface ) /                       &
[4022]1586                          ( uv_abs(k) + 1E-5_wp ) 
1587!
1588!--             Check if critical Richardson number is exceeded. Further, check
1589!--             if there is a maxium in the wind profile in order to detect also
1590!--             ABL heights in the stable boundary layer.
1591                IF ( zi_local == 0.0_wp  .AND.                                 &
1592                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
[3347]1593                   zi_local = zu(k)           
1594             ENDDO
1595!
1596!--          Assure that the minimum local boundary-layer depth is at least at
1597!--          the second vertical grid level.
1598             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1599             
1600          ENDDO
1601       
1602       ENDIF
1603!
1604!--    Do the same at the north and south boundaries.
1605       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
1606       
[4022]1607          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
1608                                         nxr - nxl + 1
1609       
[3347]1610          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
1611       
1612          DO  i = nxl, nxr
[4168]1613             k_surface = topo_top_ind(j,i,0)
[3347]1614 
1615             IF ( humidity )  THEN
1616                vpt_surface = pt(k_surface+2,j,i) *                            &
1617                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1618                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1619             ELSE
1620                vpt_surface = pt(k_surface+2,j,i)
1621                vpt_col  = pt(:,j,i)
1622             ENDIF
1623         
[4022]1624             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
1625                               MERGE( v(:,j+1,i), v(:,j,i),                    &
1626                               bc_dirichlet_s )**2 )
1627!
1628!--          Determine index of the maximum wind speed
1629             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1630         
[3347]1631             zi_local = 0.0_wp
1632             DO  k = k_surface+1, nzt               
1633                ri_bulk = zu(k) * g / vpt_surface *                            &
1634                       ( vpt_col(k) - vpt_surface ) /                          &
[4022]1635                       ( uv_abs(k) + 1E-5_wp ) 
1636!
1637!--             Check if critical Richardson number is exceeded. Further, check
1638!--             if there is a maxium in the wind profile in order to detect also
1639!--             ABL heights in the stable boundary layer.                       
1640                IF ( zi_local == 0.0_wp  .AND.                                 &
1641                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
[3347]1642                   zi_local = zu(k)   
1643             ENDDO
1644             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1645         
1646          ENDDO
1647         
1648       ENDIF
1649   
1650#if defined( __parallel )
1651       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
1652                           comm2d, ierr )
[4022]1653       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
1654                           num_boundary_gp_non_cyclic,                         &
1655                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
[3347]1656#else
1657       zi_ribulk = zi_l
[4022]1658       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
[3347]1659#endif
[4022]1660       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
[3347]1661!
1662!--    Finally, check if boundary layer depth is not below the any topography.
1663!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
1664!--    lower level of the sponge layer, as well as to adjust the synthetic
1665!--    turbulence generator accordingly. If Rayleigh damping would be applied
1666!--    near buildings, etc., this would spoil the simulation results.
[4168]1667       topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
[3347]1668       
1669#if defined( __parallel )
[4022]1670       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
[3347]1671                           comm2d, ierr )
1672#else
1673       topo_max     = topo_max_l
1674#endif
[4022]1675!        zi_ribulk = MAX( zi_ribulk, topo_max )
[3937]1676       
[4022]1677    END SUBROUTINE nesting_offl_calc_zi
[3347]1678   
1679   
1680!------------------------------------------------------------------------------!
1681! Description:
1682!------------------------------------------------------------------------------!
1683!> Adjust the height where the rayleigh damping starts, i.e. the lower level
1684!> of the sponge layer.
1685!------------------------------------------------------------------------------!
1686    SUBROUTINE adjust_sponge_layer
1687
1688       INTEGER(iwp) :: k   !< loop index in z-direction
1689   
1690       REAL(wp) ::  rdh    !< updated Rayleigh damping height
1691 
1692   
1693       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
1694            rayleigh_damping_factor > 0.0_wp )  THEN
1695!
1696!--       Update Rayleigh-damping height and re-calculate height-depending
1697!--       damping coefficients.
1698!--       Assure that rayleigh damping starts well above the boundary layer.   
1699          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
1700                     0.8_wp * zu(nzt), rayleigh_damping_height )
1701!
1702!--       Update Rayleigh damping factor
1703          DO  k = nzb+1, nzt
1704             IF ( zu(k) >= rdh )  THEN
1705                rdf(k) = rayleigh_damping_factor *                             &
1706                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
1707                                        / ( zu(nzt) - rdh ) )                  &
1708                                          )**2
1709             ENDIF
1710          ENDDO
1711          rdf_sc = rdf
1712       
1713       ENDIF
1714
1715    END SUBROUTINE adjust_sponge_layer
1716   
1717!------------------------------------------------------------------------------!
1718! Description:
1719! ------------
1720!> Performs consistency checks
1721!------------------------------------------------------------------------------!
1722    SUBROUTINE nesting_offl_check_parameters 
1723!
[4169]1724!--    Check if offline nesting is applied in nested child domain.
[3579]1725       IF ( nesting_offline  .AND.  child_domain )  THEN
1726          message_string = 'Offline nesting is only applicable in root model.'
[4169]1727          CALL message( 'offline_nesting_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
[4226]1728       ENDIF
[3347]1729
1730    END SUBROUTINE nesting_offl_check_parameters
1731   
1732!------------------------------------------------------------------------------!
1733! Description:
1734! ------------
1735!> Reads the parameter list nesting_offl_parameters
1736!------------------------------------------------------------------------------!
1737    SUBROUTINE nesting_offl_parin 
1738       
1739       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
1740
1741
1742       NAMELIST /nesting_offl_parameters/   nesting_offline
1743
1744       line = ' '
1745
1746!
1747!--    Try to find stg package
1748       REWIND ( 11 )
1749       line = ' '
1750       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
1751          READ ( 11, '(A)', END=20 )  line
1752       ENDDO
1753       BACKSPACE ( 11 )
1754
1755!
1756!--    Read namelist
1757       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
1758
1759       GOTO 20
1760
1761 10    BACKSPACE( 11 )
1762       READ( 11 , '(A)') line
1763       CALL parin_fail_message( 'nesting_offl_parameters', line )
1764
1765 20    CONTINUE
1766
1767
1768    END SUBROUTINE nesting_offl_parin
1769
1770!------------------------------------------------------------------------------!
1771! Description:
1772! ------------
1773!> Writes information about offline nesting into HEADER file
1774!------------------------------------------------------------------------------!
1775    SUBROUTINE nesting_offl_header ( io )
1776
1777       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
1778
1779       WRITE ( io, 1 )
1780       IF ( nesting_offline )  THEN
1781          WRITE ( io, 3 )
1782       ELSE
1783          WRITE ( io, 2 )
1784       ENDIF
1785
17861 FORMAT (//' Offline nesting in COSMO model:'/                                &
1787              ' -------------------------------'/)
17882 FORMAT (' --> No offlince nesting is used (default) ')
17893 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
1790
1791    END SUBROUTINE nesting_offl_header 
1792
1793!------------------------------------------------------------------------------!
1794! Description:
1795! ------------
1796!> Allocate arrays used to read boundary data from NetCDF file and initialize
1797!> boundary data.
1798!------------------------------------------------------------------------------!
1799    SUBROUTINE nesting_offl_init
[4226]1800           
[3737]1801       INTEGER(iwp) ::  n !< running index for chemical species
[3347]1802
[4227]1803!
1804!--    Get time_utc_init from origin_date_time
1805       CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
[3347]1806
1807!--    Allocate arrays for geostrophic wind components. Arrays will
[3404]1808!--    incorporate 2 time levels in order to interpolate in between.
1809       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
1810       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
[3347]1811!
[4125]1812!--    Allocate arrays for reading left/right boundary values. Arrays will
1813!--    incorporate 2  time levels in order to interpolate in between. If the core has
1814!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1815!--    access. Dummy arrays will be allocated with dimension length zero.
[3347]1816       IF ( bc_dirichlet_l )  THEN
1817          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
1818          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
1819          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
1820          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
1821          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
[3737]1822          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
1823                                          1:UBOUND( chem_species, 1 )) )
[4125]1824       ELSE
1825          ALLOCATE( nest_offl%u_left(1:1,1:1,1:1)  )
1826          ALLOCATE( nest_offl%v_left(1:1,1:1,1:1)  )
1827          ALLOCATE( nest_offl%w_left(1:1,1:1,1:1)  )
1828          IF ( humidity )       ALLOCATE( nest_offl%q_left(1:1,1:1,1:1)  )
1829          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(1:1,1:1,1:1)  )
1830          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(1:1,1:1,1:1,     &
1831                                          1:UBOUND( chem_species, 1 )) )
[3347]1832       ENDIF
1833       IF ( bc_dirichlet_r )  THEN
1834          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
1835          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1836          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1837          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1838          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
[3737]1839          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
1840                                          1:UBOUND( chem_species, 1 )) )
[4125]1841       ELSE
1842          ALLOCATE( nest_offl%u_right(1:1,1:1,1:1)  )
1843          ALLOCATE( nest_offl%v_right(1:1,1:1,1:1)  )
1844          ALLOCATE( nest_offl%w_right(1:1,1:1,1:1)  )
1845          IF ( humidity )       ALLOCATE( nest_offl%q_right(1:1,1:1,1:1)  )
1846          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(1:1,1:1,1:1)  )
1847          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(1:1,1:1,1:1,    &
1848                                          1:UBOUND( chem_species, 1 )) )
[3347]1849       ENDIF
[4125]1850!
1851!--    Allocate arrays for reading north/south boundary values. Arrays will
1852!--    incorporate 2  time levels in order to interpolate in between. If the core has
1853!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1854!--    access. Dummy arrays will be allocated with dimension length zero.
[3347]1855       IF ( bc_dirichlet_n )  THEN
1856          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1857          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1858          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1859          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1860          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
[3737]1861          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
1862                                          1:UBOUND( chem_species, 1 )) )
[4125]1863       ELSE
1864          ALLOCATE( nest_offl%u_north(1:1,1:1,1:1)  )
1865          ALLOCATE( nest_offl%v_north(1:1,1:1,1:1)  )
1866          ALLOCATE( nest_offl%w_north(1:1,1:1,1:1)  )
1867          IF ( humidity )       ALLOCATE( nest_offl%q_north(1:1,1:1,1:1)  )
1868          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(1:1,1:1,1:1)  )
1869          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(1:1,1:1,1:1,    &
1870                                          1:UBOUND( chem_species, 1 )) )
[3347]1871       ENDIF
1872       IF ( bc_dirichlet_s )  THEN
1873          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1874          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1875          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1876          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1877          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
[3737]1878          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
1879                                          1:UBOUND( chem_species, 1 )) )
[4125]1880       ELSE
1881          ALLOCATE( nest_offl%u_south(1:1,1:1,1:1)  )
1882          ALLOCATE( nest_offl%v_south(1:1,1:1,1:1)  )
1883          ALLOCATE( nest_offl%w_south(1:1,1:1,1:1)  )
1884          IF ( humidity )       ALLOCATE( nest_offl%q_south(1:1,1:1,1:1)  )
1885          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(1:1,1:1,1:1)  )
1886          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(1:1,1:1,1:1,    &
1887                                          1:UBOUND( chem_species, 1 )) )
[3347]1888       ENDIF
[4125]1889!
1890!--    Allocate arrays for reading data at the top boundary. In contrast to the
1891!--    lateral boundaries, every core reads these data so that no dummy
1892!--    arrays need to be allocated.
[3347]1893       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1894       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1895       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1896       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1897       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
[3737]1898       IF ( air_chemistry )  ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
1899                                       1:UBOUND( chem_species, 1 )) )
[3347]1900!
[3737]1901!--    For chemical species, create the names of the variables. This is necessary
1902!--    to identify the respective variable and write it onto the correct array
1903!--    in the chem_species datatype.
1904       IF ( air_chemistry )  THEN
1905          ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
1906          ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
1907          ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
1908          ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
1909          ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
1910         
1911          ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
1912          ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
1913          ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
1914          ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
1915          ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
1916!
1917!--       Initialize flags that indicate whether the variable is on file or
1918!--       not. Please note, this is only necessary for chemistry variables.
1919          nest_offl%chem_from_file_l(:) = .FALSE.
1920          nest_offl%chem_from_file_n(:) = .FALSE.
1921          nest_offl%chem_from_file_r(:) = .FALSE.
1922          nest_offl%chem_from_file_s(:) = .FALSE.
1923          nest_offl%chem_from_file_t(:) = .FALSE.
1924         
1925          DO  n = 1, UBOUND( chem_species, 1 )
1926             nest_offl%var_names_chem_l(n) = nest_offl%char_l //               &
1927                                                  TRIM(chem_species(n)%name)
1928             nest_offl%var_names_chem_n(n) = nest_offl%char_n //               &
1929                                                  TRIM(chem_species(n)%name)
1930             nest_offl%var_names_chem_r(n) = nest_offl%char_r //               &
1931                                                  TRIM(chem_species(n)%name)
1932             nest_offl%var_names_chem_s(n) = nest_offl%char_s //               &
1933                                                  TRIM(chem_species(n)%name)
1934             nest_offl%var_names_chem_t(n) = nest_offl%char_t //               &
1935                                                  TRIM(chem_species(n)%name)
1936          ENDDO
1937       ENDIF
1938!
[4226]1939!--    Before initial data input is initiated, check if dynamic input file is
1940!--    present.
1941       IF ( .NOT. input_pids_dynamic )  THEN
1942          message_string = 'nesting_offline = .TRUE. requires dynamic '  //    &
1943                            'input file ' //                                   &
1944                            TRIM( input_file_dynamic ) // TRIM( coupling_char )
1945          CALL message( 'nesting_offl_init', 'PA0546', 1, 2, 0, 6, 0 )
1946       ENDIF
1947!
[3347]1948!--    Read COSMO data at lateral and top boundaries
[4226]1949       CALL nesting_offl_input
[3347]1950!
[4169]1951!--    Check if sufficient time steps are provided to cover the entire
1952!--    simulation. Note, dynamic input is only required for the 3D simulation,
1953!--    not for the soil/wall spinup. However, as the spinup time is added
1954!--    to the end_time, this must be considered here.
[4226]1955       IF ( end_time - spinup_time >                                           &
1956            nest_offl%time(nest_offl%nt-1) - time_utc_init )  THEN
1957          message_string = 'end_time of the simulation exceeds the ' //        &
1958                           'time dimension in the dynamic input file.'
1959          CALL message( 'nesting_offl_init', 'PA0183', 1, 2, 0, 6, 0 ) 
[4169]1960       ENDIF
[4226]1961
1962       IF ( nest_offl%time(0) /= time_utc_init )  THEN
1963          message_string = 'Offline nesting: time dimension must start at ' // &
1964                           ' time_utc_init.'
1965          CALL message( 'nesting_offl_init', 'PA0676', 1, 2, 0, 6, 0 )
1966       ENDIF
[4169]1967!
[3891]1968!--    Initialize boundary data. Please note, do not initialize boundaries in
1969!--    case of restart runs. This case the boundaries are already initialized
1970!--    and the boundary data from file would be on the wrong time level. 
1971       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1972          IF ( bc_dirichlet_l )  THEN
1973             u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
1974             v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
1975             w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
1976             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                  &
1977                                      nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
1978             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                  &
1979                                      nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
1980             IF ( air_chemistry )  THEN
1981                DO  n = 1, UBOUND( chem_species, 1 )
1982                   IF( nest_offl%chem_from_file_l(n) )  THEN
1983                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) =             &
1984                                      nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
1985                   ENDIF
1986                ENDDO
1987             ENDIF
[3737]1988          ENDIF
[3891]1989          IF ( bc_dirichlet_r )  THEN
1990             u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
1991             v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
1992             w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
1993             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =               &
1994                                      nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
1995             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =               &
1996                                      nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
1997             IF ( air_chemistry )  THEN
1998                DO  n = 1, UBOUND( chem_species, 1 )
1999                   IF( nest_offl%chem_from_file_r(n) )  THEN
2000                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) =          &
2001                                      nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
2002                   ENDIF
2003                ENDDO
2004             ENDIF
[3737]2005          ENDIF
[3891]2006          IF ( bc_dirichlet_s )  THEN
2007             u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
2008             v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
2009             w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
2010             IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                  &
2011                                      nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
2012             IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                  &
2013                                      nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
2014             IF ( air_chemistry )  THEN
2015                DO  n = 1, UBOUND( chem_species, 1 )
2016                   IF( nest_offl%chem_from_file_s(n) )  THEN
2017                      chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) =             &
2018                                      nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
2019                   ENDIF
2020                ENDDO
2021             ENDIF
[3737]2022          ENDIF
[3891]2023          IF ( bc_dirichlet_n )  THEN
2024             u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
2025             v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
2026             w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
2027             IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =               &
2028                                      nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
2029             IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =               &
2030                                      nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
2031             IF ( air_chemistry )  THEN
2032                DO  n = 1, UBOUND( chem_species, 1 )
2033                   IF( nest_offl%chem_from_file_n(n) )  THEN
2034                      chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =          &
2035                                      nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
2036                   ENDIF
2037                ENDDO
2038             ENDIF
[3737]2039          ENDIF
[3891]2040!         
2041!--       Initialize geostrophic wind components. Actually this is already done in
2042!--       init_3d_model when initializing_action = 'inifor', however, in speical
2043!--       case of user-defined initialization this will be done here again, in
2044!--       order to have a consistent initialization.
2045          ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
2046          vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
2047!         
2048!--       Set bottom and top boundary condition for geostrophic wind components
2049          ug(nzt+1) = ug(nzt)
2050          vg(nzt+1) = vg(nzt)
2051          ug(nzb)   = ug(nzb+1)
2052          vg(nzb)   = vg(nzb+1)
2053       ENDIF     
[3347]2054!
2055!--    After boundary data is initialized, mask topography at the
2056!--    boundaries for the velocity components.
2057       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
2058       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
2059       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
[3891]2060!
2061!--    Initial calculation of the boundary layer depth from the prescribed
2062!--    boundary data. This is requiered for initialize the synthetic turbulence
2063!--    generator correctly.
[4022]2064       CALL nesting_offl_calc_zi
[3347]2065       
2066!
[3891]2067!--    After boundary data is initialized, ensure mass conservation. Not
2068!--    necessary in restart runs.
2069       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2070          CALL nesting_offl_mass_conservation
2071       ENDIF
[3347]2072
2073    END SUBROUTINE nesting_offl_init
2074   
2075!------------------------------------------------------------------------------!
2076! Description:
2077!------------------------------------------------------------------------------!
2078!> Interpolation function, used to interpolate boundary data in time.
2079!------------------------------------------------------------------------------!
2080    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
2081
2082       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
2083       REAL(wp)            :: var_t1              !< boundary value at t1
2084       REAL(wp)            :: var_t2              !< boundary value at t2
2085       REAL(wp)            :: fac                 !< interpolation factor
2086
2087       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
2088
2089    END FUNCTION interpolate_in_time
2090
2091
2092
2093 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.