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

Last change on this file since 3939 was 3937, checked in by suehring, 5 years ago

Move initialization of STG behind initialization of offline nesting; Bugfix in STG in case of very early restart; calculation of scaling parameters used for parametrization of synthetic turbulence profiles improved; in offline nesting, set boundary value at upper-left and upper-south grid points for u- and v-component, respectively

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