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

Last change on this file since 4282 was 4273, checked in by monakurppa, 5 years ago

Add logical switched nesting_chem and nesting_offline_chem

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