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

Last change on this file since 4079 was 4079, checked in by suehring, 2 years ago

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

  • Property svn:keywords set to Id
File size: 61.8 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 4079 2019-07-09 18:04:41Z suehring $
27! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
28!   though these won't enter the numerical solution. However, due to the mass
29!   conservation these values might some up to very large values which will
30!   occur in the run-control file
31! - Bugfix in offline nesting of chemical species
32! - Do not set Neumann conditions for TKE and passive scalar
33!
34! 4022 2019-06-12 11:52:39Z suehring
35! Detection of boundary-layer depth in stable boundary layer on basis of
36! boundary data improved
37! Routine for boundary-layer depth calculation renamed and made public
38!
39! 3987 2019-05-22 09:52:13Z kanani
40! Introduce alternative switch for debug output during timestepping
41!
42! 3964 2019-05-09 09:48:32Z suehring
43! Ensure that veloctiy term in calculation of bulk Richardson number does not
44! become zero
45!
46! 3937 2019-04-29 15:09:07Z suehring
47! Set boundary conditon on upper-left and upper-south grid point for the u- and
48! v-component, respectively.
49!
50! 3891 2019-04-12 17:52:01Z suehring
51! Bugfix, do not overwrite lateral and top boundary data in case of restart
52! runs.
53!
54! 3885 2019-04-11 11:29:34Z kanani
55! Changes related to global restructuring of location messages and introduction
56! of additional debug messages
57!
58!
59! Do local data exchange for chemistry variables only when boundary data is 
60! coming from dynamic file
61!
62! 3737 2019-02-12 16:57:06Z suehring
63! Introduce mesoscale nesting for chemical species
64!
65! 3705 2019-01-29 19:56:39Z suehring
66! Formatting adjustments
67!
68! 3704 2019-01-29 19:51:41Z suehring
69! Check implemented for offline nesting in child domain
70!
71! 3413 2018-10-24 10:28:44Z suehring
72! Keyword ID set
73!
74! 3404 2018-10-23 13:29:11Z suehring
75! Consider time-dependent geostrophic wind components in offline nesting
76!
77!
78! Initial Revision:
79! - separate offline nesting from large_scale_nudging_mod
80! - revise offline nesting, adjust for usage of synthetic turbulence generator
81! - adjust Rayleigh damping depending on the time-depending atmospheric
82!   conditions
83!
84!
85! Description:
86! ------------
87!> Offline nesting in larger-scale models. Boundary conditions for the simulation
88!> are read from NetCDF file and are prescribed onto the respective arrays.
89!> Further, a mass-flux correction is performed to maintain the mass balance.
90!--------------------------------------------------------------------------------!
91 MODULE nesting_offl_mod
92
93    USE arrays_3d,                                                             &
94        ONLY:  dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v,      &
95               v_init, vg, w, zu, zw
96                 
97    USE chem_modules,                                                          &
98        ONLY:  chem_species
99
100    USE control_parameters,                                                    &
101        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
102               bc_dirichlet_s, dt_3d, dz, constant_diffusion,                  &
103               debug_output_timestep,                                          &
104               humidity,                                                       &
105               initializing_actions,                                           &
106               message_string, nesting_offline, neutral, passive_scalar,       &
107               rans_mode, rans_tke_e, time_since_reference_point, volume_flow
108               
109    USE cpulog,                                                                &
110        ONLY:  cpu_log, log_point
111
112    USE grid_variables
113
114    USE indices,                                                               &
115        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
116               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
117
118    USE kinds
119
120    USE netcdf_data_input_mod,                                                 &
121        ONLY:  nest_offl
122       
123    USE pegrid
124
125    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
126                                     !< bulk Richardson number of 0.25
127   
128    SAVE
129    PRIVATE
130!
131!-- Public subroutines
132    PUBLIC nesting_offl_bc,                                                    &
133           nesting_offl_calc_zi,                                               &
134           nesting_offl_check_parameters,                                      &
135           nesting_offl_header,                                                &
136           nesting_offl_init,                                                  &
137           nesting_offl_mass_conservation,                                     &
138           nesting_offl_parin
139!
140!-- Public variables
141    PUBLIC zi_ribulk   
142
143    INTERFACE nesting_offl_bc
144       MODULE PROCEDURE nesting_offl_bc
145    END INTERFACE nesting_offl_bc
146   
147    INTERFACE nesting_offl_calc_zi
148       MODULE PROCEDURE nesting_offl_calc_zi
149    END INTERFACE nesting_offl_calc_zi
150   
151    INTERFACE nesting_offl_check_parameters
152       MODULE PROCEDURE nesting_offl_check_parameters
153    END INTERFACE nesting_offl_check_parameters
154   
155    INTERFACE nesting_offl_header
156       MODULE PROCEDURE nesting_offl_header
157    END INTERFACE nesting_offl_header
158   
159    INTERFACE nesting_offl_init
160       MODULE PROCEDURE nesting_offl_init
161    END INTERFACE nesting_offl_init
162           
163    INTERFACE nesting_offl_mass_conservation
164       MODULE PROCEDURE nesting_offl_mass_conservation
165    END INTERFACE nesting_offl_mass_conservation
166   
167    INTERFACE nesting_offl_parin
168       MODULE PROCEDURE nesting_offl_parin
169    END INTERFACE nesting_offl_parin
170
171 CONTAINS
172
173
174!------------------------------------------------------------------------------!
175! Description:
176! ------------
177!> In this subroutine a constant mass within the model domain is guaranteed.
178!> Larger-scale models may be based on a compressible equation system, which is
179!> not consistent with PALMs incompressible equation system. In order to avoid
180!> a decrease or increase of mass during the simulation, non-divergent flow
181!> through the lateral and top boundaries is compensated by the vertical wind
182!> component at the top boundary.
183!------------------------------------------------------------------------------!
184    SUBROUTINE nesting_offl_mass_conservation
185
186       IMPLICIT NONE
187
188       INTEGER(iwp) ::  i !< grid index in x-direction
189       INTEGER(iwp) ::  j !< grid index in y-direction
190       INTEGER(iwp) ::  k !< grid index in z-direction
191
192       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
193       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
194       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
195
196
197       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
198
199       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
200       
201       volume_flow   = 0.0_wp
202       volume_flow_l = 0.0_wp
203
204       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
205
206       IF ( bc_dirichlet_l )  THEN
207          i = nxl
208          DO  j = nys, nyn
209             DO  k = nzb+1, nzt
210                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
211                                     * MERGE( 1.0_wp, 0.0_wp,                  &
212                                              BTEST( wall_flags_0(k,j,i), 1 ) )
213             ENDDO
214          ENDDO
215       ENDIF
216       IF ( bc_dirichlet_r )  THEN
217          i = nxr+1
218          DO  j = nys, nyn
219             DO  k = nzb+1, nzt
220                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
221                                     * MERGE( 1.0_wp, 0.0_wp,                  &
222                                              BTEST( wall_flags_0(k,j,i), 1 ) )
223             ENDDO
224          ENDDO
225       ENDIF
226       IF ( bc_dirichlet_s )  THEN
227          j = nys
228          DO  i = nxl, nxr
229             DO  k = nzb+1, nzt
230                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
231                                     * MERGE( 1.0_wp, 0.0_wp,                  &
232                                              BTEST( wall_flags_0(k,j,i), 2 ) )
233             ENDDO
234          ENDDO
235       ENDIF
236       IF ( bc_dirichlet_n )  THEN
237          j = nyn+1
238          DO  i = nxl, nxr
239             DO  k = nzb+1, nzt
240                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
241                                     * MERGE( 1.0_wp, 0.0_wp,                  &
242                                              BTEST( wall_flags_0(k,j,i), 2 ) )
243             ENDDO
244          ENDDO
245       ENDIF
246!
247!--    Top boundary
248       k = nzt
249       DO  i = nxl, nxr
250          DO  j = nys, nyn
251             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
252          ENDDO
253       ENDDO
254
255#if defined( __parallel )
256       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
257       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
258                           comm2d, ierr )
259#else
260       volume_flow = volume_flow_l
261#endif
262
263       w_correct = SUM( volume_flow ) * d_area_t
264
265       DO  i = nxl, nxr
266          DO  j = nys, nyn
267             DO  k = nzt, nzt + 1
268                w(k,j,i) = w(k,j,i) + w_correct
269             ENDDO
270          ENDDO
271       ENDDO
272       
273       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
274
275       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
276
277    END SUBROUTINE nesting_offl_mass_conservation
278
279
280!------------------------------------------------------------------------------!
281! Description:
282! ------------
283!> Set the lateral and top boundary conditions in case the PALM domain is
284!> nested offline in a mesoscale model. Further, average boundary data and
285!> determine mean profiles, further used for correct damping in the sponge
286!> layer.
287!------------------------------------------------------------------------------!
288    SUBROUTINE nesting_offl_bc                     
289
290       IMPLICIT NONE
291
292       INTEGER(iwp) ::  i !< running index x-direction
293       INTEGER(iwp) ::  j !< running index y-direction
294       INTEGER(iwp) ::  k !< running index z-direction
295       INTEGER(iwp) ::  n !< running index for chemical species
296
297       REAL(wp) ::  fac_dt   !< interpolation factor
298       
299       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
300       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
301       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
302       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
303       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
304       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
305       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
306       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
307       
308
309       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
310
311       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
312!
313!--    Set mean profiles, derived from boundary data, to zero
314       pt_ref   = 0.0_wp
315       q_ref    = 0.0_wp
316       u_ref    = 0.0_wp
317       v_ref    = 0.0_wp
318       
319       pt_ref_l = 0.0_wp
320       q_ref_l  = 0.0_wp
321       u_ref_l  = 0.0_wp
322       v_ref_l  = 0.0_wp
323!
324!--    Determine interpolation factor and limit it to 1. This is because
325!--    t+dt can slightly exceed time(tind_p) before boundary data is updated
326!--    again.
327       fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind)  &
328                + dt_3d ) /                                                    &
329           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
330       fac_dt = MIN( 1.0_wp, fac_dt )
331!
332!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
333!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
334!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
335!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
336!--    Please note, at the left (for u) and south (for v) boundary, values
337!--    for u and v are set also at i/j=-1, since these values are used in
338!--    boundary_conditions() to restore prognostic values.
339!--    Further, sum up data to calculate mean profiles from boundary data,
340!--    used for Rayleigh damping.
341       IF ( bc_dirichlet_l )  THEN
342
343          DO  j = nys, nyn
344             DO  k = nzb+1, nzt
345                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
346                                                nest_offl%u_left(1,k,j),       &
347                                                fac_dt ) *                     &
348                             MERGE( 1.0_wp, 0.0_wp,                            &
349                                    BTEST( wall_flags_0(k,j,0), 1 ) )
350                u(k,j,-1) = u(k,j,0)
351             ENDDO
352             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
353          ENDDO
354
355          DO  j = nys, nyn
356             DO  k = nzb+1, nzt-1
357                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
358                                                 nest_offl%w_left(1,k,j),      &
359                                                 fac_dt ) *                    &
360                            MERGE( 1.0_wp, 0.0_wp,                             &
361                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
362             ENDDO
363             w(nzt,j,-1) = w(nzt-1,j,-1)
364          ENDDO
365
366          DO  j = nysv, nyn
367             DO  k = nzb+1, nzt
368                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
369                                                 nest_offl%v_left(1,k,j),      &
370                                                 fac_dt ) *                    &
371                               MERGE( 1.0_wp, 0.0_wp,                          &
372                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
373             ENDDO
374             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
375          ENDDO
376
377          IF ( .NOT. neutral )  THEN
378             DO  j = nys, nyn
379                DO  k = nzb+1, nzt
380                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
381                                                     nest_offl%pt_left(1,k,j), &
382                                                     fac_dt )
383 
384                ENDDO
385                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
386             ENDDO
387          ENDIF
388
389          IF ( humidity )  THEN
390             DO  j = nys, nyn
391                DO  k = nzb+1, nzt
392                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
393                                                    nest_offl%q_left(1,k,j),   &
394                                                    fac_dt )
395 
396                ENDDO
397                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
398             ENDDO
399          ENDIF
400         
401          IF ( air_chemistry )  THEN
402             DO  n = 1, UBOUND( chem_species, 1 )
403                IF ( nest_offl%chem_from_file_l(n) )  THEN                   
404                   DO  j = nys, nyn
405                      DO  k = nzb+1, nzt
406                         chem_species(n)%conc(k,j,-1) = interpolate_in_time(   &
407                                                  nest_offl%chem_left(0,k,j,n),&
408                                                  nest_offl%chem_left(1,k,j,n),&
409                                                  fac_dt                   )
410                      ENDDO
411                   ENDDO
412                ENDIF
413             ENDDO
414          ENDIF
415
416       ENDIF
417
418       IF ( bc_dirichlet_r )  THEN
419
420          DO  j = nys, nyn
421             DO  k = nzb+1, nzt
422                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
423                                                    nest_offl%u_right(1,k,j),  &
424                                                    fac_dt ) *                 &
425                             MERGE( 1.0_wp, 0.0_wp,                            &
426                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
427             ENDDO
428             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
429          ENDDO
430          DO  j = nys, nyn
431             DO  k = nzb+1, nzt-1
432                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
433                                                    nest_offl%w_right(1,k,j),  &
434                                                    fac_dt ) *                 &
435                             MERGE( 1.0_wp, 0.0_wp,                            &
436                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
437             ENDDO
438             w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
439          ENDDO
440
441          DO  j = nysv, nyn
442             DO  k = nzb+1, nzt
443                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
444                                                    nest_offl%v_right(1,k,j),  &
445                                                    fac_dt ) *                 &
446                             MERGE( 1.0_wp, 0.0_wp,                            &
447                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
448             ENDDO
449             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
450          ENDDO
451
452          IF ( .NOT. neutral )  THEN
453             DO  j = nys, nyn
454                DO  k = nzb+1, nzt
455                   pt(k,j,nxr+1) = interpolate_in_time(                        &
456                                                  nest_offl%pt_right(0,k,j),   &
457                                                  nest_offl%pt_right(1,k,j),   &
458                                                  fac_dt )
459                ENDDO
460                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
461             ENDDO
462          ENDIF
463
464          IF ( humidity )  THEN
465             DO  j = nys, nyn
466                DO  k = nzb+1, nzt
467                   q(k,j,nxr+1) = interpolate_in_time(                         &
468                                                  nest_offl%q_right(0,k,j),    &
469                                                  nest_offl%q_right(1,k,j),    &
470                                                  fac_dt ) 
471 
472                ENDDO
473                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
474             ENDDO
475          ENDIF
476         
477          IF ( air_chemistry )  THEN
478             DO  n = 1, UBOUND( chem_species, 1 )
479                IF ( nest_offl%chem_from_file_r(n) )  THEN 
480                   DO  j = nys, nyn
481                      DO  k = nzb+1, nzt
482                         chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
483                                                 nest_offl%chem_right(0,k,j,n),&
484                                                 nest_offl%chem_right(1,k,j,n),&
485                                                 fac_dt                       )
486                      ENDDO
487                   ENDDO
488                ENDIF
489             ENDDO
490          ENDIF
491
492       ENDIF
493
494       IF ( bc_dirichlet_s )  THEN
495
496          DO  i = nxl, nxr
497             DO  k = nzb+1, nzt
498                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
499                                                nest_offl%v_south(1,k,i),      &
500                                                fac_dt ) *                     &
501                           MERGE( 1.0_wp, 0.0_wp,                              &
502                                  BTEST( wall_flags_0(k,0,i), 2 ) )
503                v(k,-1,i) = v(k,0,i) 
504             ENDDO
505             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
506          ENDDO
507
508          DO  i = nxl, nxr
509             DO  k = nzb+1, nzt-1
510                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
511                                                 nest_offl%w_south(1,k,i),     &
512                                                 fac_dt ) *                    &
513                           MERGE( 1.0_wp, 0.0_wp,                              &
514                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
515             ENDDO
516             w(nzt,-1,i) = w(nzt-1,-1,i)
517          ENDDO
518
519          DO  i = nxlu, nxr
520             DO  k = nzb+1, nzt
521                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
522                                                 nest_offl%u_south(1,k,i),     &
523                                                 fac_dt ) *                    &
524                           MERGE( 1.0_wp, 0.0_wp,                              &
525                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
526             ENDDO
527             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
528          ENDDO
529
530          IF ( .NOT. neutral )  THEN
531             DO  i = nxl, nxr
532                DO  k = nzb+1, nzt
533                   pt(k,-1,i) = interpolate_in_time(                           &
534                                                 nest_offl%pt_south(0,k,i),    &
535                                                 nest_offl%pt_south(1,k,i),    &
536                                                 fac_dt )
537 
538                ENDDO
539                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
540             ENDDO
541          ENDIF
542
543          IF ( humidity )  THEN
544             DO  i = nxl, nxr
545                DO  k = nzb+1, nzt
546                   q(k,-1,i) = interpolate_in_time(                            &
547                                                 nest_offl%q_south(0,k,i),     &
548                                                 nest_offl%q_south(1,k,i),     &
549                                                 fac_dt )
550 
551                ENDDO
552                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
553             ENDDO
554          ENDIF
555         
556          IF ( air_chemistry )  THEN
557             DO  n = 1, UBOUND( chem_species, 1 )
558                IF ( nest_offl%chem_from_file_s(n) )  THEN 
559                   DO  i = nxl, nxr
560                      DO  k = nzb+1, nzt
561                         chem_species(n)%conc(k,-1,i) = interpolate_in_time(   &
562                                                 nest_offl%chem_south(0,k,i,n),&
563                                                 nest_offl%chem_south(1,k,i,n),&
564                                                 fac_dt                    )
565                      ENDDO
566                   ENDDO
567                ENDIF
568             ENDDO
569          ENDIF
570
571       ENDIF
572
573       IF ( bc_dirichlet_n )  THEN
574
575          DO  i = nxl, nxr
576             DO  k = nzb+1, nzt
577                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
578                                                    nest_offl%v_north(1,k,i),  &
579                                                    fac_dt ) *                 &
580                               MERGE( 1.0_wp, 0.0_wp,                          &
581                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
582             ENDDO
583             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
584          ENDDO
585          DO  i = nxl, nxr
586             DO  k = nzb+1, nzt-1
587                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
588                                                    nest_offl%w_north(1,k,i),  &
589                                                    fac_dt ) *                 &
590                               MERGE( 1.0_wp, 0.0_wp,                          &
591                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
592             ENDDO
593             w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
594          ENDDO
595
596          DO  i = nxlu, nxr
597             DO  k = nzb+1, nzt
598                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
599                                                    nest_offl%u_north(1,k,i),  &
600                                                    fac_dt ) *                 &
601                               MERGE( 1.0_wp, 0.0_wp,                          &
602                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
603
604             ENDDO
605             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
606          ENDDO
607
608          IF ( .NOT. neutral )  THEN
609             DO  i = nxl, nxr
610                DO  k = nzb+1, nzt
611                   pt(k,nyn+1,i) = interpolate_in_time(                        &
612                                                    nest_offl%pt_north(0,k,i), &
613                                                    nest_offl%pt_north(1,k,i), &
614                                                    fac_dt )
615 
616                ENDDO
617                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
618             ENDDO
619          ENDIF
620
621          IF ( humidity )  THEN
622             DO  i = nxl, nxr
623                DO  k = nzb+1, nzt
624                   q(k,nyn+1,i) = interpolate_in_time(                         &
625                                                    nest_offl%q_north(0,k,i),  &
626                                                    nest_offl%q_north(1,k,i),  &
627                                                    fac_dt )
628 
629                ENDDO
630                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
631             ENDDO
632          ENDIF
633         
634          IF ( air_chemistry )  THEN
635             DO  n = 1, UBOUND( chem_species, 1 )
636                IF ( nest_offl%chem_from_file_n(n) )  THEN 
637                   DO  i = nxl, nxr
638                      DO  k = nzb+1, nzt
639                         chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
640                                                 nest_offl%chem_north(0,k,i,n),&
641                                                 nest_offl%chem_north(1,k,i,n),&
642                                                 fac_dt                       )
643                      ENDDO
644                   ENDDO
645                ENDIF
646             ENDDO
647          ENDIF
648
649       ENDIF
650!
651!--    Top boundary
652       DO  i = nxlu, nxr
653          DO  j = nys, nyn
654             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
655                                                 nest_offl%u_top(1,j,i),       &
656                                                 fac_dt ) *                    &
657                           MERGE( 1.0_wp, 0.0_wp,                              &
658                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
659             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
660          ENDDO
661       ENDDO
662!
663!--    For left boundary set boundary condition for u-component also at top
664!--    grid point.
665!--    Note, this has no effect on the numeric solution, only for data output.
666       IF ( bc_dirichlet_l )  u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
667
668       DO  i = nxl, nxr
669          DO  j = nysv, nyn
670             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
671                                                 nest_offl%v_top(1,j,i),       &
672                                                 fac_dt ) *                    &
673                           MERGE( 1.0_wp, 0.0_wp,                              &
674                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
675             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
676          ENDDO
677       ENDDO
678!
679!--    For south boundary set boundary condition for v-component also at top
680!--    grid point.
681!--    Note, this has no effect on the numeric solution, only for data output.
682       IF ( bc_dirichlet_s )  v(nzt+1,nys,:) = v(nzt+1,nysv,:)
683
684       DO  i = nxl, nxr
685          DO  j = nys, nyn
686             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
687                                               nest_offl%w_top(1,j,i),         &
688                                               fac_dt ) *                      &
689                           MERGE( 1.0_wp, 0.0_wp,                              &
690                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
691             w(nzt+1,j,i) = w(nzt,j,i)
692          ENDDO
693       ENDDO
694
695
696       IF ( .NOT. neutral )  THEN
697          DO  i = nxl, nxr
698             DO  j = nys, nyn
699                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
700                                                     nest_offl%pt_top(1,j,i),  &
701                                                     fac_dt )
702                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
703             ENDDO
704          ENDDO
705       ENDIF
706
707       IF ( humidity )  THEN
708          DO  i = nxl, nxr
709             DO  j = nys, nyn
710                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
711                                                    nest_offl%q_top(1,j,i),    &
712                                                    fac_dt )
713                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
714             ENDDO
715          ENDDO
716       ENDIF
717       
718       IF ( air_chemistry )  THEN
719          DO  n = 1, UBOUND( chem_species, 1 )
720             IF ( nest_offl%chem_from_file_t(n) )  THEN 
721                DO  i = nxl, nxr
722                   DO  j = nys, nyn
723                      chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
724                                              nest_offl%chem_top(0,j,i,n),   &
725                                              nest_offl%chem_top(1,j,i,n),   &
726                                              fac_dt                       )
727                   ENDDO
728                ENDDO
729             ENDIF
730          ENDDO
731       ENDIF
732!
733!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
734!--    passive scalar, dissipation, and chemical species if required
735       IF ( rans_mode  .AND.  rans_tke_e )  THEN
736          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
737          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
738          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
739          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
740       ENDIF
741!        IF ( .NOT. constant_diffusion )  THEN
742!           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
743!           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
744!           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
745!           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
746!           e(nzt+1,:,:) = e(nzt,:,:)
747!        ENDIF
748!        IF ( passive_scalar )  THEN
749!           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
750!           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
751!           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
752!           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
753!        ENDIF
754
755       CALL exchange_horiz( u, nbgp )
756       CALL exchange_horiz( v, nbgp )
757       CALL exchange_horiz( w, nbgp )
758       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
759       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
760       IF ( air_chemistry )  THEN
761          DO  n = 1, UBOUND( chem_species, 1 )
762!
763!--          Do local exchange only when necessary, i.e. when data is coming
764!--          from dynamic file.
765             IF ( nest_offl%chem_from_file_t(n) )                              &
766                CALL exchange_horiz( chem_species(n)%conc, nbgp )
767          ENDDO
768       ENDIF
769!
770!--    Set top boundary condition at all horizontal grid points, also at the
771!--    lateral boundary grid points.
772       w(nzt+1,:,:) = w(nzt,:,:)       
773!
774!--    In case of Rayleigh damping, where the profiles u_init, v_init
775!--    q_init and pt_init are still used, update these profiles from the
776!--    averaged boundary data.
777!--    But first, average these data.
778#if defined( __parallel )
779       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
780                           comm2d, ierr )
781       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
782                           comm2d, ierr )
783       IF ( humidity )  THEN
784          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
785                              comm2d, ierr )
786       ENDIF
787       IF ( .NOT. neutral )  THEN
788          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
789                              comm2d, ierr )
790       ENDIF
791#else
792       u_ref  = u_ref_l
793       v_ref  = v_ref_l
794       IF ( humidity )       q_ref  = q_ref_l
795       IF ( .NOT. neutral )  pt_ref = pt_ref_l
796#endif
797!
798!--    Average data. Note, reference profiles up to nzt are derived from lateral
799!--    boundaries, at the model top it is derived from the top boundary. Thus,
800!--    number of input data is different from nzb:nzt compared to nzt+1.
801!--    Derived from lateral boundaries.
802       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
803                                               KIND = wp ) 
804       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
805                                               KIND = wp )
806       IF ( humidity )                                                         &
807          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
808                                                          ( ny + 1 + nx + 1 ), &
809                                                    KIND = wp )
810       IF ( .NOT. neutral )                                                    &
811          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
812                                                          ( ny + 1 + nx + 1 ), &
813                                              KIND = wp )
814!
815!--    Derived from top boundary.   
816       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
817       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
818       IF ( humidity )                                                         &
819          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
820                                                KIND = wp )
821       IF ( .NOT. neutral )                                                    &
822          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
823                                                KIND = wp )
824!
825!--    Write onto init profiles, which are used for damping
826       u_init = u_ref
827       v_init = v_ref
828       IF ( humidity      )  q_init  = q_ref
829       IF ( .NOT. neutral )  pt_init = pt_ref
830!
831!--    Set bottom boundary condition
832       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
833       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
834!
835!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
836!--    Therefore, calculate boundary-layer depth first.
837       CALL nesting_offl_calc_zi
838       CALL adjust_sponge_layer
839   
840!
841!--    Update geostrophic wind components from dynamic input file.
842       DO  k = nzb+1, nzt
843          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
844                                       fac_dt )
845          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
846                                       fac_dt )
847       ENDDO
848       ug(nzt+1) = ug(nzt)
849       vg(nzt+1) = vg(nzt)
850   
851       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
852
853       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'end' )
854
855
856    END SUBROUTINE nesting_offl_bc
857
858!------------------------------------------------------------------------------!
859! Description:
860!------------------------------------------------------------------------------!
861!> Calculates the boundary-layer depth from the boundary data, according to
862!> bulk-Richardson criterion.
863!------------------------------------------------------------------------------!
864    SUBROUTINE nesting_offl_calc_zi
865       
866       USE basic_constants_and_equations_mod,                                  &
867           ONLY:  g
868       
869       USE kinds
870       
871       USE surface_mod,                                                        &
872           ONLY:  get_topography_top_index, get_topography_top_index_ji
873
874       IMPLICIT NONE
875
876       INTEGER(iwp) :: i                            !< loop index in x-direction
877       INTEGER(iwp) :: j                            !< loop index in y-direction
878       INTEGER(iwp) :: k                            !< loop index in z-direction
879       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
880       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
881       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
882       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
883   
884       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
885       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
886       REAL(wp) ::  topo_max                !< maximum topography level in model domain
887       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
888       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
889       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
890       REAL(wp) ::  zi_local                !< local boundary-layer depth 
891   
892       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
893       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
894
895       
896!
897!--    Calculate mean boundary-layer height from boundary data.
898!--    Start with the left and right boundaries.
899       zi_l      = 0.0_wp
900       num_boundary_gp_non_cyclic_l = 0
901       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
902!
903!--       Sum-up and store number of boundary grid points used for averaging
904!--       ABL depth
905          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
906                                         nxr - nxl + 1
907!
908!--       Determine index along x. Please note, index indicates boundary
909!--       grid point for scalars.
910          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
911   
912          DO  j = nys, nyn
913!
914!--          Determine topography top index at current (j,i) index
915             k_surface = get_topography_top_index_ji( j, i, 's' )
916!
917!--          Pre-compute surface virtual temperature. Therefore, use 2nd
918!--          prognostic level according to Heinze et al. (2017).
919             IF ( humidity )  THEN
920                vpt_surface = pt(k_surface+2,j,i) *                            &
921                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
922                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
923             ELSE
924                vpt_surface = pt(k_surface+2,j,i)
925                vpt_col     = pt(:,j,i)
926             ENDIF
927!
928!--          Calculate local boundary layer height from bulk Richardson number,
929!--          i.e. the height where the bulk Richardson number exceeds its
930!--          critical value of 0.25 (according to Heinze et al., 2017).
931!--          Note, no interpolation of u- and v-component is made, as both
932!--          are mainly mean inflow profiles with very small spatial variation.
933!--          Add a safety factor in case the velocity term becomes zero. This
934!--          may happen if overhanging 3D structures are directly located at
935!--          the boundary, where velocity inside the building is zero
936!--          (k_surface is the index of the lowest upward-facing surface).
937             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
938                                      bc_dirichlet_l )**2 +                   &
939                               v(:,j,i)**2 )
940!
941!--          Determine index of the maximum wind speed                               
942             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
943
944             zi_local = 0.0_wp
945             DO  k = k_surface+1, nzt
946                ri_bulk = zu(k) * g / vpt_surface *                            &
947                          ( vpt_col(k) - vpt_surface ) /                       &
948                          ( uv_abs(k) + 1E-5_wp ) 
949!
950!--             Check if critical Richardson number is exceeded. Further, check
951!--             if there is a maxium in the wind profile in order to detect also
952!--             ABL heights in the stable boundary layer.
953                IF ( zi_local == 0.0_wp  .AND.                                 &
954                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
955                   zi_local = zu(k)           
956             ENDDO
957!
958!--          Assure that the minimum local boundary-layer depth is at least at
959!--          the second vertical grid level.
960             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
961             
962          ENDDO
963       
964       ENDIF
965!
966!--    Do the same at the north and south boundaries.
967       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
968       
969          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
970                                         nxr - nxl + 1
971       
972          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
973       
974          DO  i = nxl, nxr
975             k_surface = get_topography_top_index_ji( j, i, 's' )
976 
977             IF ( humidity )  THEN
978                vpt_surface = pt(k_surface+2,j,i) *                            &
979                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
980                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
981             ELSE
982                vpt_surface = pt(k_surface+2,j,i)
983                vpt_col  = pt(:,j,i)
984             ENDIF
985         
986             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
987                               MERGE( v(:,j+1,i), v(:,j,i),                    &
988                               bc_dirichlet_s )**2 )
989!
990!--          Determine index of the maximum wind speed
991             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
992         
993             zi_local = 0.0_wp
994             DO  k = k_surface+1, nzt               
995                ri_bulk = zu(k) * g / vpt_surface *                            &
996                       ( vpt_col(k) - vpt_surface ) /                          &
997                       ( uv_abs(k) + 1E-5_wp ) 
998!
999!--             Check if critical Richardson number is exceeded. Further, check
1000!--             if there is a maxium in the wind profile in order to detect also
1001!--             ABL heights in the stable boundary layer.                       
1002                IF ( zi_local == 0.0_wp  .AND.                                 &
1003                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
1004                   zi_local = zu(k)   
1005             ENDDO
1006             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1007         
1008          ENDDO
1009         
1010       ENDIF
1011   
1012#if defined( __parallel )
1013       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
1014                           comm2d, ierr )
1015       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
1016                           num_boundary_gp_non_cyclic,                         &
1017                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1018#else
1019       zi_ribulk = zi_l
1020       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
1021#endif
1022       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
1023!
1024!--    Finally, check if boundary layer depth is not below the any topography.
1025!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
1026!--    lower level of the sponge layer, as well as to adjust the synthetic
1027!--    turbulence generator accordingly. If Rayleigh damping would be applied
1028!--    near buildings, etc., this would spoil the simulation results.
1029       topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
1030       
1031#if defined( __parallel )
1032       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
1033                           comm2d, ierr )
1034#else
1035       topo_max     = topo_max_l
1036#endif
1037!        zi_ribulk = MAX( zi_ribulk, topo_max )
1038       
1039    END SUBROUTINE nesting_offl_calc_zi
1040   
1041   
1042!------------------------------------------------------------------------------!
1043! Description:
1044!------------------------------------------------------------------------------!
1045!> Adjust the height where the rayleigh damping starts, i.e. the lower level
1046!> of the sponge layer.
1047!------------------------------------------------------------------------------!
1048    SUBROUTINE adjust_sponge_layer
1049       
1050       USE arrays_3d,                                                          &
1051           ONLY:  rdf, rdf_sc, zu
1052       
1053       USE basic_constants_and_equations_mod,                                  &
1054           ONLY:  pi
1055           
1056       USE control_parameters,                                                 &
1057           ONLY:  rayleigh_damping_height, rayleigh_damping_factor
1058       
1059       USE kinds
1060
1061       IMPLICIT NONE
1062
1063       INTEGER(iwp) :: k   !< loop index in z-direction
1064   
1065       REAL(wp) ::  rdh    !< updated Rayleigh damping height
1066 
1067   
1068       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
1069            rayleigh_damping_factor > 0.0_wp )  THEN
1070!
1071!--       Update Rayleigh-damping height and re-calculate height-depending
1072!--       damping coefficients.
1073!--       Assure that rayleigh damping starts well above the boundary layer.   
1074          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
1075                     0.8_wp * zu(nzt), rayleigh_damping_height )
1076!
1077!--       Update Rayleigh damping factor
1078          DO  k = nzb+1, nzt
1079             IF ( zu(k) >= rdh )  THEN
1080                rdf(k) = rayleigh_damping_factor *                             &
1081                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
1082                                        / ( zu(nzt) - rdh ) )                  &
1083                                          )**2
1084             ENDIF
1085          ENDDO
1086          rdf_sc = rdf
1087       
1088       ENDIF
1089
1090    END SUBROUTINE adjust_sponge_layer
1091   
1092!------------------------------------------------------------------------------!
1093! Description:
1094! ------------
1095!> Performs consistency checks
1096!------------------------------------------------------------------------------!
1097    SUBROUTINE nesting_offl_check_parameters
1098   
1099       USE control_parameters,                                                    &
1100        ONLY:  child_domain, message_string, nesting_offline
1101
1102       IMPLICIT NONE
1103!
1104!--    Perform checks
1105       IF ( nesting_offline  .AND.  child_domain )  THEN
1106          message_string = 'Offline nesting is only applicable in root model.'
1107          CALL message( 'stg_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
1108       ENDIF
1109
1110
1111    END SUBROUTINE nesting_offl_check_parameters
1112   
1113!------------------------------------------------------------------------------!
1114! Description:
1115! ------------
1116!> Reads the parameter list nesting_offl_parameters
1117!------------------------------------------------------------------------------!
1118    SUBROUTINE nesting_offl_parin
1119
1120       IMPLICIT NONE
1121       
1122       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
1123
1124
1125       NAMELIST /nesting_offl_parameters/   nesting_offline
1126
1127       line = ' '
1128
1129!
1130!--    Try to find stg package
1131       REWIND ( 11 )
1132       line = ' '
1133       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
1134          READ ( 11, '(A)', END=20 )  line
1135       ENDDO
1136       BACKSPACE ( 11 )
1137
1138!
1139!--    Read namelist
1140       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
1141
1142       GOTO 20
1143
1144 10    BACKSPACE( 11 )
1145       READ( 11 , '(A)') line
1146       CALL parin_fail_message( 'nesting_offl_parameters', line )
1147
1148 20    CONTINUE
1149
1150
1151    END SUBROUTINE nesting_offl_parin
1152
1153!------------------------------------------------------------------------------!
1154! Description:
1155! ------------
1156!> Writes information about offline nesting into HEADER file
1157!------------------------------------------------------------------------------!
1158    SUBROUTINE nesting_offl_header ( io )
1159
1160       IMPLICIT NONE
1161
1162       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
1163
1164       WRITE ( io, 1 )
1165       IF ( nesting_offline )  THEN
1166          WRITE ( io, 3 )
1167       ELSE
1168          WRITE ( io, 2 )
1169       ENDIF
1170
11711 FORMAT (//' Offline nesting in COSMO model:'/                                &
1172              ' -------------------------------'/)
11732 FORMAT (' --> No offlince nesting is used (default) ')
11743 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
1175
1176    END SUBROUTINE nesting_offl_header
1177
1178!------------------------------------------------------------------------------!
1179! Description:
1180! ------------
1181!> Allocate arrays used to read boundary data from NetCDF file and initialize
1182!> boundary data.
1183!------------------------------------------------------------------------------!
1184    SUBROUTINE nesting_offl_init
1185   
1186       USE netcdf_data_input_mod,                                              &
1187           ONLY:  netcdf_data_input_offline_nesting 
1188
1189       IMPLICIT NONE
1190       
1191       INTEGER(iwp) ::  n !< running index for chemical species
1192
1193
1194!--    Allocate arrays for geostrophic wind components. Arrays will
1195!--    incorporate 2 time levels in order to interpolate in between.
1196       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
1197       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
1198!
1199!--    Allocate arrays for reading boundary values. Arrays will incorporate 2
1200!--    time levels in order to interpolate in between.
1201       IF ( bc_dirichlet_l )  THEN
1202          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
1203          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
1204          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
1205          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
1206          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
1207          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
1208                                          1:UBOUND( chem_species, 1 )) )
1209       ENDIF
1210       IF ( bc_dirichlet_r )  THEN
1211          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
1212          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1213          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1214          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1215          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1216          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
1217                                          1:UBOUND( chem_species, 1 )) )
1218       ENDIF
1219       IF ( bc_dirichlet_n )  THEN
1220          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1221          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1222          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1223          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1224          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1225          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
1226                                          1:UBOUND( chem_species, 1 )) )
1227       ENDIF
1228       IF ( bc_dirichlet_s )  THEN
1229          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1230          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1231          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1232          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1233          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1234          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
1235                                          1:UBOUND( chem_species, 1 )) )
1236       ENDIF
1237       
1238       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1239       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1240       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1241       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1242       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1243       IF ( air_chemistry )  ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
1244                                       1:UBOUND( chem_species, 1 )) )
1245!
1246!--    For chemical species, create the names of the variables. This is necessary
1247!--    to identify the respective variable and write it onto the correct array
1248!--    in the chem_species datatype.
1249       IF ( air_chemistry )  THEN
1250          ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
1251          ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
1252          ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
1253          ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
1254          ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
1255         
1256          ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
1257          ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
1258          ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
1259          ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
1260          ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
1261!
1262!--       Initialize flags that indicate whether the variable is on file or
1263!--       not. Please note, this is only necessary for chemistry variables.
1264          nest_offl%chem_from_file_l(:) = .FALSE.
1265          nest_offl%chem_from_file_n(:) = .FALSE.
1266          nest_offl%chem_from_file_r(:) = .FALSE.
1267          nest_offl%chem_from_file_s(:) = .FALSE.
1268          nest_offl%chem_from_file_t(:) = .FALSE.
1269         
1270          DO  n = 1, UBOUND( chem_species, 1 )
1271             nest_offl%var_names_chem_l(n) = nest_offl%char_l //               &
1272                                                  TRIM(chem_species(n)%name)
1273             nest_offl%var_names_chem_n(n) = nest_offl%char_n //               &
1274                                                  TRIM(chem_species(n)%name)
1275             nest_offl%var_names_chem_r(n) = nest_offl%char_r //               &
1276                                                  TRIM(chem_species(n)%name)
1277             nest_offl%var_names_chem_s(n) = nest_offl%char_s //               &
1278                                                  TRIM(chem_species(n)%name)
1279             nest_offl%var_names_chem_t(n) = nest_offl%char_t //               &
1280                                                  TRIM(chem_species(n)%name)
1281          ENDDO
1282       ENDIF
1283!
1284!--    Read COSMO data at lateral and top boundaries
1285       CALL netcdf_data_input_offline_nesting
1286!
1287!--    Initialize boundary data. Please note, do not initialize boundaries in
1288!--    case of restart runs. This case the boundaries are already initialized
1289!--    and the boundary data from file would be on the wrong time level. 
1290       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1291          IF ( bc_dirichlet_l )  THEN
1292             u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
1293             v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
1294             w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
1295             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                  &
1296                                      nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
1297             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                  &
1298                                      nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
1299             IF ( air_chemistry )  THEN
1300                DO  n = 1, UBOUND( chem_species, 1 )
1301                   IF( nest_offl%chem_from_file_l(n) )  THEN
1302                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) =             &
1303                                      nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
1304                   ENDIF
1305                ENDDO
1306             ENDIF
1307          ENDIF
1308          IF ( bc_dirichlet_r )  THEN
1309             u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
1310             v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
1311             w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
1312             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =               &
1313                                      nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
1314             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =               &
1315                                      nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
1316             IF ( air_chemistry )  THEN
1317                DO  n = 1, UBOUND( chem_species, 1 )
1318                   IF( nest_offl%chem_from_file_r(n) )  THEN
1319                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) =          &
1320                                      nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
1321                   ENDIF
1322                ENDDO
1323             ENDIF
1324          ENDIF
1325          IF ( bc_dirichlet_s )  THEN
1326             u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
1327             v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
1328             w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
1329             IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                  &
1330                                      nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
1331             IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                  &
1332                                      nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
1333             IF ( air_chemistry )  THEN
1334                DO  n = 1, UBOUND( chem_species, 1 )
1335                   IF( nest_offl%chem_from_file_s(n) )  THEN
1336                      chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) =             &
1337                                      nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
1338                   ENDIF
1339                ENDDO
1340             ENDIF
1341          ENDIF
1342          IF ( bc_dirichlet_n )  THEN
1343             u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
1344             v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
1345             w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
1346             IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =               &
1347                                      nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
1348             IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =               &
1349                                      nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
1350             IF ( air_chemistry )  THEN
1351                DO  n = 1, UBOUND( chem_species, 1 )
1352                   IF( nest_offl%chem_from_file_n(n) )  THEN
1353                      chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =          &
1354                                      nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
1355                   ENDIF
1356                ENDDO
1357             ENDIF
1358          ENDIF
1359!         
1360!--       Initialize geostrophic wind components. Actually this is already done in
1361!--       init_3d_model when initializing_action = 'inifor', however, in speical
1362!--       case of user-defined initialization this will be done here again, in
1363!--       order to have a consistent initialization.
1364          ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
1365          vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
1366!         
1367!--       Set bottom and top boundary condition for geostrophic wind components
1368          ug(nzt+1) = ug(nzt)
1369          vg(nzt+1) = vg(nzt)
1370          ug(nzb)   = ug(nzb+1)
1371          vg(nzb)   = vg(nzb+1)
1372       ENDIF     
1373!
1374!--    After boundary data is initialized, mask topography at the
1375!--    boundaries for the velocity components.
1376       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1377       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1378       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1379!
1380!--    Initial calculation of the boundary layer depth from the prescribed
1381!--    boundary data. This is requiered for initialize the synthetic turbulence
1382!--    generator correctly.
1383       CALL nesting_offl_calc_zi
1384       
1385!
1386!--    After boundary data is initialized, ensure mass conservation. Not
1387!--    necessary in restart runs.
1388       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1389          CALL nesting_offl_mass_conservation
1390       ENDIF
1391
1392    END SUBROUTINE nesting_offl_init
1393   
1394!------------------------------------------------------------------------------!
1395! Description:
1396!------------------------------------------------------------------------------!
1397!> Interpolation function, used to interpolate boundary data in time.
1398!------------------------------------------------------------------------------!
1399    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
1400       
1401       USE kinds
1402
1403       IMPLICIT NONE
1404
1405       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
1406       REAL(wp)            :: var_t1              !< boundary value at t1
1407       REAL(wp)            :: var_t2              !< boundary value at t2
1408       REAL(wp)            :: fac                 !< interpolation factor
1409
1410       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
1411
1412    END FUNCTION interpolate_in_time
1413
1414
1415
1416 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.