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

Last change on this file since 3981 was 3964, checked in by suehring, 5 years ago

In a nested child domain, distinguish between soil moisture and temperature initialization in case the parent domain is initialized via the dynamic input file; in the offline nesting, add a safety factor for the calculation of bulk Richardson number in order to avoid division by zero which can potentially happen if 3D buildings are located directly at the lateral model boundaries

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