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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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