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

Last change on this file since 3301 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

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