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

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

Enable netcdf parallel input for lateral boundary conditions in dynamic input file

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