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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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