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

Last change on this file since 4168 was 4168, checked in by suehring, 23 months ago

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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