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

Last change on this file since 3031 was 3026, checked in by schwenkel, 6 years ago

Changed the name specific humidity to mixing ratio

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