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

Last change on this file since 3842 was 3737, checked in by suehring, 5 years ago

Enable initialization of chemistry variables via dynamic input file; enable mesoscale offline nesting of chemistry variables if data is in dynamic input file

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