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

Last change on this file since 3873 was 3858, checked in by suehring, 5 years ago

Avoid local communication in offline nesting when it is not necessary

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