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

Last change on this file since 4012 was 3987, checked in by kanani, 5 years ago

clean up location, debug and error messages

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