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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

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