source: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90 @ 3182

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

  • Property svn:keywords set to Id
File size: 78.2 KB
Line 
1!> @file large_scale_forcing_nudging_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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! * Adjustment to new Inifor version:
23!   - No vertical interpolation/extrapolation of lateral boundary data required
24!     any more (Inifor can treat grid stretching now
25!   - Revise initialization in case of COSMO forcing
26! * Rename variables and subroutines for offline nesting
27!
28! Former revisions:
29! -----------------
30! $Id: large_scale_forcing_nudging_mod.f90 3182 2018-07-27 13:36:03Z suehring $
31! Error messages revised
32!
33! 3045 2018-05-28 07:55:41Z Giersch
34! Error messages revised
35!
36! 3026 2018-05-22 10:30:53Z schwenkel
37! Changed the name specific humidity to mixing ratio, since we are computing
38! mixing ratios.
39!
40! 2970 2018-04-13 15:09:23Z suehring
41! Bugfix in old large-scale forcing mode
42!
43! 2938 2018-03-27 15:52:42Z suehring
44! Further improvements for nesting in larger-scale model
45!
46! 2863 2018-03-08 11:36:25Z suehring
47! Corrected "Former revisions" section
48!
49! 2696 2017-12-14 17:12:51Z kanani
50! Change in file header (GPL part)
51! Forcing with larger-scale models implemented (MS)
52!
53! 2342 2017-08-08 11:00:43Z boeske
54! fixed check if surface forcing data is available until end of simulation
55!
56! 2320 2017-07-21 12:47:43Z suehring
57! initial revision
58!
59! Description:
60! ------------
61!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
62!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
63!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
64!> to given profiles on a relaxation timescale tnudge.
65!> Profiles are read in from NUDGING_DATA.
66!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
67!> @todo: Revise reading of ASCII-files
68!> @todo: Remove unused variables and control flags
69!> @todo: Revise large-scale facing of surface variables
70!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
71!--------------------------------------------------------------------------------!
72 MODULE lsf_nudging_mod
73
74    USE arrays_3d,                                                             &
75        ONLY:  dzw, e, diss, heatflux_input_conversion, pt, pt_init, q,        &
76               q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs,       &
77               waterflux_input_conversion, zu, zw                 
78
79    USE control_parameters,                                                    &
80        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
81               bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
82               constant_heatflux, constant_waterflux,                          &
83               data_output_pr, dt_3d, end_time,                                &
84               humidity, initializing_actions, intermediate_timestep_count,    &
85               ibc_pt_b, ibc_q_b,                                              &
86               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
87               lsf_exception, message_string, nesting_offline, neutral,        &
88               nudging, passive_scalar, pt_surface, ocean, q_surface,          &
89               surface_heatflux, surface_pressure, surface_waterflux,          &
90               topography, use_subsidence_tendencies
91
92    USE grid_variables
93
94    USE indices,                                                               &
95        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
96               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
97
98    USE kinds
99
100    USE netcdf_data_input_mod,                                                 &
101        ONLY:  nest_offl
102
103    USE pegrid
104
105    USE surface_mod,                                                           &
106        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
107
108    USE statistics,                                                            &
109        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
110
111    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
112    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
113
114    REAL(wp) ::  d_area_t
115
116    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
117    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qnudge      !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
118    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tnudge      !< vertical profile of nudging time scale interpolated to vertical grid (nudging) 
119    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_lpt  !< temperature tendency due to large scale advection (large scale forcing)
120    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_q    !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
121    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_lpt  !< temperature tendency due to subsidence/ascent (large scale forcing)
122    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_q    !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
123    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug_vert     !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
124    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  unudge      !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
125    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vnudge      !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
126    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg_vert     !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
127    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wnudge      !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
128    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wsubs_vert  !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
129
130    REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_surf      !< time-dependent surface sensible heat flux (large scale forcing)
131    REAL(wp), DIMENSION(:), ALLOCATABLE ::  timenudge     !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
132    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_surf     !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
133    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_vert     !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
134
135    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tmp_tnudge    !< current nudging time scale
136
137    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_surf        !< time-dependent surface pressure (large scale forcing)
138    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_surf       !< time-dependent surface temperature (large scale forcing)
139    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_surf     !< time-dependent surface latent heat flux (large scale forcing)
140    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_surf        !< time-dependent surface water vapor mixing ratio (large scale forcing)
141
142    SAVE
143    PRIVATE
144!
145!-- Public subroutines
146    PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init,  &
147           lsf_nudging_check_parameters, nudge_init,                           &
148           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
149           lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         & 
150           nudge, nudge_ref
151           
152!
153!-- Public variables
154    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
155           td_sub_q, time_vert
156
157
158    INTERFACE ls_advec
159       MODULE PROCEDURE ls_advec
160       MODULE PROCEDURE ls_advec_ij
161    END INTERFACE ls_advec
162
163    INTERFACE nudge
164       MODULE PROCEDURE nudge
165       MODULE PROCEDURE nudge_ij
166    END INTERFACE nudge
167
168 CONTAINS
169
170
171!------------------------------------------------------------------------------!
172! Description:
173! ------------
174!> In this subroutine a constant mass within the model domain is guaranteed.
175!> Larger-scale models may be based on a compressible equation system, which is
176!> not consistent with PALMs incompressible equation system. In order to avoid
177!> a decrease or increase of mass during the simulation, non-divergent flow
178!> through the lateral and top boundaries is compensated by the vertical wind
179!> component at the top boundary.
180!------------------------------------------------------------------------------!
181    SUBROUTINE lsf_nesting_offline_mass_conservation
182
183       USE control_parameters,                                                 &
184           ONLY:  volume_flow
185
186       IMPLICIT NONE
187
188       INTEGER(iwp) ::  i !< grid index in x-direction
189       INTEGER(iwp) ::  j !< grid index in y-direction
190       INTEGER(iwp) ::  k !< grid index in z-direction
191
192       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
193       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
194
195       volume_flow   = 0.0_wp
196       volume_flow_l = 0.0_wp
197
198       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
199
200       IF ( bc_dirichlet_l )  THEN
201          i = nxl
202          DO  j = nys, nyn
203             DO  k = nzb+1, nzt
204                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
205                                     * MERGE( 1.0_wp, 0.0_wp,                  &
206                                              BTEST( wall_flags_0(k,j,i), 1 ) )
207             ENDDO
208          ENDDO
209       ENDIF
210       IF ( bc_dirichlet_r )  THEN
211          i = nxr+1
212          DO  j = nys, nyn
213             DO  k = nzb+1, nzt
214                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
215                                     * MERGE( 1.0_wp, 0.0_wp,                  &
216                                              BTEST( wall_flags_0(k,j,i), 1 ) )
217             ENDDO
218          ENDDO
219       ENDIF
220       IF ( bc_dirichlet_s )  THEN
221          j = nys
222          DO  i = nxl, nxr
223             DO  k = nzb+1, nzt
224                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
225                                     * MERGE( 1.0_wp, 0.0_wp,                  &
226                                              BTEST( wall_flags_0(k,j,i), 2 ) )
227             ENDDO
228          ENDDO
229       ENDIF
230       IF ( bc_dirichlet_n )  THEN
231          j = nyn+1
232          DO  i = nxl, nxr
233             DO  k = nzb+1, nzt
234                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
235                                     * MERGE( 1.0_wp, 0.0_wp,                  &
236                                              BTEST( wall_flags_0(k,j,i), 2 ) )
237             ENDDO
238          ENDDO
239       ENDIF
240!
241!--    Top boundary
242       k = nzt
243       DO  i = nxl, nxr
244          DO  j = nys, nyn
245             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
246          ENDDO
247       ENDDO
248
249#if defined( __parallel )
250       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
251       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
252                           comm2d, ierr )
253#else
254       volume_flow = volume_flow_l
255#endif
256
257       w_correct = SUM( volume_flow ) * d_area_t
258
259       DO  i = nxl, nxr
260          DO  j = nys, nyn
261             DO  k = nzt, nzt + 1
262                w(k,j,i) = w(k,j,i) + w_correct
263             ENDDO
264          ENDDO
265       ENDDO
266
267    END SUBROUTINE lsf_nesting_offline_mass_conservation
268
269
270!------------------------------------------------------------------------------!
271! Description:
272! ------------
273!> Set the lateral and top boundary conditions in case the PALM domain is
274!> nested offline in a mesoscale model.
275!------------------------------------------------------------------------------!
276    SUBROUTINE lsf_nesting_offline
277
278       USE control_parameters,                                                 &
279           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
280                  bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
281                  rans_tke_e, time_since_reference_point                     
282
283       IMPLICIT NONE
284
285       INTEGER(iwp) ::  i !< running index x-direction
286       INTEGER(iwp) ::  j !< running index y-direction
287       INTEGER(iwp) ::  k !< running index z-direction
288       INTEGER(iwp) ::  t !< running index for time levels
289
290       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
291       REAL(wp) ::  t_ref   !< time past since last reference step
292     
293!
294!--    Calculate time interval of forcing data       
295       ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
296                            nest_offl%time(nest_offl%tind) )
297!
298!--    Calculate reziproke time past since last reference step. Please note,
299!--    the time coordinate is still not updated, so that the actual time need
300!--    to be incremented by dt_3d. Moreover, note that the simulation time
301!--    passed since simulation start is time_since_reference_point, not
302!--    simulated_time!
303       t_ref = time_since_reference_point + dt_3d -                            &
304                                            nest_offl%time(nest_offl%tind)
305                                           
306       IF ( bc_dirichlet_l )  THEN
307
308          DO  j = nys, nyn
309             DO  k = nzb+1, nzt
310                u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
311                       ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
312                         MERGE( 1.0_wp, 0.0_wp,                                &
313                                BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
314             ENDDO
315          ENDDO
316
317          DO  j = nys, nyn
318             DO  k = nzb+1, nzt-1
319                w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
320                       ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
321                         MERGE( 1.0_wp, 0.0_wp,                                &
322                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
323             ENDDO
324          ENDDO
325
326          DO  j = nysv, nyn
327             DO  k = nzb+1, nzt
328                v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
329                       ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
330                         MERGE( 1.0_wp, 0.0_wp,                                &
331                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
332             ENDDO
333          ENDDO
334
335          IF ( .NOT. neutral )  THEN
336             DO  j = nys, nyn
337                DO  k = nzb+1, nzt
338                   pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
339                                                                   t_ref   *   &
340                       ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
341 
342                ENDDO
343             ENDDO
344          ENDIF
345
346          IF ( humidity )  THEN
347             DO  j = nys, nyn
348                DO  k = nzb+1, nzt
349                   q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
350                                                                 t_ref     *   &
351                       ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
352 
353                ENDDO
354             ENDDO
355          ENDIF
356
357       ENDIF
358
359       IF ( bc_dirichlet_r )  THEN
360
361          DO  j = nys, nyn
362             DO  k = nzb+1, nzt
363                u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
364                      ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
365                         MERGE( 1.0_wp, 0.0_wp,                                 &
366                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
367
368             ENDDO
369          ENDDO
370          DO  j = nys, nyn
371             DO  k = nzb+1, nzt-1
372                w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
373                      ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
374                         MERGE( 1.0_wp, 0.0_wp,                                 &
375                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
376             ENDDO
377          ENDDO
378
379          DO  j = nysv, nyn
380             DO  k = nzb+1, nzt
381                v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
382                      ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
383                         MERGE( 1.0_wp, 0.0_wp,                                 &
384                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
385             ENDDO
386          ENDDO
387
388          IF ( .NOT. neutral )  THEN
389             DO  j = nys, nyn
390                DO  k = nzb+1, nzt
391                   pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
392                                                                    t_ref   *  &
393                     ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
394 
395                ENDDO
396             ENDDO
397          ENDIF
398
399          IF ( humidity )  THEN
400             DO  j = nys, nyn
401                DO  k = nzb+1, nzt
402                   q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
403                                                                    t_ref   *  &
404                       ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
405 
406                ENDDO
407             ENDDO
408          ENDIF
409
410       ENDIF
411
412       IF ( bc_dirichlet_s )  THEN
413
414          DO  i = nxl, nxr
415             DO  k = nzb+1, nzt
416                v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
417                      ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
418                         MERGE( 1.0_wp, 0.0_wp,                                 &
419                                BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
420             ENDDO
421          ENDDO
422
423          DO  i = nxl, nxr
424             DO  k = nzb+1, nzt-1
425                w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
426                        ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
427                           MERGE( 1.0_wp, 0.0_wp,                                &
428                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
429             ENDDO
430          ENDDO
431
432          DO  i = nxlu, nxr
433             DO  k = nzb+1, nzt
434                u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
435                        ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
436                           MERGE( 1.0_wp, 0.0_wp,                                &
437                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
438             ENDDO
439          ENDDO
440
441          IF ( .NOT. neutral )  THEN
442             DO  i = nxl, nxr
443                DO  k = nzb+1, nzt
444                   pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
445                                                                    t_ref   *  &
446                     ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
447 
448                ENDDO
449             ENDDO
450          ENDIF
451
452          IF ( humidity )  THEN
453             DO  i = nxl, nxr
454                DO  k = nzb+1, nzt
455                   q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
456                                                                    t_ref   *  &
457                       ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
458 
459                ENDDO
460             ENDDO
461          ENDIF
462
463       ENDIF
464
465       IF ( bc_dirichlet_n )  THEN
466
467          DO  i = nxl, nxr
468             DO  k = nzb+1, nzt
469                v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
470                        ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
471                           MERGE( 1.0_wp, 0.0_wp,                                 &
472                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
473             ENDDO
474          ENDDO
475          DO  i = nxl, nxr
476             DO  k = nzb+1, nzt-1
477                w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
478                        ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
479                           MERGE( 1.0_wp, 0.0_wp,                                &
480                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
481             ENDDO
482          ENDDO
483
484          DO  i = nxlu, nxr
485             DO  k = nzb+1, nzt
486                u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
487                        ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
488                           MERGE( 1.0_wp, 0.0_wp,                                &
489                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
490
491             ENDDO
492          ENDDO
493
494          IF ( .NOT. neutral )  THEN
495             DO  i = nxl, nxr
496                DO  k = nzb+1, nzt
497                   pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
498                                                                    t_ref   *  &
499                     ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
500 
501                ENDDO
502             ENDDO
503          ENDIF
504
505          IF ( humidity )  THEN
506             DO  i = nxl, nxr
507                DO  k = nzb+1, nzt
508                   q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
509                                                                    t_ref   *  &
510                       ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
511 
512                ENDDO
513             ENDDO
514          ENDIF
515
516       ENDIF
517!
518!--    Top boundary.
519       DO  i = nxlu, nxr
520          DO  j = nys, nyn
521             u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
522                        ( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) *  &
523                           MERGE( 1.0_wp, 0.0_wp,                              &
524                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
525          ENDDO
526       ENDDO
527
528       DO  i = nxl, nxr
529          DO  j = nysv, nyn
530             v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
531                        ( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) *  &
532                           MERGE( 1.0_wp, 0.0_wp,                              &
533                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
534          ENDDO
535       ENDDO
536
537       DO  i = nxl, nxr
538          DO  j = nys, nyn
539             w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
540                        ( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) *  &
541                           MERGE( 1.0_wp, 0.0_wp,                              &
542                                  BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
543          ENDDO
544       ENDDO
545
546
547       IF ( .NOT. neutral )  THEN
548          DO  i = nxl, nxr
549             DO  j = nys, nyn
550                pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
551                        ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
552             ENDDO
553          ENDDO
554       ENDIF
555
556       IF ( humidity )  THEN
557          DO  i = nxl, nxr
558             DO  j = nys, nyn
559                q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
560                        ( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
561             ENDDO
562          ENDDO
563       ENDIF
564!
565!--    At the edges( left-south, left-north, right-south and right-north) set
566!--    data on ghost points.
567       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_s )  THEN
568          DO  i = 1, nbgp
569             u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
570             w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
571             IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
572             IF ( humidity      )  q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
573          ENDDO
574          DO  i = 1, nbgp+1
575             v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
576          ENDDO
577       ENDIF
578       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_n )  THEN
579          DO  i = 1, nbgp
580             u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
581             v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
582             w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
583             IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
584             IF ( humidity      )  q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
585          ENDDO
586       ENDIF
587       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_s )  THEN
588          DO  i = 1, nbgp
589             u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
590             w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
591             IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
592             IF ( humidity      )  q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
593          ENDDO
594          DO  i = 1, nbgp+1
595             v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
596          ENDDO
597       ENDIF
598       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_n )  THEN
599          DO  i = 1, nbgp
600             u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
601             v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
602             w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
603             IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
604             IF ( humidity      )  q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
605          ENDDO
606       ENDIF
607!
608!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
609!--    passive scalar, dissipation, and chemical species if required
610       IF ( rans_mode  .AND.  rans_tke_e )  THEN
611          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
612          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
613          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
614          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
615       ENDIF
616       IF ( .NOT. constant_diffusion )  THEN
617          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
618          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
619          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
620          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
621          e(nzt+1,:,:) = e(nzt,:,:)
622       ENDIF
623       IF ( passive_scalar )  THEN
624          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
625          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
626          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
627          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
628       ENDIF
629
630
631
632       CALL exchange_horiz( u, nbgp )
633       CALL exchange_horiz( v, nbgp )
634       CALL exchange_horiz( w, nbgp )
635       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
636       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
637
638!
639!--    Set surface pressure. Please note, time-dependent surface
640!--    pressure would require changes in anelastic approximation and
641!--    treatment of fluxes.
642!--    For the moment, comment this out!
643!      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
644!                                                      ddt_lsf * t_ref *       &
645!                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
646!                                    - nest_offl%surface_pressure(nest_offl%tind) )
647
648    END SUBROUTINE lsf_nesting_offline
649
650!------------------------------------------------------------------------------!
651! Description:
652! ------------
653!> @todo Missing subroutine description.
654!------------------------------------------------------------------------------!
655    SUBROUTINE lsf_nudging_check_parameters
656
657       IMPLICIT NONE
658!
659!--    Check nudging and large scale forcing from external file
660       IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
661          message_string = 'Nudging requires large_scale_forcing = .T.. &'//   &
662                        'Surface fluxes and geostrophic wind should be &'//    &
663                        'prescribed in file LSF_DATA'
664          CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
665       ENDIF
666
667       IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.              &
668                                          bc_ns /= 'cyclic' ) )  THEN
669          message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
670                        'the usage of large scale forcing from external file.'
671          CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
672       ENDIF
673
674       IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
675          message_string = 'The usage of large scale forcing from external &'//&
676                        'file LSF_DATA requires humidity = .T..'
677          CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
678       ENDIF
679
680       IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
681          message_string = 'The usage of large scale forcing from external &'// &
682                        'file LSF_DATA is not implemented for passive scalars'
683          CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
684       ENDIF
685
686       IF ( large_scale_forcing  .AND.  topography /= 'flat'                   &
687                              .AND.  .NOT.  lsf_exception )  THEN
688          message_string = 'The usage of large scale forcing from external &'//&
689                        'file LSF_DATA is not implemented for non-flat topography'
690          CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
691       ENDIF
692
693       IF ( large_scale_forcing  .AND.  ocean )  THEN
694          message_string = 'The usage of large scale forcing from external &'//&
695                        'file LSF_DATA is not implemented for ocean runs'
696          CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
697       ENDIF
698
699    END SUBROUTINE lsf_nudging_check_parameters
700
701!------------------------------------------------------------------------------!
702! Description:
703! ------------
704!> Check data output of profiles for land surface model
705!------------------------------------------------------------------------------!
706    SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit,    &
707                                                 dopr_unit )
708 
709       USE profil_parameter
710
711       IMPLICIT NONE
712   
713       CHARACTER (LEN=*) ::  unit      !<
714       CHARACTER (LEN=*) ::  variable  !<
715       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
716 
717       INTEGER(iwp) ::  user_pr_index !<   
718       INTEGER(iwp) ::  var_count     !<
719
720       SELECT CASE ( TRIM( variable ) )
721       
722
723          CASE ( 'td_lsa_lpt' )
724             IF (  .NOT.  large_scale_forcing )  THEN
725                message_string = 'data_output_pr = ' //                        &
726                                 TRIM( data_output_pr(var_count) ) //          &
727                                 ' is not implemented for ' //                 &
728                                 'large_scale_forcing = .FALSE.'
729                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
730                               1, 2, 0, 6, 0 )
731             ELSE
732                dopr_index(var_count) = 81
733                dopr_unit             = 'K/s'
734                unit                  = 'K/s'
735                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
736             ENDIF
737
738          CASE ( 'td_lsa_q' )
739             IF (  .NOT.  large_scale_forcing )  THEN
740                message_string = 'data_output_pr = ' //                        &
741                                 TRIM( data_output_pr(var_count) ) //          &
742                                 ' is not implemented for ' //                 &
743                                 'large_scale_forcing = .FALSE.'
744                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
745                               1, 2, 0, 6, 0 )
746             ELSE
747                dopr_index(var_count) = 82
748                dopr_unit             = 'kg/kgs'
749                unit                  = 'kg/kgs'
750                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
751             ENDIF
752          CASE ( 'td_sub_lpt' )
753             IF (  .NOT.  large_scale_forcing )  THEN
754                message_string = 'data_output_pr = ' //                        &
755                                 TRIM( data_output_pr(var_count) ) //          &
756                                 ' is not implemented for ' //                 &
757                                 'large_scale_forcing = .FALSE.'
758                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
759                               1, 2, 0, 6, 0 )
760             ELSE
761                dopr_index(var_count) = 83
762                dopr_unit             = 'K/s'
763                unit                  = 'K/s'
764                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
765             ENDIF
766
767          CASE ( 'td_sub_q' )
768             IF (  .NOT.  large_scale_forcing )  THEN
769                message_string = 'data_output_pr = ' //                        &
770                                 TRIM( data_output_pr(var_count) ) //          &
771                                 ' is not implemented for ' //                 &
772                                 'large_scale_forcing = .FALSE.'
773                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
774                               1, 2, 0, 6, 0 )
775             ELSE
776                dopr_index(var_count) = 84
777                dopr_unit             = 'kg/kgs'
778                unit                  = 'kg/kgs'
779                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
780             ENDIF
781
782          CASE ( 'td_nud_lpt' )
783             IF (  .NOT.  nudging )  THEN
784                message_string = 'data_output_pr = ' //                        &
785                                 TRIM( data_output_pr(var_count) ) //          &
786                                 ' is not implemented for ' //                 &
787                                 'nudging = .FALSE.'
788                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
789                               1, 2, 0, 6, 0 )
790             ELSE
791                dopr_index(var_count) = 85
792                dopr_unit             = 'K/s'
793                unit                  = 'K/s'
794                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
795             ENDIF
796
797          CASE ( 'td_nud_q' )
798             IF (  .NOT.  nudging )  THEN
799                message_string = 'data_output_pr = ' //                        &
800                                 TRIM( data_output_pr(var_count) ) //          &
801                                 ' is not implemented for ' //                 &
802                                 'nudging = .FALSE.'
803                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
804                               1, 2, 0, 6, 0 )
805             ELSE
806                dopr_index(var_count) = 86
807                dopr_unit             = 'kg/kgs'
808                unit                  = 'kg/kgs'
809                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
810             ENDIF
811
812          CASE ( 'td_nud_u' )
813             IF (  .NOT.  nudging )  THEN
814                message_string = 'data_output_pr = ' //                        &
815                                 TRIM( data_output_pr(var_count) ) //          &
816                                 ' is not implemented for ' //                 &
817                                 'nudging = .FALSE.'
818                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
819                               1, 2, 0, 6, 0 )
820             ELSE
821                dopr_index(var_count) = 87
822                dopr_unit             = 'm/s2'
823                unit                  = 'm/s2'
824                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
825             ENDIF
826
827          CASE ( 'td_nud_v' )
828             IF (  .NOT.  nudging )  THEN
829                message_string = 'data_output_pr = ' //                        &
830                                 TRIM( data_output_pr(var_count) ) //          &
831                                 ' is not implemented for ' //                 &
832                                 'nudging = .FALSE.'
833                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
834                               1, 2, 0, 6, 0 )
835             ELSE
836                dopr_index(var_count) = 88
837                dopr_unit             = 'm/s2'
838                unit                  = 'm/s2'
839                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
840             ENDIF
841
842
843          CASE DEFAULT
844             unit = 'illegal'
845   
846       END SELECT
847
848    END SUBROUTINE lsf_nudging_check_data_output_pr
849
850!------------------------------------------------------------------------------!
851! Description:
852! ------------
853!> @todo Missing subroutine description.
854!------------------------------------------------------------------------------!
855    SUBROUTINE lsf_nudging_header ( io )
856
857       IMPLICIT NONE
858
859       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
860
861       WRITE ( io, 1 )
862       IF ( large_scale_forcing )  THEN
863          WRITE ( io, 3 )
864          WRITE ( io, 4 )
865
866          IF ( large_scale_subsidence )  THEN
867             IF ( .NOT. use_subsidence_tendencies )  THEN
868                WRITE ( io, 5 )
869             ELSE
870                WRITE ( io, 6 )
871             ENDIF
872          ENDIF
873
874          IF ( bc_pt_b == 'dirichlet' )  THEN
875             WRITE ( io, 12 )
876          ELSEIF ( bc_pt_b == 'neumann' )  THEN
877             WRITE ( io, 13 )
878          ENDIF
879
880          IF ( bc_q_b == 'dirichlet' )  THEN
881             WRITE ( io, 14 )
882          ELSEIF ( bc_q_b == 'neumann' )  THEN
883             WRITE ( io, 15 )
884          ENDIF
885
886          WRITE ( io, 7 )
887          IF ( nudging )  THEN
888             WRITE ( io, 10 )
889          ENDIF
890       ELSE
891          WRITE ( io, 2 )
892          WRITE ( io, 11 )
893       ENDIF
894       IF ( large_scale_subsidence )  THEN
895          WRITE ( io, 8 )
896          WRITE ( io, 9 )
897       ENDIF
898
899
9001 FORMAT (//' Large scale forcing and nudging:'/ &
901              ' -------------------------------'/)
9022 FORMAT (' --> No large scale forcing from external is used (default) ')
9033 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
9044 FORMAT ('     - large scale advection tendencies ')
9055 FORMAT ('     - large scale subsidence velocity w_subs ')
9066 FORMAT ('     - large scale subsidence tendencies ')
9077 FORMAT ('     - and geostrophic wind components ug and vg')
9088 FORMAT (' --> Large-scale vertical motion is used in the ', &
909                  'prognostic equation(s) for')
9109 FORMAT ('     the scalar(s) only')
91110 FORMAT (' --> Nudging is used')
91211 FORMAT (' --> No nudging is used (default) ')
91312 FORMAT ('     - prescribed surface values for temperature')
91413 FORMAT ('     - prescribed surface fluxes for temperature')
91514 FORMAT ('     - prescribed surface values for humidity')
91615 FORMAT ('     - prescribed surface fluxes for humidity')
917
918    END SUBROUTINE lsf_nudging_header
919
920!------------------------------------------------------------------------------!
921! Description:
922! ------------
923!> @todo Missing subroutine description.
924!------------------------------------------------------------------------------!
925    SUBROUTINE lsf_init
926
927       USE control_parameters,                                                 &
928           ONLY:  bc_lr_cyc, bc_ns_cyc
929
930       USE netcdf_data_input_mod,                                              &
931           ONLY:  netcdf_data_input_lsf 
932
933       IMPLICIT NONE
934
935       CHARACTER(100) ::  chmess      !<
936       CHARACTER(1)   ::  hash        !<
937
938       INTEGER(iwp) ::  ierrn         !<
939       INTEGER(iwp) ::  finput = 90   !<
940       INTEGER(iwp) ::  k             !<
941       INTEGER(iwp) ::  nt            !<
942       INTEGER(iwp) ::  t             !< running index for time levels
943
944       REAL(wp) ::  fac               !<
945       REAL(wp) ::  highheight        !<
946       REAL(wp) ::  highug_vert       !<
947       REAL(wp) ::  highvg_vert       !<
948       REAL(wp) ::  highwsubs_vert    !<
949       REAL(wp) ::  lowheight         !<
950       REAL(wp) ::  lowug_vert        !<
951       REAL(wp) ::  lowvg_vert        !<
952       REAL(wp) ::  lowwsubs_vert     !<
953       REAL(wp) ::  high_td_lsa_lpt   !<
954       REAL(wp) ::  low_td_lsa_lpt    !<
955       REAL(wp) ::  high_td_lsa_q     !<
956       REAL(wp) ::  low_td_lsa_q      !<
957       REAL(wp) ::  high_td_sub_lpt   !<
958       REAL(wp) ::  low_td_sub_lpt    !<
959       REAL(wp) ::  high_td_sub_q     !<
960       REAL(wp) ::  low_td_sub_q      !<
961       REAL(wp) ::  r_dummy           !<
962
963       IF ( nesting_offline )  THEN
964!
965!--       Allocate arrays for geostrophic wind components. Arrays will
966!--       incorporate 2 time levels in order to interpolate in between. Please
967!--       note, forcing using geostrophic wind components is only required in
968!--       case of cyclic boundary conditions.
969          IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
970             ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
971             ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
972          ENDIF
973!
974!--       Allocate arrays for reading boundary values. Arrays will incorporate 2
975!--       time levels in order to interpolate in between.
976          IF ( bc_dirichlet_l )  THEN
977             ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
978             ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
979             ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
980             IF ( humidity )      ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
981             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
982          ENDIF
983          IF ( bc_dirichlet_r )  THEN
984             ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
985             ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
986             ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
987             IF ( humidity )      ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
988             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
989          ENDIF
990          IF ( bc_dirichlet_n )  THEN
991             ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
992             ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
993             ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
994             IF ( humidity )      ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
995             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
996          ENDIF
997          IF ( bc_dirichlet_s )  THEN
998             ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
999             ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1000             ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1001             IF ( humidity )      ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1002             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1003          ENDIF
1004         
1005          ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1006          ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1007          ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1008          IF ( humidity )      ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1009          IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1010
1011!
1012!--       Read COSMO data at lateral and top boundaries
1013          CALL netcdf_data_input_lsf
1014!
1015!--       Write COSMO data at lateral and top boundaries
1016          CALL lsf_nesting_offline
1017!
1018!--       After 3D data is initialized, ensure mass conservation
1019          CALL lsf_nesting_offline_mass_conservation
1020!
1021!--       Initialize surface pressure. Please note, time-dependent surface
1022!--       pressure would require changes in anelastic approximation and
1023!--       treatment of fluxes.
1024!--       For the moment, comment this out!
1025!         surface_pressure = nest_offl%surface_pressure(0)
1026
1027       ELSE
1028
1029          ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
1030                    qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
1031                    td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
1032                    td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
1033                    time_vert(0:nlsf), time_surf(0:nlsf),                      &
1034                    ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
1035                    wsubs_vert(nzb:nzt+1,0:nlsf) )
1036
1037          p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
1038          shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
1039          td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
1040          time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
1041          wsubs_vert = 0.0_wp
1042
1043!
1044!--       Array for storing large scale forcing and nudging tendencies at each
1045!--       timestep for data output
1046          ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
1047          sums_ls_l = 0.0_wp
1048
1049          ngp_sums_ls = (nz+2)*6
1050
1051          OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
1052                 FORM='FORMATTED', IOSTAT=ierrn )
1053
1054          IF ( ierrn /= 0 )  THEN
1055             message_string = 'file LSF_DATA does not exist'
1056             CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
1057          ENDIF
1058
1059          ierrn = 0
1060!
1061!--       First three lines of LSF_DATA contain header
1062          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1063          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1064          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1065
1066          IF ( ierrn /= 0 )  THEN
1067             message_string = 'errors in file LSF_DATA'
1068             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1069          ENDIF
1070
1071!
1072!--       Surface values are read in
1073          nt     = 0
1074          ierrn = 0
1075
1076          DO WHILE ( time_surf(nt) < end_time )
1077             nt = nt + 1
1078             READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
1079                                                qsws_surf(nt), pt_surf(nt),    &
1080                                                q_surf(nt), p_surf(nt)
1081
1082             IF ( ierrn /= 0 )  THEN
1083               WRITE ( message_string, * ) 'No time dependent surface ' //     &
1084                                 'variables in & LSF_DATA for end of run found'
1085
1086                CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
1087             ENDIF
1088          ENDDO
1089
1090          IF ( time_surf(1) > end_time )  THEN
1091             WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
1092                                         '&LSF_DATA set in after end of ' ,        &
1093                                         'simulation - lsf_surf is set to FALSE'
1094             CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
1095             lsf_surf = .FALSE.
1096          ENDIF
1097
1098!
1099!--       Go to the end of the list with surface variables
1100          DO WHILE ( ierrn == 0 )
1101             READ ( finput, *, IOSTAT = ierrn ) r_dummy
1102          ENDDO
1103
1104!
1105!--       Profiles of ug, vg and w_subs are read in (large scale forcing)
1106
1107          nt = 0
1108          DO WHILE ( time_vert(nt) < end_time )
1109             nt = nt + 1
1110             hash = "#"
1111             ierrn = 1 ! not zero
1112!
1113!--          Search for the next line consisting of "# time",
1114!--          from there onwards the profiles will be read
1115             DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
1116                READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
1117                IF ( ierrn < 0 )  THEN
1118                   WRITE( message_string, * ) 'No time dependent vertical profiles',&
1119                                    ' in & LSF_DATA for end of run found'
1120                   CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
1121                ENDIF
1122             ENDDO
1123
1124             IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
1125
1126             READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
1127                                              lowwsubs_vert, low_td_lsa_lpt,    &
1128                                              low_td_lsa_q, low_td_sub_lpt,     &
1129                                              low_td_sub_q
1130             IF ( ierrn /= 0 )  THEN
1131                message_string = 'errors in file LSF_DATA'
1132                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1133             ENDIF
1134
1135             READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
1136                                              highvg_vert, highwsubs_vert,     &
1137                                              high_td_lsa_lpt, high_td_lsa_q,  &
1138                                              high_td_sub_lpt, high_td_sub_q
1139         
1140             IF ( ierrn /= 0 )  THEN
1141                message_string = 'errors in file LSF_DATA'
1142                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1143             ENDIF
1144
1145
1146             DO  k = nzb, nzt+1
1147                IF ( highheight < zu(k) )  THEN
1148                   lowheight      = highheight
1149                   lowug_vert     = highug_vert
1150                   lowvg_vert     = highvg_vert
1151                   lowwsubs_vert  = highwsubs_vert
1152                   low_td_lsa_lpt = high_td_lsa_lpt
1153                   low_td_lsa_q   = high_td_lsa_q
1154                   low_td_sub_lpt = high_td_sub_lpt
1155                   low_td_sub_q   = high_td_sub_q
1156
1157                   ierrn = 0
1158                   READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
1159                                                    highvg_vert, highwsubs_vert,&
1160                                                    high_td_lsa_lpt,            &
1161                                                    high_td_lsa_q,              &
1162                                                    high_td_sub_lpt, high_td_sub_q
1163
1164                   IF ( ierrn /= 0 )  THEN
1165                      WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',  &
1166                           'is higher than the maximum height in LSF_DATA ',   &
1167                           'which is ', lowheight, 'm. Interpolation on PALM ',&
1168                           'grid is not possible.'
1169                      CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
1170                   ENDIF
1171
1172                ENDIF
1173
1174!
1175!--             Interpolation of prescribed profiles in space
1176                fac = (highheight-zu(k))/(highheight - lowheight)
1177
1178                ug_vert(k,nt)    = fac * lowug_vert                            &
1179                                   + ( 1.0_wp - fac ) * highug_vert
1180                vg_vert(k,nt)    = fac * lowvg_vert                            &
1181                                   + ( 1.0_wp - fac ) * highvg_vert
1182                wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
1183                                   + ( 1.0_wp - fac ) * highwsubs_vert
1184
1185                td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
1186                                   + ( 1.0_wp - fac ) * high_td_lsa_lpt
1187                td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
1188                                   + ( 1.0_wp - fac ) * high_td_lsa_q
1189                td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
1190                                   + ( 1.0_wp - fac ) * high_td_sub_lpt
1191                td_sub_q(k,nt)   = fac * low_td_sub_q                          &
1192                                   + ( 1.0_wp - fac ) * high_td_sub_q
1193
1194             ENDDO
1195
1196          ENDDO
1197
1198!
1199!--       Large scale vertical velocity has to be zero at the surface
1200          wsubs_vert(nzb,:) = 0.0_wp
1201   
1202          IF ( time_vert(1) > end_time )  THEN
1203             WRITE ( message_string, * ) 'Time dependent large scale profile ',&
1204                                'forcing from&LSF_DATA sets in after end of ' ,&
1205                                'simulation - lsf_vert is set to FALSE'
1206             CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
1207             lsf_vert = .FALSE.
1208          ENDIF
1209
1210          CLOSE( finput )
1211
1212       ENDIF
1213
1214    END SUBROUTINE lsf_init
1215
1216!------------------------------------------------------------------------------!
1217! Description:
1218! ------------
1219!> @todo Missing subroutine description.
1220!------------------------------------------------------------------------------!
1221    SUBROUTINE ls_forcing_surf ( time )
1222
1223       IMPLICIT NONE
1224
1225       INTEGER(iwp) ::  nt                     !<
1226
1227       REAL(wp)             :: dum_surf_flux  !<
1228       REAL(wp)             :: fac            !<
1229       REAL(wp), INTENT(in) :: time           !<
1230
1231!
1232!--    Interpolation in time of LSF_DATA at the surface
1233       nt = 1
1234       DO WHILE ( time > time_surf(nt) )
1235          nt = nt + 1
1236       ENDDO
1237       IF ( time /= time_surf(nt) )  THEN
1238         nt = nt - 1
1239       ENDIF
1240
1241       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
1242
1243       IF ( ibc_pt_b == 0 )  THEN
1244!
1245!--       In case of Dirichlet boundary condition shf must not
1246!--       be set - it is calculated via MOST in prandtl_fluxes
1247          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1248
1249       ELSEIF ( ibc_pt_b == 1 )  THEN
1250!
1251!--       In case of Neumann boundary condition pt_surface is needed for
1252!--       calculation of reference density
1253          dum_surf_flux = ( shf_surf(nt) + fac *                               &
1254                            ( shf_surf(nt+1) - shf_surf(nt) )                  &
1255                          ) * heatflux_input_conversion(nzb)
1256!
1257!--       Save surface sensible heat flux on default, natural and urban surface
1258!--       type, if required
1259          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%shf(:) = dum_surf_flux
1260          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%shf(:)    = dum_surf_flux
1261          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%shf(:)    = dum_surf_flux
1262
1263          pt_surface    = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1264
1265       ENDIF
1266
1267       IF ( ibc_q_b == 0 )  THEN
1268!
1269!--       In case of Dirichlet boundary condition qsws must not
1270!--       be set - it is calculated via MOST in prandtl_fluxes
1271          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
1272
1273       ELSEIF ( ibc_q_b == 1 )  THEN
1274          dum_surf_flux = ( qsws_surf(nt) + fac *                              &
1275                             ( qsws_surf(nt+1) - qsws_surf(nt) )               &
1276                             ) * waterflux_input_conversion(nzb)
1277!
1278!--       Save surface sensible heat flux on default, natural and urban surface
1279!--       type, if required
1280          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%qsws(:) = dum_surf_flux
1281          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%qsws(:)    = dum_surf_flux
1282          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%qsws(:)    = dum_surf_flux
1283
1284       ENDIF
1285!
1286!--    Surface heat- and waterflux will be written later onto surface elements
1287       IF ( .NOT.  neutral  .AND.  constant_heatflux  .AND.                    &
1288            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1289             surface_heatflux = shf_surf(1)
1290       ENDIF
1291       IF ( humidity  .AND.  constant_waterflux  .AND.                         &
1292            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1293             surface_waterflux = qsws_surf(1)
1294       ENDIF
1295
1296       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
1297
1298    END SUBROUTINE ls_forcing_surf
1299
1300
1301
1302
1303!------------------------------------------------------------------------------!
1304! Description:
1305! ------------
1306!> @todo Missing subroutine description.
1307!------------------------------------------------------------------------------!
1308    SUBROUTINE ls_forcing_vert ( time )
1309
1310
1311       IMPLICIT NONE
1312
1313       INTEGER(iwp) ::  nt                     !<
1314
1315       REAL(wp)             ::  fac           !<
1316       REAL(wp), INTENT(in) ::  time          !<
1317
1318!
1319!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
1320       nt = 1
1321       DO WHILE ( time > time_vert(nt) )
1322          nt = nt + 1
1323       ENDDO
1324       IF ( time /= time_vert(nt) )  THEN
1325         nt = nt - 1
1326       ENDIF
1327
1328       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1329
1330       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
1331       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
1332
1333       IF ( large_scale_subsidence )  THEN
1334          w_subs = wsubs_vert(:,nt)                                            &
1335                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
1336       ENDIF
1337
1338    END SUBROUTINE ls_forcing_vert
1339
1340
1341!------------------------------------------------------------------------------!
1342! Description:
1343! ------------
1344!> Call for all grid points
1345!------------------------------------------------------------------------------!
1346    SUBROUTINE ls_advec ( time, prog_var )
1347     
1348
1349       IMPLICIT NONE
1350
1351       CHARACTER (LEN=*) ::  prog_var   !<
1352
1353       REAL(wp), INTENT(in)  :: time    !<
1354       REAL(wp) :: fac                  !< 
1355
1356       INTEGER(iwp) ::  i               !<
1357       INTEGER(iwp) ::  j               !<
1358       INTEGER(iwp) ::  k               !<
1359       INTEGER(iwp) ::  nt               !<
1360
1361!
1362!--    Interpolation in time of LSF_DATA
1363       nt = 1
1364       DO WHILE ( time > time_vert(nt) )
1365          nt = nt + 1
1366       ENDDO
1367       IF ( time /= time_vert(nt) )  THEN
1368         nt = nt - 1
1369       ENDIF
1370
1371       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1372
1373!
1374!--    Add horizontal large scale advection tendencies of pt and q
1375       SELECT CASE ( prog_var )
1376
1377          CASE ( 'pt' )
1378
1379             DO  i = nxl, nxr
1380                DO  j = nys, nyn
1381                   DO  k = nzb+1, nzt
1382                      tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac *     &
1383                                    ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
1384                                        MERGE( 1.0_wp, 0.0_wp,                 &
1385                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1386                   ENDDO
1387                ENDDO
1388             ENDDO
1389
1390          CASE ( 'q' )
1391
1392             DO  i = nxl, nxr
1393                DO  j = nys, nyn
1394                   DO  k = nzb+1, nzt
1395                      tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac *       &
1396                                    ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *    &
1397                                        MERGE( 1.0_wp, 0.0_wp,                 &
1398                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1399                   ENDDO
1400                ENDDO
1401             ENDDO
1402
1403       END SELECT
1404
1405!
1406!--    Subsidence of pt and q with prescribed subsidence tendencies
1407       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1408
1409          SELECT CASE ( prog_var )
1410
1411             CASE ( 'pt' )
1412
1413                DO  i = nxl, nxr
1414                   DO  j = nys, nyn
1415                      DO  k = nzb+1, nzt
1416                         tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *  &
1417                                     ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
1418                                        MERGE( 1.0_wp, 0.0_wp,                 &
1419                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1420                      ENDDO
1421                   ENDDO
1422                ENDDO
1423 
1424             CASE ( 'q' )
1425
1426                DO  i = nxl, nxr
1427                   DO  j = nys, nyn
1428                      DO  k = nzb+1, nzt
1429                         tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *    &
1430                                       ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
1431                                        MERGE( 1.0_wp, 0.0_wp,                 &
1432                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1433                      ENDDO
1434                   ENDDO
1435                ENDDO
1436
1437          END SELECT
1438
1439       ENDIF
1440
1441    END SUBROUTINE ls_advec
1442
1443
1444!------------------------------------------------------------------------------!
1445! Description:
1446! ------------
1447!> Call for grid point i,j
1448!------------------------------------------------------------------------------!
1449    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
1450
1451       IMPLICIT NONE
1452
1453       CHARACTER (LEN=*) ::  prog_var   !<
1454
1455       REAL(wp), INTENT(in)  :: time    !<
1456       REAL(wp) :: fac                  !<
1457
1458       INTEGER(iwp) ::  i               !<
1459       INTEGER(iwp) ::  j               !<
1460       INTEGER(iwp) ::  k               !<
1461       INTEGER(iwp) ::  nt               !<
1462
1463!
1464!--    Interpolation in time of LSF_DATA
1465       nt = 1
1466       DO WHILE ( time > time_vert(nt) )
1467          nt = nt + 1
1468       ENDDO
1469       IF ( time /= time_vert(nt) )  THEN
1470         nt = nt - 1
1471       ENDIF
1472
1473       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1474
1475!
1476!--    Add horizontal large scale advection tendencies of pt and q
1477       SELECT CASE ( prog_var )
1478
1479          CASE ( 'pt' )
1480
1481             DO  k = nzb+1, nzt
1482                tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt)                   &
1483                             + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
1484                                        MERGE( 1.0_wp, 0.0_wp,                 &
1485                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1486             ENDDO
1487
1488          CASE ( 'q' )
1489
1490             DO  k = nzb+1, nzt
1491                tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt)                     &
1492                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *  &
1493                                        MERGE( 1.0_wp, 0.0_wp,                 &
1494                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1495             ENDDO
1496
1497       END SELECT
1498
1499!
1500!--    Subsidence of pt and q with prescribed profiles
1501       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1502
1503          SELECT CASE ( prog_var )
1504
1505             CASE ( 'pt' )
1506
1507                DO  k = nzb+1, nzt
1508                   tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *        &
1509                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) *   &
1510                                        MERGE( 1.0_wp, 0.0_wp,                 &
1511                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1512                ENDDO
1513 
1514             CASE ( 'q' )
1515
1516                DO  k = nzb+1, nzt
1517                   tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *          &
1518                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) *       &
1519                                        MERGE( 1.0_wp, 0.0_wp,                 &
1520                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1521                ENDDO
1522
1523          END SELECT
1524
1525       ENDIF
1526
1527    END SUBROUTINE ls_advec_ij
1528
1529
1530!------------------------------------------------------------------------------!
1531! Description:
1532! ------------
1533!> @todo Missing subroutine description.
1534!------------------------------------------------------------------------------!
1535    SUBROUTINE nudge_init
1536
1537       IMPLICIT NONE
1538
1539
1540       INTEGER(iwp) ::  finput = 90  !<
1541       INTEGER(iwp) ::  ierrn        !<
1542       INTEGER(iwp) ::  k            !<
1543       INTEGER(iwp) ::  nt            !<
1544
1545       CHARACTER(1) ::  hash     !<
1546
1547       REAL(wp) ::  highheight   !<
1548       REAL(wp) ::  highqnudge   !<
1549       REAL(wp) ::  highptnudge  !<
1550       REAL(wp) ::  highunudge   !<
1551       REAL(wp) ::  highvnudge   !<
1552       REAL(wp) ::  highwnudge   !<
1553       REAL(wp) ::  hightnudge   !<
1554
1555       REAL(wp) ::  lowheight    !<
1556       REAL(wp) ::  lowqnudge    !<
1557       REAL(wp) ::  lowptnudge   !<
1558       REAL(wp) ::  lowunudge    !<
1559       REAL(wp) ::  lowvnudge    !<
1560       REAL(wp) ::  lowwnudge    !<
1561       REAL(wp) ::  lowtnudge    !<
1562
1563       REAL(wp) ::  fac          !<
1564
1565       ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
1566                 tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge),  &
1567                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
1568
1569       ALLOCATE( tmp_tnudge(nzb:nzt) )
1570
1571       ALLOCATE( timenudge(0:ntnudge) )
1572
1573       ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
1574       vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
1575!
1576!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
1577       tmp_tnudge = 21600.0_wp
1578
1579       nt = 0
1580       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
1581              FORM='FORMATTED', IOSTAT=ierrn )
1582
1583       IF ( ierrn /= 0 )  THEN
1584          message_string = 'file NUDGING_DATA does not exist'
1585          CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
1586       ENDIF
1587
1588       ierrn = 0
1589
1590 rloop:DO
1591          nt = nt + 1
1592          hash = "#"
1593          ierrn = 1 ! not zero
1594!
1595!--       Search for the next line consisting of "# time",
1596!--       from there onwards the profiles will be read
1597          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
1598         
1599            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
1600            IF ( ierrn < 0 )  EXIT rloop
1601
1602          ENDDO
1603
1604          ierrn = 0
1605          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge,   &
1606                                           lowvnudge, lowwnudge , lowptnudge, &
1607                                           lowqnudge
1608
1609          IF ( ierrn /= 0 )  THEN
1610             message_string = 'errors in file NUDGING_DATA'
1611             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1612          ENDIF
1613
1614          ierrn = 0
1615          READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge,   &
1616                                           highvnudge, highwnudge , highptnudge, &
1617                                           highqnudge
1618
1619          IF ( ierrn /= 0 )  THEN
1620             message_string = 'errors in file NUDGING_DATA'
1621             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1622          ENDIF
1623
1624          DO  k = nzb, nzt+1
1625             DO WHILE ( highheight < zu(k) )
1626                lowheight  = highheight
1627                lowtnudge  = hightnudge
1628                lowunudge  = highunudge
1629                lowvnudge  = highvnudge
1630                lowwnudge  = highwnudge
1631                lowptnudge = highptnudge
1632                lowqnudge  = highqnudge
1633 
1634                ierrn = 0
1635                READ ( finput, *, IOSTAT=ierrn )  highheight , hightnudge ,    &
1636                                                  highunudge , highvnudge ,    &
1637                                                  highwnudge , highptnudge,    &
1638                                                  highqnudge
1639                IF (ierrn /= 0 )  THEN
1640                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ',  &
1641                        'higher than the maximum height in NUDING_DATA which ',&
1642                        'is ', lowheight, 'm. Interpolation on PALM ',         &
1643                        'grid is not possible.'
1644                   CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
1645                ENDIF
1646             ENDDO
1647
1648!
1649!--          Interpolation of prescribed profiles in space
1650
1651             fac = ( highheight - zu(k) ) / ( highheight - lowheight )
1652
1653             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
1654             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
1655             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
1656             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
1657             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
1658             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
1659          ENDDO
1660
1661       ENDDO rloop
1662
1663       CLOSE ( finput )
1664
1665!
1666!--    Overwrite initial profiles in case of nudging
1667       IF ( nudging )  THEN
1668          pt_init = ptnudge(:,1)
1669          u_init  = unudge(:,1)
1670          v_init  = vnudge(:,1)
1671          IF ( humidity  )  THEN ! is passive_scalar correct???
1672             q_init = qnudge(:,1)
1673          ENDIF
1674
1675          WRITE( message_string, * ) 'Initial profiles of u, v, pt and q ',    &
1676                                     'from NUDGING_DATA are used.'
1677          CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
1678       ENDIF
1679
1680
1681    END SUBROUTINE nudge_init
1682
1683!------------------------------------------------------------------------------!
1684! Description:
1685! ------------
1686!> @todo Missing subroutine description.
1687!------------------------------------------------------------------------------!
1688    SUBROUTINE calc_tnudge ( time )
1689
1690       IMPLICIT NONE
1691
1692
1693       REAL(wp) ::  dtm         !<
1694       REAL(wp) ::  dtp         !<
1695       REAL(wp) ::  time        !<
1696
1697       INTEGER(iwp) ::  k   !<
1698       INTEGER(iwp) ::  nt  !<
1699
1700       nt = 1
1701       DO WHILE ( time > timenudge(nt) )
1702         nt = nt+1
1703       ENDDO
1704       IF ( time /= timenudge(1) ) THEN
1705         nt = nt-1
1706       ENDIF
1707
1708       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1709       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1710
1711       DO  k = nzb, nzt
1712          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
1713       ENDDO
1714
1715    END SUBROUTINE calc_tnudge
1716
1717!------------------------------------------------------------------------------!
1718! Description:
1719! ------------
1720!> Call for all grid points
1721!------------------------------------------------------------------------------!
1722    SUBROUTINE nudge ( time, prog_var )
1723
1724       IMPLICIT NONE
1725
1726       CHARACTER (LEN=*) ::  prog_var  !<
1727
1728       REAL(wp) ::  tmp_tend    !<
1729       REAL(wp) ::  dtm         !<
1730       REAL(wp) ::  dtp         !<
1731       REAL(wp) ::  time        !<
1732
1733       INTEGER(iwp) ::  i  !<
1734       INTEGER(iwp) ::  j  !<
1735       INTEGER(iwp) ::  k  !<
1736       INTEGER(iwp) ::  nt  !<
1737
1738
1739       nt = 1
1740       DO WHILE ( time > timenudge(nt) )
1741         nt = nt+1
1742       ENDDO
1743       IF ( time /= timenudge(1) ) THEN
1744         nt = nt-1
1745       ENDIF
1746
1747       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1748       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1749
1750       SELECT CASE ( prog_var )
1751
1752          CASE ( 'u' )
1753
1754             DO  i = nxl, nxr
1755                DO  j = nys, nyn
1756
1757                   DO  k = nzb+1, nzt
1758
1759                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
1760                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1761
1762                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1763                                        MERGE( 1.0_wp, 0.0_wp,                 &
1764                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1765
1766                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *             &
1767                                     weight_substep(intermediate_timestep_count)
1768                   ENDDO
1769                 
1770                   sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1771 
1772                ENDDO
1773            ENDDO
1774
1775          CASE ( 'v' )
1776
1777             DO  i = nxl, nxr
1778                DO  j = nys, nyn
1779
1780                   DO  k = nzb+1, nzt
1781
1782                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
1783                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1784
1785                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1786                                        MERGE( 1.0_wp, 0.0_wp,                 &
1787                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1788
1789                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *             &
1790                                     weight_substep(intermediate_timestep_count)
1791                   ENDDO
1792                 
1793                   sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1794
1795                ENDDO
1796            ENDDO
1797
1798          CASE ( 'pt' )
1799
1800             DO  i = nxl, nxr
1801                DO  j = nys, nyn
1802
1803                   DO  k = nzb+1, nzt
1804
1805                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
1806                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1807
1808                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1809                                        MERGE( 1.0_wp, 0.0_wp,                 &
1810                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1811
1812                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *             &
1813                                     weight_substep(intermediate_timestep_count)
1814                   ENDDO
1815
1816                   sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1817
1818                ENDDO
1819            ENDDO
1820
1821          CASE ( 'q' )
1822
1823             DO  i = nxl, nxr
1824                DO  j = nys, nyn
1825
1826                   DO  k = nzb+1, nzt
1827
1828                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
1829                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1830
1831                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1832                                        MERGE( 1.0_wp, 0.0_wp,                 &
1833                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1834
1835                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *             &
1836                                     weight_substep(intermediate_timestep_count)
1837                   ENDDO
1838                 
1839                   sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1840
1841                ENDDO
1842            ENDDO
1843
1844          CASE DEFAULT
1845             message_string = 'unknown prognostic variable "' // prog_var // '"'
1846             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1847
1848       END SELECT
1849
1850    END SUBROUTINE nudge
1851
1852
1853!------------------------------------------------------------------------------!
1854! Description:
1855! ------------
1856!> Call for grid point i,j
1857!------------------------------------------------------------------------------!
1858
1859    SUBROUTINE nudge_ij( i, j, time, prog_var )
1860
1861       IMPLICIT NONE
1862
1863
1864       CHARACTER (LEN=*) ::  prog_var  !<
1865
1866       REAL(wp) ::  tmp_tend    !<
1867       REAL(wp) ::  dtm         !<
1868       REAL(wp) ::  dtp         !<
1869       REAL(wp) ::  time        !<
1870
1871       INTEGER(iwp) ::  i  !<
1872       INTEGER(iwp) ::  j  !<
1873       INTEGER(iwp) ::  k  !<
1874       INTEGER(iwp) ::  nt  !<
1875
1876
1877       nt = 1
1878       DO WHILE ( time > timenudge(nt) )
1879         nt = nt+1
1880       ENDDO
1881       IF ( time /= timenudge(1) )  THEN
1882         nt = nt-1
1883       ENDIF
1884
1885       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1886       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1887
1888       SELECT CASE ( prog_var )
1889
1890          CASE ( 'u' )
1891
1892             DO  k = nzb+1, nzt
1893
1894                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
1895                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1896
1897                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1898                                        MERGE( 1.0_wp, 0.0_wp,                 &
1899                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1900
1901                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
1902                                 * weight_substep(intermediate_timestep_count)
1903             ENDDO
1904
1905             sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1906
1907          CASE ( 'v' )
1908
1909             DO  k = nzb+1, nzt
1910
1911                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
1912                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1913
1914                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1915                                        MERGE( 1.0_wp, 0.0_wp,                 &
1916                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1917
1918                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
1919                                 * weight_substep(intermediate_timestep_count)
1920             ENDDO
1921
1922             sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1923
1924          CASE ( 'pt' )
1925
1926             DO  k = nzb+1, nzt
1927
1928                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
1929                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1930
1931                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1932                                        MERGE( 1.0_wp, 0.0_wp,                 &
1933                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1934
1935                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
1936                                 * weight_substep(intermediate_timestep_count)
1937             ENDDO
1938
1939             sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1940
1941
1942          CASE ( 'q' )
1943
1944             DO  k = nzb+1, nzt
1945
1946                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
1947                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1948
1949                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1950                                        MERGE( 1.0_wp, 0.0_wp,                 &
1951                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1952
1953                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
1954                                 * weight_substep(intermediate_timestep_count)
1955             ENDDO
1956
1957             sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1958
1959          CASE DEFAULT
1960             message_string = 'unknown prognostic variable "' // prog_var // '"'
1961             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1962
1963       END SELECT
1964
1965
1966    END SUBROUTINE nudge_ij
1967
1968
1969!------------------------------------------------------------------------------!
1970! Description:
1971! ------------
1972!> @todo Missing subroutine description.
1973!------------------------------------------------------------------------------!
1974    SUBROUTINE nudge_ref ( time )
1975
1976       IMPLICIT NONE
1977
1978       INTEGER(iwp) ::  nt                    !<
1979
1980       REAL(wp)             ::  fac           !<
1981       REAL(wp), INTENT(in) ::  time          !<
1982
1983!
1984!--    Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
1985!--    needed for correct upper boundary conditions for pt and q and in case that
1986!      large scale subsidence as well as scalar Rayleigh-damping are used
1987       nt = 1
1988       DO WHILE ( time > time_vert(nt) )
1989          nt = nt + 1
1990       ENDDO
1991       IF ( time /= time_vert(nt) )  THEN
1992        nt = nt - 1
1993       ENDIF
1994
1995       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1996
1997       pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
1998       q_init  = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
1999       u_init  = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
2000       v_init  = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
2001
2002    END SUBROUTINE nudge_ref
2003
2004
2005 END MODULE lsf_nudging_mod
Note: See TracBrowser for help on using the repository browser.