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

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

Synthetic turbulence generator: Revise bias correction of imposed perturbations (correction via volume flow can create instabilities in case the mean volume flow is close to zero); Introduce lower limits in calculation of coefficient matrix, else the calculation may become numerically unstable; Impose perturbations every timestep, even though no new set of perturbations is generated in case dt_stg_call /= dt_3d; Implement a gradual decrease of Reynolds stress and length scales above ABL height (within 1 length scale above ABL depth to 1/10) rather than an abrupt decline; Bugfix in non-nested case: use ABL height for parametrized turbulence; Offline nesting: Rename subroutine for ABL height calculation and make it public; Change top boundary condition for pressure back to Neumann

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