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

Last change on this file since 3896 was 3891, checked in by suehring, 5 years ago

Bugfix in initialization of turbulence generator in case of restart runs and offline nesting

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